# Epigenetics and early domestication: differences in hypothalamic DNA methylation between red junglefowl divergently selected for high or low fear of humans

Johan Bélteky, Beatrix Agnvall, Lejla Bektic, Andrey Höglund, Per Jensen, Carlos Guerrero-Bosagna

<figure><img src="image_1.png" /></figure>

## To cite this version:

Johan Bélteky, Beatrix Agnvall, Lejla Bektic, Andrey Höglund, Per Jensen, et al., Epigenetics and early domestication: differences in hypothalamic DNA methylation between red junglefowl divergently selected for high or low fear of humans. Genetics Selection Evolution, 2018, 50 (1), pp.13-10.1186/s12711-018-0384-z. http://hal-02405220

HAL Id: hal-02405220

https://hal.science/hal-02405220

Submitted on 11 Dec 2019

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Epigenetics and early domestication: differences in hypothalamic DNA methylation between red junglefowl divergently selected for high or low fear of humans

Johan Bélteky ID, Beatrix Agnvall, Lejla Bektic, Andrey Höglund, Per Jensen and Carlos Guerrero-Bosagna*

**Abstract**

**Background:** Domestication of animals leads to large phenotypic alterations within a short evolutionary time-period. Such alterations are caused by genomic variations, yet the prevalence of modified traits is higher than expected if they were caused only by classical genetics and mutations. Epigenetic mechanisms may also be important in driving domesticated phenotypes such as behavior traits. Gene expression can be modulated epigenetically by mechanisms such as DNA methylation, resulting in modifications that are not only variable and susceptible to environmental stimuli, but also sometimes transgenerationally stable. To study such mechanisms in early domestication, we used as model two selected lines of red junglefowl (ancestors of modern chickens) that were bred for either high or low fear of humans over five generations, and investigated differences in hypothalamic DNA methylation between the two populations.

**Results:** Twenty-two 1-kb windows were differentially methylated between the two selected lines at *p* < 0.05 after false discovery rate correction. The annotated functions of the genes within these windows indicated epigenetic regulation of metabolic and signaling pathways, which agrees with the changes in gene expression that were previously reported for the same tissue and animals.

**Conclusions:** Our results show that selection for an important domestication-related behavioral trait such as tameness can cause divergent epigenetic patterns within only five generations, and that these changes could have an important role in chicken domestication.

**Background**

Animal domestication is a process that occurred in an evolutionary short timespan—less than 10,000 years [1]—and substantially altered several traits such as size, behavior, coat color, and physiological and morphological features in all domesticated species. Similar phenotypic variations tend to reoccur even in widely unrelated species, which are referred to as the ‘domestication phenotype’ [1–3]. This domestication process has altered the genetic structure that underlies these traits and has created distinct patterns that differentiate domestic animals

from their wild counterparts [4–7]. Similar phenotypic changes are also seen in recent experimental domestication processes such as in the silver fox [8]. Selection for tameness, or decreased fear of humans, was undoubtedly strong during the early period of domestication and has been proposed to be the major factor that drove the domesticated phenotype, which may, accordingly, have developed partly as correlated selection responses [8–11].

Accumulating evidence shows that phenotypic variation in behavior and other traits within a species or breed can be influenced by epigenetic factors [12], in addition to genetic factors. For example, variation in personality of the great tit is statistically associated with DNA

<footnote>*Correspondence: carbo@ifm.liu.se
AVIAN Behavioural Physiology and Genomics Group, IFM Biology,
Linköping University, 581 83 Linköping, Sweden</footnote>
methylation levels at a dopamine receptor gene [13]. Similarly, in rat pups, the DNA methylation status of the hippocampal glucocorticoid receptor gene is affected by maternal behavior [14]. In domestic chickens, differences in DNA methylation are related to susceptibility to disease [15], immune responses [16], growth, and metabolism [17]. Moreover, epigenetic differences (e.g. in red blood cells) can emerge simply after individuals are subjected to different rearing conditions [18].

DNA methylation, i.e. the addition of a methyl group to cytosine in CpG dinucleotides, has the potential to affect gene expression. Although DNA methylation patterns are generally maintained after cell division, they can sometimes be modified by the action of external stimuli [14, 19]. Environmentally-altered DNA methylation patterns can be transmitted through the germ line [20] and be stable in somatic tissues over generations [21–25]. Somatic epigenetic differences, whether shaped by the environment or intentionally or unintentionally selected, can in turn affect phenotypic traits. Thus, epigenetic mechanisms could be an important factor in the rapid phenotypic changes that occur during domestication. This is supported, for example, by evidence of significant hypermethylation in purebred dogs compared to wolves [26], in domestic compared to wild worms [27], and in domesticated White Leghorn chickens compared to the ancestral red junglefowl (RJF) [28]. Interestingly, in the example of Darwin’s finches, DNA methylation in blood was shown to be a better indicator of evolutionary phenotypic diversification than overall genetic changes (copy number variations) [29].

The chicken, the world's most produced species for food, was domesticated from the RJF (*Gallus gallus*) approximately 8000 years ago and is today one of the most phenotypically diverse species of vertebrates [30–32]. With both wild and domestic chickens being available, they have been used in many comparative studies that evaluated domestication effects [33–35]. These studies have, for example, shown that domesticated chickens have variable coloration patterns, grow faster, reproduce more efficiently, and differ in their social behavior and general activity compared to RJF. In the current experiment, we attempted to experimentally recreate the early domestication process of chickens by selecting RJF bidirectionally for low or high levels of fear of humans [36]. Five generations of repeated selection generated significant phenotypic variation that affected size, behavior, and metabolism [11, 36, 37], as well as transcriptomic differences in both hypothalamic and frontal cortex tissues [38, 39]. Hence, domesticated phenotypes developed in a few generations (five) as side effects of continuous selection for high or low fear of humans.

In this study, DNA samples from the hypothalamus of the fifth generation of a lineage of RJF selected for high or low fear of humans were used to perform methylated DNA immunoprecipitation (MeDIP) sequencing, in order to identify differentially-methylated (DM) regions. Previous studies on the same populations found significant differences in brain gene expression as a consequence of selection [38]. Hence, our hypothesis was that epigenetic differences, measured as differential DNA methylation, would be observed between individuals selected for high or low fear of humans over five generations. We analyzed whether these differences mainly targeted genes and genomic regions that are relevant for tameness and stress.

# Methods

## Animals and sampling

Two unrelated populations of RJF were used to generate a parental population (P0), from which two selection lines were bred for high or low fear of humans for five generations (S5). For a detailed description of the breeding scheme, selection, and housing conditions of animals, see [36, 37]. Briefly, after an outbreeding scheme that lasted for two generations, the P0 generation consisted of approximately 70 birds. The P0 generation was divided into two groups composed of animals with high or low levels of fearfulness towards humans based on scores recorded from a standardized test of fear of humans, described in detail in [36]. Individuals were then bidirectionally selected based on scores of fear of humans, resulting in two lines, i.e. one with a high fear of humans (HFH) and one with a low fear of humans (LFH). Approximately 50 animals per selection line were hatched in each generation from 5 to 10 families, and animals from the two lines were housed and reared together upon hatching in order to standardize rearing conditions. All animals received food and water ad libitum.

Animals from generation S5 were sacrificed at the age of 48 weeks by rapid decapitation. Brains were dissected from the skull and snap-frozen in liquid nitrogen within 10 min, as reported in [38]. To access the hypothalamus, which is located in the inferior region of the brain, the whole brain was dissected into smaller parts. The region between the two optic lobes was rotated upside down to reveal a white butterfly-like pattern, with the thalamus/hypothalamus centered as a small, red protruding structure. In total, 12 S5 individuals were used for this study: six HFH and six LFH birds. Each group consisted of three males and three females, which were chosen randomly from each experimental group. These same 12 individuals had previously been used in studies that examined brain transcriptomic changes [38, 39].
## DNA isolation

DNA was extracted from hypothalamus tissue using an All-Prep RNA/DNA kit (Qiagen), following the manufacturer's instructions. In short, approximately 20 to 30 mg of tissue were homogenized with 600 µL of Buffer RTL Plus using FastPrep®-24 (MP Biomedicals). DNA was separated from RNA using AllPrep DNA spin columns, which were kept on ice while RNA was purified for gene expression analysis [38]. 350 µL of Buffer AW1 were added to each spin column, which was then centrifuged for 15 s at 17,000 g. The supernatant was discarded and a mixture of 20 µL of proteinase K (20 mg/mL) and 60 µL of Buffer AW1 was added to each sample. The samples were then incubated for 5 min at room temperature, cleaned with 350 µL of Buffer AW1, centrifuged for 15 s at full speed (17,000 g), and centrifuged again with 500 µL of Buffer AW2 for 2 min at full speed. Columns were dried by an additional centrifugation at 17,000 g for 1 min and then placed in 1.5 mL microcentrifuge tubes. For elution of DNA, 50 µL of Buffer EB were added to the spin column, followed by a 10-min incubation at room temperature and centrifugation for 1 min at 8000 g. DNA concentration and purity were measured using a NanoDrop® ND-2000c (ThermoFisher Scientific).

## MeDIP

DNA methylation analysis was performed through enrichment of the methylated fraction of the genome by immunoprecipitation with an anti-methyl-cytosine antibody (MeDIP) [40], followed by next-generation sequencing (MeDIP-seq) [41]. MeDIP-seq is a validated method for determining DNA methylation that is widely used in epigenetic research, and shows excellent performance for genome coverage compared to equivalent methods [42]. The MeDIP procedure was carried out according to a protocol that was previously optimized in our laboratory using chicken DNA [43]. From each sample, 4 µg DNA were used for the MeDIP capture. DNA was diluted in H₂O to a total volume of 8.0 µL, and sonicated at “high” setting for six 30-s intervals using a Bioruptor® Standard sonicator (Diagenode). Samples were then run on a 2% agarose gel for confirmation of fragment lengths. Samples were cleaned from excessively long fragments using a PCR purification kit (Qiagen), following the manufacturer's instructions, and sample concentrations were measured using a NanoDrop® ND-2000c (ThermoFisher Scientific). The sonicated DNA was diluted with 450 µL of TE buffer (10 mM Tris HCl, pH 7.5; 1 mM EDTA) and heat-denatured at 95 °C for 10 min, followed by cooling on ice for 5 min. To each sample, 51 µL of 10× IP buffer (100 mM NaPhosphate, pH 7; 5 M NaCl, 250 µL Triton-X 100) were added, followed by the addition of 10 µg of antibody (monoclonal mouse anti 5-methylcytosine (2 µg/µL), 5-mC, Diagenode). Then, samples were

incubated at 4 °C for 2 h on a rotating platform. Agarose beads (Protein A/G Plus-Agarose, Santa Cruz Biotechnology) were washed before use by centrifuging 50 µL of bead suspension in a 1.5 mL microcentrifuge tube for 2 min at 6000 rpm at 4 °C. The supernatant was removed and 1 mL of PBS-BSA 0.1% solution (1 mL 1× PBS + 2 µL 50 mg/ml BSA) was added. The mixture was incubated for 5 min at 4 °C on a rotating platform and then centrifuged for 2 min at 6000 rpm and 4 °C. The supernatant was subsequently removed, and the cleaning steps were repeated three more times. After cleaning, 50 µL of 1× IP buffer were added to the washed beads, and the DNA-antibody mixture was transferred to the bead mixture. The solution was incubated for 2 h at 4 °C on a rotating platform. The beads and captured DNA-antibodies were washed by centrifugation of the mixture for 2 min at 6000 rpm and 4 °C. The supernatant was removed and 1 mL of 1× IP buffer was added. The mixture was incubated for 5 min at 4 °C on a rotating platform, followed by centrifugation for 2 min at 6000 rpm and 4 °C, and removal of the supernatant. This procedure was repeated three times. To digest the beads, 210 µL of digestion buffer and 20 µL of Proteinase K (20 mg/mL) were added, followed by incubation for 2 h at 55 °C on a rotating platform. The DNA was cleaned from the bead fragments by filtering through Pierce™ Spin Cups—Paper Filter (ThermoFisher Scientific) while centrifuging at max speed for 30 s. The flow-through was collected and 3 µL of glycerin (5 mg/mL) were added. DNA was precipitated by adding 20 µL of 5 M NaCl and 750 µL ethanol, both ice cold. The solution was mixed and incubated for 30 min on ice before centrifugation at 14,000 rpm for 30 min at 4 °C. The supernatant was carefully removed and samples dried in a heating block at 50 °C for 5 min. Samples were resuspended in 30 µL of H₂O and heated on a heating block at 50 °C for 5 min before measuring DNA concentrations on a Nanodrop. The DNA samples were then used for whole-genome amplification using a WGA2 kit (Sigma-Aldrich), following the manufacturer's instructions. WGA samples were cleaned using a QIAquick PCR Purification Kit (Qiagen) following the manufacturer's instructions, and used for next-generation sequencing.

## Next-generation sequencing

The 12 samples of DNA that were extracted from hypothalamus tissue from six HRH and six LFH individuals and enriched for the methylated fraction were sequenced. Sequencing was performed on IonProton chips (Ion Torrent Systems, Inc) at the National Genomics Infrastructure (NGI), using the Ion Fragment Library kit according to the manufacturer's protocol. Data generated from the sequencing was processed using the Torrent Suite (version 5.0.2) software (ThermoFisher Scientific), which
is a complete analysis solution for Ion Torrent data that includes signal processing, base calling, trimming of low-quality reads, as well as alignment to a reference genome. After quality checks, the data were aligned to the Galgal4 chicken reference genome (International Chicken Genome Consortium). Two samples of DNA enriched for the methylated fraction were loaded per Ion PI chip. Thus, we used six chips for 12 samples. The raw sequencing data were uploaded to Array Express (http://www.ebi.ac.uk/arrayexpress/) under accession number E-MTAB-6407.

## Statistical analyses

Pairwise sample comparisons of methylated regions between selection groups, sexes or the interaction between sex and selection were performed in R (version 3.3.1) using edgeR, integrated in the Bioconductor (release 3.3) package MEDIPS [44], by dividing the genome into 1-kb windows, starting at the first position on chromosome 1, and trimming away duplicated sequences. Benjamini-Hochberg correction via R's p-adjust function was used to adjust for false discovery rate (FDR) in the multiple testing, with a cut-off of *p* < 0.1 used for significance. The cut-off was used to detect changes in a larger number of windows, allowing for a better overview of the effects of selection. In order to validate the results, random resampling was performed by randomizing samples into two random groups of six samples each, regardless of HFH or LFH selection, and running the entire MEDIPS analysis as above. The randomized groups were balanced regarding sex, to avoid sex differences in the results. Twenty replicated runs were carried out.

DM regions that covered an annotated gene region or its promoter, defined as starting 7.25 kb upstream of the transcriptional start site [28], were included in the set of regions for functional analyses. The software WebGestalt2017 (updated 1/27/2017) was used for gene ontology enrichment and KEGG pathway analysis on gene symbols obtained via Ensembl [45]. WebGestalt uses a hypergeometric distribution for significance estimation and combines it with Bonferroni-Hochberg adjustment of *p* values. All known chicken genes in the chicken genome from Ensembl were used as a base for enrichment analysis. In order to identify gene network modules that were over-represented in the DM regions, the web-based ConsensusPathDB tool was used [46] with standard settings and allowing for intermediate nodes.

## Results

### Differentially-methylated regions

A total of 430 million sequencing reads were generated, of which more than 99% were aligned to the chicken

reference genome, resulting in a 9× genome coverage per sample. Each sample generated approximately 8.5 Gb of sequence data, with an average read length of 160 bp.

A total of 990,000 windows were analyzed by edgeR, which covered 10.6 million CpGs in the chicken autosomal genome. By sorting windows by the estimated log-transformed fold change between the HFH and LFH lines, approximately 50% of all windows were found to have a negative fold change (i.e. more methylation in the LFH than HFH samples). Moreover, the distribution of CpGs in windows with positive or negative fold changes occurred in equal proportions, which indicated that there was no bias for CpG enrichment in either of the two selection groups (Table 1). Using MEDIPS on sequence counts, 51,048 significantly differentially methylated windows were detected at *p* < 0.05, with only 22 that remained significant at *p* < 0.1 adjusted by FDR (Fig. 1 and Table 2). Nine of the 22 windows that displayed significant differential methylation were hypermethylated in the LFH group. The number of CpGs within the 1-kb windows ranged from 3 to 23, yielding a CpG density of less than 3 CpG/100 bp, which is in line with previous results obtained on somatic cells when analyzing transgenerational DM regions in response to environmental challenges in rats by MeDIP [47].

The random resampling did not result in significantly differentially-methylated windows in any of the 20 replicates [see Additional file 1: Table S1].

### Sex effects

Sex-specific effects on DNA methylation, estimated when comparing combined HFH/LFH females against HFH/LFH males, were much larger than the effects of selection. Comparison of DNA methylation patterns between sexes, with line excluded as a factor, revealed

<table><caption>**Table 1 Summary statistics of sequenced methylated regions in the selection lines**</caption><thead><tr><th></th><th>Total</th><th>Positive logFC</th><th>Negative logFC</th></tr></thead><tbody><tr><td>Windows</td><td>990,461</td><td>489,962</td><td>500,499</td></tr><tr><td>CpGs in windows</td><td>10,610,935</td><td>5,174,029</td><td>5,427,885</td></tr><tr><td>Significant windows</td><td>22</td><td>13</td><td>9</td></tr><tr><td>CpGs in significantly DM windows</td><td>187</td><td>132</td><td>55</td></tr><tr><td>Significant windows males</td><td>51</td><td>21</td><td>30</td></tr><tr><td>Significant windows females</td><td>66</td><td>30</td><td>36</td></tr></tbody></table>

Total number of windows and number of CpG dinucleotides within the windows, and differences between selection objectives for both groups and within sex. A positive log-transformed fold change (logFC) indicates windows with less methylation in LFH than in HFH animals, and vice versa for negative logFC for which LFH animals have more methylated strands than HFH animals. Significant windows are those with FDR-corrected *p* < 0.1
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Manhattan plot of differentially-methylated (DM) windows between selection lines. All 990,000 windows were plotted with genomic location on the X-axis and negative log10 *p* values on the Y-axis. The red horizontal line indicates the threshold for significantly DM windows at *p* < 0.1 after FDR correction. Microchromosome labels were filtered out for readability</figcaption></figure>

many differences, with over 2500 DM windows on the autosomes and almost 2000 on the sex chromosomes W and Z [see Additional file 2: Figure S1]. A more detailed analysis of the DM regions indicated that 79% of the windows were hypermethylated in males.

Successive selection of HFH and LFH birds seemed to result in large sex-specific effects in DNA methylation differences, with 66 windows being DM between HFH and LFH females (FDR adjusted *p* < 0.1) [see Additional file 3: Table S2]. Moreover, 51 DM windows were identified between HFH and LFH males (FDR adjusted *p* < 0.1). Comparison of DM windows identified within males and within females showed overlaps at only four windows, three of which were located on chromosome 7 and covered the expressed sequence tag (EST) Gga.15462, which was previously detected in testis mRNA in RJF and White Leghorn [48]. The methylation patterns for

<table><caption>**Table 2** Differentially-methylated windows between the selection lines</caption><thead><tr><th>Chr</th><th>Window start (bp)</th><th>CpGs</th><th>logFC</th><th>Adjusted p value</th><th>Gene symbol</th><th>Gene description</th></tr></thead><tbody><tr><td>1</td><td>121,257,001</td><td>4</td><td>-1.33</td><td>9.52e-02</td><td>AP1S2</td><td>Adaptor related protein complex 1 sigma 2 subunit [Source:HGNC Symbol;Acc:HGNC:560]</td></tr><tr><td>1</td><td>142,947,001</td><td>4</td><td>-0.93</td><td>5.87e-02</td><td>NALCN</td><td>Sodium leak channel, non-selective [Source:CGNC Symbol;Acc:12660]</td></tr><tr><td>1</td><td>174,158,001</td><td>7</td><td>-0.90</td><td>9.52e-02</td><td></td><td></td></tr><tr><td>1</td><td>188,959,001</td><td>7</td><td>1.40</td><td>5.14e-02</td><td>DLG2</td><td>Discs, large homolog 2 (Drosophila) [Source:CGNC Symbol;Acc:12965]</td></tr><tr><td>1</td><td>192,474,001</td><td>5</td><td>1.35</td><td>3.27e-07</td><td>CAPN5</td><td>Calpain 5 [Source:CGNC Symbol;Acc:479]</td></tr><tr><td>2</td><td>73,558,001</td><td>3</td><td>-2.75</td><td>8.12e-02</td><td></td><td></td></tr><tr><td>2</td><td>85,911,001</td><td>12</td><td>0.80</td><td>9.52e-02</td><td></td><td></td></tr><tr><td>3</td><td>9,149,001</td><td>7</td><td>1.11</td><td>3.22e-02</td><td>EHBP1</td><td>EH domain binding protein 1 [Source:CGNC Symbol;Acc:6763]</td></tr><tr><td>3</td><td>53,750,001</td><td>7</td><td>1.31</td><td>2.73e-02</td><td>ARFGEF3</td><td>ARFGEF family member 3 [Source:CGNC Symbol;Acc:10361]</td></tr><tr><td>3</td><td>106,780,001</td><td>5</td><td>1.70</td><td>4.83e-03</td><td>MTMR9</td><td>Myotubularin related protein 9 [Source:CGNC Symbol;Acc:12476]</td></tr><tr><td>4</td><td>27,239,001</td><td>3</td><td>-0.98</td><td>9.24e-02</td><td></td><td></td></tr><tr><td>5</td><td>30,875,001</td><td>3</td><td>-0.80</td><td>3.22e-02</td><td></td><td></td></tr><tr><td>5</td><td>57,981,001</td><td>7</td><td>1.12</td><td>7.28e-02</td><td>C5h14orf166</td><td>Gallus gallus chromosome 5 open reading frame, mRNA. [Source:RefSeq mRNA;Acc:NM_205369]</td></tr><tr><td>5</td><td>58,733,001</td><td>12</td><td>0.99</td><td>7.79e-02</td><td></td><td></td></tr><tr><td>11</td><td>16,608,001</td><td>16</td><td>-1.29</td><td>5.91e-02</td><td></td><td></td></tr><tr><td>14</td><td>5,936,001</td><td>11</td><td>-0.83</td><td>8.12e-02</td><td></td><td></td></tr><tr><td>15</td><td>6,629,001</td><td>22</td><td>1.13</td><td>3.22e-02</td><td>USP30</td><td>Ubiquitin specific peptidase 30 [Source:CGNC Symbol;Acc:50370]</td></tr><tr><td>17</td><td>1,582,001</td><td>8</td><td>1.87</td><td>8.12e-02</td><td>MRPL41</td><td>Mitochondrial ribosomal protein L41 [Source:CGNC Symbol;Acc:6639]</td></tr><tr><td>17</td><td>1,582,001</td><td>8</td><td>1.87</td><td>8.12e-02</td><td>PNPLA7</td><td>Patatin-like phospholipase domain containing 7 [Source:CGNC Symbol;Acc:6632]</td></tr><tr><td>18</td><td>10,927,001</td><td>23</td><td>0.89</td><td>5.14e-02</td><td>ARMC7</td><td>Armadillo repeat containing 7 [Source:CGNC Symbol;Acc:6023]</td></tr><tr><td>19</td><td>58,001</td><td>3</td><td>-0.95</td><td>1.01e-02</td><td>DDX25</td><td>DEAD-box helicase 25 [Source:HGNC Symbol;Acc:HGNC:18698]</td></tr><tr><td>20</td><td>11,544,001</td><td>11</td><td>1.15</td><td>3.22e-02</td><td></td><td></td></tr><tr><td>23</td><td>3,692,001</td><td>5</td><td>1.38</td><td>7.28e-02</td><td>GNL2</td><td>Guanine nucleotide binding protein-like 2 (nucleolar) [Source:CGNC Symbol;Acc:1398]</td></tr></tbody></table>

Comparison between HFH and LFH red junglefowl individuals resulted in 22 differentially-methylated 1-kb windows at FDR-corrected *p* < 0.1. Of these, 13 aligned within genes or promoter regions of genes. The log fold changes (logFC) are based on the comparison between HFH and LFH birds, with HFH birds set as reference in the MEDIPS analysis
this EST region differed between selection lines and between sexes, i.e. it was highly methylated in HFH males and in LFH females. The fourth window was located on chromosome 1 in an intron of the *CAPN5* gene and was hypermethylated in HFH animals of both sexes.

## Functional annotation

Genes from the 22 windows that showed significant DM between HFH and LFH animals were extracted via Ensembl and used for gene ontology (GO) and pathway (WebGestalt) analyses. Comparison of genes associated with the DM windows between the two lines showed enrichment for a number of terms, but not after FDR correction (Table 3). The enriched terms were related to mitochondria, ion channel activity, ribosomes, and regulation of membrane potential. Due to the small size of the dataset, KEGG pathway enrichment analysis yielded no results. GO analysis of sex-specific line effects showed little overlap of enrichment in males compared to females but yielded terms related to DNA replication, the GABA

<table><caption>**Table 3** Enrichment results for gene ontology (GO) analysis of differentially-methylated regions between selection lines across and within sexes</caption><thead><tr><th>GO ID</th><th>GO term</th></tr></thead><tbody><tr><td colspan="2">Across sexes</td></tr><tr><td>GO:0005261</td><td>Cation channel activity</td></tr><tr><td>GO:0003723</td><td>RNA binding</td></tr><tr><td>GO:0003924</td><td>GTPase activity</td></tr><tr><td>GO:0005759</td><td>Mitochondrial matrix</td></tr><tr><td>GO:0006412</td><td>Translation</td></tr><tr><td>GO:0006816</td><td>Calcium ion transport</td></tr><tr><td>GO:0042254</td><td>Ribosome biogenesis</td></tr><tr><td>GO:0042391</td><td>Regulation of membrane potential</td></tr><tr><td colspan="2">Female</td></tr><tr><td>GO:0044456</td><td>Synapse part</td></tr><tr><td>GO:0043197</td><td>Dendric spine</td></tr><tr><td>GO:0098878</td><td>Neurotransmitter complex</td></tr><tr><td>GO:0022803</td><td>Passive transmembrane transporter activity</td></tr><tr><td>GO:0016597</td><td>Amino acid binding</td></tr><tr><td>GO:0008066</td><td>Glutamate receptor activity</td></tr><tr><td colspan="2">Male</td></tr><tr><td>GO:1902710</td><td>GABA receptor complex</td></tr><tr><td>GO:0045211</td><td>Postsynaptic membrane</td></tr><tr><td>GO:0005657</td><td>Replication fork</td></tr><tr><td>GO:0034776</td><td>Response to histamine</td></tr><tr><td>GO:0006312</td><td>Mitotic recombination</td></tr><tr><td>GO:0004520</td><td>Endodeoxyribonuclease activity</td></tr><tr><td>GO:0005254</td><td>Chloride channel activity</td></tr><tr><td>GO:0008094</td><td>DNA-dependent ATPase activity</td></tr></tbody></table>

Annotations for differentially-methylated windows were used to search for enrichment of terms in the groups, both for overall effects of selection objective and sex-specific changes

receptor complex and chloride channel activity in males and terms related to transmembrane transporter activity, the synapse part, and neurotransmitter complex in females.

## Transcriptional effect

We examined the overlap between the DM windows identified here and the genomic regions that included genes that were previously reported as differentially expressed (DE) in the hypothalamus between the same selection groups and generation [38]. Genomic regions located 7.25 kb upstream from DE transcripts were evaluated in order to include promoter regions [28]. No direct overlaps were detected between DM and DE regions, even when including DM regions with a fold change higher than 1.5 and *p* < 0.05 before FDR correction. ConsensusPathDB analysis of these regions revealed that the DM genes *MRPL41* and *ARMC7* (Table 2) were significantly connected with the DE genes *MAEA* and *ANKRD1* that were detected by Bélteky et al. [38] through one intermediate node. DM windows within males and within females were also compared to previously reported DE genes, but neither sex showed any overlap between expression and methylation changes [38].

## Discussion

In this experiment, we studied the effects of bidirectional selection for tameness on the DNA methylome of the hypothalamus in RJF. After five generations of selection for high and low fear of humans, 22 genomic regions with differential DNA methylation in the hypothalamus were detected. The DM regions targeted several genes with functions that may be related to previously reported phenotypic selection responses [11, 36], such as metabolism and signaling. However, the number of significant windows before correcting for multiple-testing was much larger than 22 windows, with over 51,000 differently methylated windows observed at *p* > 0.05. This large reduction in significant windows after correcting for FDR may indicate that extensive variation exists within the groups for methylation patterns in regions that are not strongly affected by the selection process. Such naturally occurring inter-individual variation has already been reported in humans [49].

Among the 22 DM windows, none overlapped with DE genes reported for the same tissue of birds from the same populations, which agrees with other studies on RJF and domesticated White Leghorn in which no overlap between DM promoters and gene expression of the same genes was observed. This suggests that other epigenetic mechanisms are involved in regulating gene expression [28], or that DNA methylation changes in regions outside
of promoter and gene regions, i.e. in intergenic regions, could affect transcription in a trans-acting fashion [50, 51], or that the DM regions affect genes that are not activated by any pathways at the life stage studied here. Furthermore, gene expression may be affected primarily by methylation changes of single CpG sites within a CpG-dense region, which may not be detected by MeDIP-seq [52].

Although we observed no direct overlap between DM and DE genes, GO analysis of genes included in DM windows revealed functions similar to those of the DE genes in the same group of animals. Interestingly, some DM regions were related to behavioral and metabolic pathways. Further research should confirm whether and how these pathways are affected by selection for tameness. However, based on the genes targeted by the DM regions between HFH and LFH animals, we can derive some hints. For example, *calpain 5 (CAPN5)*, which encodes a cytosolic cysteine protease that acts on signaling-related molecules and is involved in cell differentiation and proliferation [53], was reported to interact with nuclear receptors and thus to impact metabolism [54]. Mutations in the gene *sodium leak channel non-selective (NALCN)*, which is involved in the control of neuronal excitability, can lead to speech impairment and intellectual disability in humans [55, 56]. The gene *DLG2*, which encodes a scaffold protein belonging to the membrane-associated guanylate kinase (MAGUK) family and is active at post-synaptic sites, is associated with neurodevelopmental disorders, schizophrenia and cognition in humans [57, 58]. Downregulation of the mitochondrial gene *deubiquitinase USP30* is reported to enhance degradation of damage mitochondria in human neurons, which is beneficial during neurological disorders [59]. The fact that several of the genes related to DM regions are associated with neuronal functions suggests that molecular mechanisms involved in the behavioral changes have emerged between HFH and LFH birds [60].

We compared the DM regions detected here with DM probes that were previously detected in the hypothalamus by Nätt et al. [28] between RJF and the domesticated White Leghorn, in order to test whether the same genes with a modified methylation status were found in the current experiment and during the domestication of White Leghorn. We did not identify any common genes between these two selection processes, which could be explained by the fact that the driver of domestication of White Leghorn was mainly selection for egg size, while in our selection lines the main driver was tameness. This suggests that different selection pressures generate distinctive sets of epigenetic changes, which in turn are related to specific phenotypic traits. This concept can also apply to the transcriptome, since we observed

no concordance between previously published DE genes from the same birds as those used in this study and DE genes observed in the hypothalamus of RJF and White Leghorn [28, 38].

A possible explanation for the lack of overlap between DM regions in the current and previous studies of the chicken hypothalamus is that the observed differences may be related to differential genetic drift in HFH and LFH populations [28]. However, previous studies on gene expression differences in the same populations, along with an unselected line of RJF, indicated that the expression changes that occurred between the parental generation and the fifth selected generations resulted from the imposed selection and not from genetic drift [38, 39]. These results, together with the random resampling performed in the current experiment, strengthen the argument that the epigenetic changes observed are a consequence of the artificial selection imposed within a rather short period of time (5 generations), and not of genetic drift.

In addition to the DM regions detected between the HFH and LFH lines in this experiment, we found numerous sex-specific changes, which are consistent with previously reported results [61]. However, in our study, differences were larger than observed previously, because it covered the entire genome, in contrast to only the promoter regions in a previous study that analyzed DNA methylation using promoter arrays [61]. The largest sex-specific DM regions were associated to the genes *ZFR* on chromosome 1 and *MHM* on chromosome Z, which were both highly methylated in males. This large number of sex-specific differences in DNA-methylation suggests that the two sexes respond differently to selection pressures related to behavioral traits, in agreement with previous findings [60], but raising some interesting questions about underlying mechanisms. An interesting finding was that many of these sex-specific DNA methylation changes are on autosomes, which suggests that DNA methylation of autosomes plays a more important role than expected in gender-specific characteristics in the brain of vertebrates.

## Conclusions

We detected 22 DM regions by comparing hypothalamic DNA from RJF selected during five generations for high versus low fear of humans. Functional annotation of the genes associated with these DM regions showed that they are related to, for example, cellular metabolism and neural signaling, similar to what was previously reported in terms of gene expression differences for the same animals. Our results suggest that bidirectional selection for tameness involves epigenetic factors that can even differ in a sex-specific manner. Observation of divergent
DNA methylation patterns in the hypothalamus after only five generations of artificial selection highlights the importance of epigenetic mechanisms, in addition to genetic composition, in evolutionary phenotypic variation that emerges in response to selection pressures. Future research should delve into the molecular mechanism involved in the emergence of somatic epigenetic differences during selection. Two options could explain this. One is that epigenetic differences are linked to genetic differences emerging during selection, as suggested by Verhulst et al. [13]. Another possibility is that selection on phenotypes would concomitantly select specific germ line epigenomes. Different germ line epigenomes in divergent selection lines would then influence somatic epigenomes later on in the ontogeny of each individual. Interestingly, we detected several sex-specific epigenetic changes on the autosomes, which raises the question about whether epigenetic differences have a role in gender-specific behavioral responses that are unrelated to sex chromosomes.

## Additional files

**Additional file 1: Table S1.** Randomized resampling to validate results, with distribution of HFH and LFH animals in each randomization. Three females and three males were selected at random from the pool of samples, while also allowing for resampling of the same individual more than once. Two groups were generated and analyzed using MEDIPS and the same parameters as those set for the treatment groups. The procedure was repeated 20 times. HFH—High Fear of Humans, LFH—Low Fear of Humans.

**Additional file 2: Figure S1.** Differentially-methylated (DM) windows between males and females. Differences in methylation between sexes for each window in the genome were visualized via a Manhattan plot. The red horizontal line indicates the threshold for significantly DM windows at *p* < 0.1 after FDR correction. Microchromosome labels have been filtered out for readability.

**Additional file 3: Table S2.** Sex-specific differentially-methylated (DM) windows. The MEDIPS analysis yielded 66 DM windows in males and 66 windows in females at *p* < 0.1. The log fold changes (logFC) were obtained by using highly fearful birds as reference in the MEDIPS analysis.

## Authors' contributions

All authors contributed extensively to the work presented in this paper. JB, BA and PJ designed the experiment. BA supervised breeding of the selection lines. JB, LB, AH and CGB analyzed and interpreted the sequencing data. JB administered the experiment and wrote the first version of the manuscript. JB and CGB supervised the final version. All authors read and approved the final manuscript.

## Acknowledgements

The authors would like to acknowledge support of the National Genomics Infrastructure (NGI)/Uppsala Genome Center and UPPMAX for providing assistance in massive parallel sequencing and computational infrastructure. Work performed at NGI/Uppsala Genome Center has been funded by RFI/VR and Science for Life Laboratory, Sweden. We also extend thanks to Petros Batakis for sample preparations, and special thanks to Markus Jöngren for help in developing the breeding programme.

## Competing interests

The authors declare that they have no competing interests.

## Availability of data and materials

The raw sequencing data have been uploaded to Array Express (http://www.ebi.ac.uk/arrayexpress/) under accession number E-MTAB-6407.

## Consent for publication

Not applicable.

## Ethics approval and consent to participate

The experiments reported in this paper were carried out in accordance with ethical guidelines, and were approved by the Linköping Animal Ethics Committee, license no 122-10.

## Funding

The work was supported by grants from the research council Formas, Vetenskapsrådet, and the European Research Council, ERC (Grant No. 322206 GENEWELL).

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Received: 15 August 2017 Accepted: 9 March 2018

Published online: 02 April 2018

## References

1. Zeder MA. The domestication of animals. J Anthropol Res. 2012;68:161–90.

2. Price EO. Animal domestication and behavior. Wallingford: CABI Publishing; 2002.

3. Sánchez-Villagra MR, Geiger M, Schneider RA. The taming of the neural crest: a developmental perspective on the origins of morphological covariation in domesticated mammals. R Soc Open Sci. 2016;3:160107.

4. Rubin CJ, Zody MC, Eriksson J, Meadows JR, Sherwood E, Webster MT, et al. Whole-genome resequencing reveals loci under selection during chicken domestication. Nature. 2010;464:587–91.

5. Andersson L, Georges M. Domestic-animal genomics: deciphering the genetics of complex traits. Nat Rev Genet. 2004;5:202–12.

6. Heyne HO, Lautenschlager S, Nelson R, Besnier F, Rotival M, Cagan A, et al. Genetic influences on brain gene expression in rats selected for tameness and aggression. Genetics. 2014;198:1277–90.

7. Kukekova AV, Johnson JL, Teiling C, Li L, Oskina IN, Kharlamova AV, et al. Sequence comparison of prefrontal cortical brain transcriptome from a tame and an aggressive silver fox (*Vulpes vulpes*). BMC Genomics. 2011;12:482.

8. Trut L, Oskina I, Kharlamova A. Animal evolution during domestication: the domesticated fox as a model. BioEssays. 2009;31:349–60.

9. The Belyaev DK, Wilhelmine E. Key 1978 invitational lecture. Destabilizing selection as a factor in domestication. J Hered. 1979;70:301–8.

10. Albert FW, Shchepina O, Winter C, Rompler H, Teupser D, Palme R, et al. Phenotypic differences in behavior, physiology and neurochemistry between rats selected for tameness and for defensive aggression towards humans. Horm Behav. 2008;53:413–21.

11. Agnvall B, Katajamaa R, Altimiras J, Jensen P. Is domestication driven by reduced fear of humans? Boldness, metabolism and serotonin levels in divergently selected red junglefowl (*Gallus gallus*). Biol Lett. 2015;11:20150509.

12. Jensen P. Adding 'epi-' to behaviour genetics: implications for animal domestication. J Exp Biol. 2015;218:32–40.

13. Verhulst EC, Mateman AC, Zwier MV, Caro SP, Verhoeven KJF, Van Oers K, et al. Evidence from pyrosequencing indicates that natural variation in animal personality is associated with *DRD4* DNA methylation. Mol Ecol. 2016;25:1801–11.
14. Weaver IC, Cervoni N, Champagne FA, D'Alessio AC, Sharma S, Seckl JR, et al. Epigenetic programming by maternal behavior. Nat Neurosci. 2004;7:847–54.

15. Tian F, Zhan F, VanderKraats ND, Hiken JF, Edwards JR, Zhang H, et al. DNMT gene expression and methylome in Marek's disease resistant and susceptible chickens prior to and following infection by MDV. Epigenetics. 2013;8:431–44.

16. Berghof TVL, Parmentier HK, Lammers A. Transgenerational epigenetic effects on innate immunity in broilers: an underestimated field to be explored? Poult Sci. 2013;92:2904–13.

17. Hu Y, Xu H, Li Z, Zheng X, Jia X, Nie Q, et al. Comparison of the genome-wide DNA methylation profiles between fast-growing and slow-growing broilers. PLoS One. 2013;8:e56411.

18. Pertille F, Brantsaeter M, Nordgreen J, Coutinho LL, Janczak AM, Jensen P, et al. DNA methylation profiles in red blood cells of adult hens correlate to their rearing conditions. J Exp Biol. 2017;220:3579–87.

19. Raynal NJ, Si J, Taby RF, Gharibyan V, Ahmed S, Jelinek J, et al. DNA methylation does not stably lock gene expression but instead serves as a molecular mark for gene silencing memory. Cancer Res. 2012;72:1170–81.

20. Guerrero-Bosagna C, Settles M, Lucker B, Skinner MK. Epigenetic transgenerational actions of vinclozolin on promoter regions of the sperm epigenome. PLoS One. 2010;5:e13100.

21. Lindqvist C, Janczak AM, Natt D, Baranowska I, Lindqvist N, Wichman A, et al. Transmission of stress-induced learning impairment and associated brain gene expression from parents to offspring in chickens. PLoS One. 2007;2:e364.

22. Goerlich VC, Natt D, Elfwing M, Macdonald B, Jensen P. Transgenerational effects of early experience on behavioral, hormonal and gene expression responses to acute stress in the precocial chicken. Horm Behav. 2012;61:711–8.

23. Gapp K, Soldado-Magraner S, Alvarez-Sánchez M, Bohacek J, Vernaz G, Shu H, et al. Early life stress in fathers improves behavioural flexibility in their offspring. Nat Commun. 2014;5:5466.

24. Franklin TB, Russig H, Weiss IC, Gräff J, Linder N, Michalon A, et al. Epigenetic transmission of the impact of early stress across generations. Biol Psychiatry. 2010;68:408–15.

25. Nätt D, Lindqvist N, Stranneheim H, Lundeberg J, Torjesen PA, Jensen P. Inheritance of acquired behaviour adaptations and brain gene expression in chickens. PLoS One. 2009;4:e6405.

26. Janowitz Koch I, Clark MM, Thompson MJ, Deere-Machemer KA, Wang J, Duarte L, et al. The concerted impact of domestication and transposon insertions on methylation patterns between dogs and grey wolves. Mol Ecol. 2016;25:1838–55.

27. Xiang H, Li X, Dai F, Xu X, Tan A, Chen L, et al. Comparative methylomics between domesticated and wild silkworms implies possible epigenetic influences on silkworm domestication. BMC Genomics. 2013;14:646.

28. Nätt D, Rubin CJ, Wright D, Johnsson M, Bélteky J, Andersson L, et al. Heritable genome-wide variation of gene expression and promoter methylation between wild and domesticated chickens. BMC Genomics. 2012;13:59.

29. Skinner MK, Gurerrero-Bosagna C, Haque MM, Nilsson EE, Koop JAH, Knutie SA, et al. Epigenetics and the evolution of Darwin's finches. Genome Biol Evol. 2014;6:1972–89.

30. Niu D, Fu Y, Luo J, Ruan H, Yu XP, Chen G, et al. The origin and genetic diversity of Chinese native chicken breeds. Biochem Genet. 2002;40:163–74.

31. Zuidhof MJ, Schneider BL, Carney VL, Korver DR, Robinson FE. Growth, efficiency, and yield of commercial broilers from 1957, 1978, and 2005. Poult Sci. 2014;93:2970–82.

32. Nicol CJ. The behavioural biology of cickens. Wallingford: CABI; 2015.

33. Schutz KE, Jensen P. Effects of resource allocation on behavioural strategies: a comparison of red junglefowl (*Gallus gallus*) and two domesticated breeds of poultry. Ethology. 2001;107:753–65.

34. Vaisanen J, Jensen P. Social versus exploration and foraging motivation in young red junglefowl (*Gallus gallus*) and White Leghorn layers. Appl Anim Behav Sci. 2003;84:139–58.

35. Lindgren I, Altimiras J. Sensitivity of organ growth to chronically low oxygen levels during incubation in red junglefowl and domesticated chicken breeds. Poult Sci. 2011;90:126–35.

36. Agnvall B, Jöngren M, Strandberg E, Jensen P. Heritability and genetic correlations of fear-related behaviour in red junglefowl—possible implications for early domestication. PLoS One. 2012;7:e35162.

37. Agnvall B, Ali A, Olby S, Jensen P. Red Junglefowl (*Gallus gallus*) selected for low fear of humans are larger, more dominant and produce larger offspring. Animal. 2014;8:1498–505.

38. Bélteky J, Agnvall B, Johnsson M, Wright D, Jensen P. Domestication and tameness: brain gene expression in red junglefowl selected for less fear of humans suggests effects on reproduction and immunology. R Soc Open Sci. 2016;3:160033.

39. Bélteky J, Agnvall B, Jensen P. Gene expression of behaviorally relevant genes in the cerebral hemisphere changes after selection for tameness in Red Junglefowl. PLoS One. 2017;12:e0177004.

40. Weber M, Davies JJ, Wittig D, Oakeley EJ, Haase M, Lam WL, et al. Chromosome-wide and promoter-specific analyses identify sites of differential DNA methylation in normal and transformed human cells. Nat Genet. 2005;37:853–62.

41. Taiwo O, Wilson GA, Morris T, Seisenberger S, Reik W, Pearce D, et al. Methylome analysis using MeDIP-seq with low DNA concentrations. Nat Protoc. 2012;7:617–36.

42. Bock C, Tomazou EM, Brinkman AB, Muller F, Simmer F, Gu H, et al. Quantitative comparison of genome-wide DNA methylation mapping technologies. Nat Biotechnol. 2010;28:1106–14.

43. Guerrero-Bosagna C, Jensen P. Optimized method for methylated DNA immuno-precipitation. MethodsX. 2015;2:432–9.

44. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.

45. McLaren W, Pritchard B, Rios D, Chen Y, Flicek P, Cunningham F. Deriving the consequences of genomic variants with the Ensembl API and SNP effect predictor. Bioinformatics. 2010;26:2069–70.

46. Kamburov A, Pentchev K, Galicka H, Wierling C, Lehrach H, Herwig R. ConsensusPathDB: toward a more complete picture of cell biology. Nucleic Acids Res. 2011;39:D712–7.

47. Skinner MK, Guerrero-Bosagna C. Role of CpG deserts in the epigenetic transgenerational inheritance of differential DNA methylation regions. BMC Genomics. 2014;15:692.

48. Savolainen P, Fitzsimmons C, Arvestad L, Andersson L, Lundeberg J. ESTs from brain and testis of White Leghorn and red junglefowl: annotation, bioinformatic classification of unknown transcripts and analysis of expression levels. Cytogenet Genome Res. 2005;111:79–87.

49. Rakyan VK, Hildmann T, Novik KL, Lewin J, Tost J, Cox AV, et al. DNA methylation profiling of the human major histocompatibility complex: a pilot study for the human epigenome project. PLoS Biol. 2004;2:e405.

50. Maunakea AK, Nagarajan RP, Bilenky M, Ballinger TJ, Souza C, Fouse SD, et al. Conserved role of intragenic DNA methylation in regulating alternative promoters. Nature. 2010;466:253–7.

51. Smith ZD, Chan MM, Mikkelsen TS, Gu H, Gnarke A, Regev A, et al. A unique regulatory phase of DNA methylation in the early mammalian embryo. Nature. 2012;484:339–44.

52. Jeong HM, Lee S, Chae H, Kim R, Kwon MJ, Oh E, et al. Efficiency of methylated DNA immunoprecipitation bisulphite sequencing for whole-genome DNA methylation analysis. Epigenomics. 2016;8:1061–77.

53. Sato K, Kawashima S. Calpain function in the modulation of signal transduction molecules. Biol Chem. 2001;382:743–51.

54. Sáez ME, Grilo A, Morón FJ, Manzano L, Martínez-Larrad MT, González-Pérez A, et al. Interaction between calpain 5, peroxisome proliferator-activated receptor-gamma and peroxisome proliferator-activated receptor-delta genes: a polygenic approach to obesity. Cardiovasc Diabetol. 2008;7:23.

55. Al-Sayed Moeenaldeen D, Al-Zaidan H, Albakheet A, Hakami H, Kenana R, Al-Yafee Y, et al. Mutations in NALCN cause an autosomal-recessive syndrome with severe hypotonia, speech impairment, and cognitive delay. Am J Hum Genet. 2013;93:721–6.

56. Köroğlu Ç, Seven M, Tolun A. Recessive truncating NALCN mutation in infantile neuroaxonal dystrophy with facial dysmorphism. J Med Genet. 2013;50:515–20.

57. Ingason A, Giegling I, Hartmann AM, Genius J, Konte B, Friedl M, et al. Expression analysis in a rat psychosis model identifies novel candidate genes validated in a large case-control sample of schizophrenia. Transl Psychiatry. 2015;5:e656.
58. Nithianantharajah J, Komiyama NH, McKechanie A, Johnstone M, Blackwood DH, Clair DS, et al. Synaptic scaffold evolution generated components of vertebrate cognitive complexity. Nat Neurosci. 2013;16:16–24.

59. Bingol B, Tea JS, Phu L, Reichelt M, Bakalarski CE, Song Q, et al. The mitochondrial deubiquitinase USP30 opposes parkin-mediated mitophagy. Nature. 2014;510:370–5.

60. Agnvall B, Jensen P. Effects of divergent selection for fear of humans on behaviour in Red Junglefowl. PLoS One. 2016;11:e0166075.

61. Nätt D, Agnvall B, Jensen P. Large sex differences in chicken behavior and brain gene expression coincide with few differences in promoter DNA-methylation. PLoS One. 2014;9:e96376.

Submit your next manuscript to BioMed Central and we will help you at every step:

* We accept pre-submission inquiries

* Our selector tool helps you to find the most relevant journal

* We provide round the clock customer support

* Convenient online submission

* Thorough peer review

* Inclusion in PubMed and all major indexing services

* Maximum visibility for your research

Submit your manuscript at
www.biomedcentral.com/submit

( ) BioMed Central