Stressful conditions are common in the environment where production animals are reared. Stress in animals is usually determined by the levels of stress-related hormones. A big challenge, however, is in determining the history of exposure of an organism to stress, because the release of stress hormones can show an acute (and recent) but not a sustained exposure to stress. Epigenetic tools provide an alternative option to evaluate past exposure to long-term stress. Chickens provide a unique model to study stress effects in the epigenome of red blood cells (RBCs), a cell type of easy access and nucleated in birds. The present study investigated whether two different rearing conditions in chickens can be identified by looking at DNA methylation patterns in their RBCs later in life. These conditions were rearing in open aviaries versus in cages, which are likely to differ regarding the amount of stress they generate. Our comparison revealed 115 genomic windows with significant changes in RBC DNA methylation between experimental groups, which were located around 53 genes and within 22 intronic regions. Our results set the ground for future detection of long-term stress in live production animals by measuring DNA methylation in a cell type of easy accessibility.
Stress in production animals generated by unsustainable production methods is a frequent issue of concern. Besides the ethical issue of inducing unnecessary stress in animals, detrimental practices in the animal production industry have consequences from a human health perspective (Rostagno, 2009). The environment where production animals are reared influences not only their later health and wellbeing but also the quality of the food originating from them (Broom, 2010). Stressful conditions to which production animals can be subjected include extreme illumination patterns (Morgan and Tromborg, 2006; Olanrewaju et al., 2006), social isolation or crowding (Goerlich et al., 2012), food restriction (Morgan and Tromborg, 2006; Savory and Lariviere, 2000), too high or too low temperature, restriction of movement, barren environments, and lack of appropriate substrates for foraging, exploration and manipulation (Morgan and Tromborg, 2006).
Stress in animals is associated with a cascade of hormonal responses (Henry, 1992). The primary physiological stress response observed is an increase in hypothalamic-pituitary–adrenal (HPA) axis activity, which results in elevated levels of glucocorticoids (Fallahsharoudi et al., 2015). Initially, increases in testosterone levels related to increased anxiety are observed (Henry, 1992). Subsequently, decreases in the noradrenaline/adrenaline ratio are observed, concomitant with increases in adrenaline, prolactin and fatty acids (Henry, 1992). In conditions of further distress, adrenocorticotropic hormone and cortisol levels will increase (Henry, 1992). Because of this plethora of hormonal changes generated by stressful conditions, stress in animals is usually determined by the levels of stress-related hormones such as cortisol and adrenaline (Ishibashi et al., 2013; Muller et al., 2013). A big challenge, however, is in determining the history of the exposure of an organism to stress, given that the release of stress hormones can show an acute (and recent) but not a sustained exposure to stressful conditions (Henry, 1992).
An alternative option to the use of hormonal measurements to evaluate past exposure to long-term stress could be to utilize epigenetic tools instead. Epigenetics involves studying how environmental exposures affect gene regulation during the lifetime of organisms. Epigenetic changes are defined as accessory chemical modifications of the DNA that regulate gene expression and are mitotically stable (Skinner et al., 2010). These modifications include DNA methylation or hydroxymethylation of nucleotides, chemical modification of histones, interaction of DNA with small RNAs, or states of chromatin condensation (Denham et al., 2014; Feil and Fraga, 2011; Teperek-Tkacz et al., 2011). Altering epigenetic states can lead to distinguishable phenotypic consequences such as changes in coat colour (Dolinoy et al., 2007) or an increased susceptibility to diseases (Guerrero-Bosagna and Skinner, 2012). A variety of model organisms have been used in epigenetic research, including laboratory rodents (Dolinoy et al., 2007; Guerrero-Bosagna et al., 2008; Susiarjo et al., 2013), flies (Seong et al., 2011), honey bees (Dickman et al., 2013; Lyko et al., 2010), plants (Cubas et al., 1999; Manning et al., 2006) and yeast (Zhang et al., 2013). However, in spite of the importance of epigenetic mechanisms in biology in general, epigenetic studies in production animals are scarce. Among production animals, chickens have been suggested as a promising model for epigenetic studies (Fresard et al., 2013). Two important reasons for this are that the chicken genome has been extensively sequenced (Rubin et al., 2010) and chickens have historically been an important model for translational research with implications for human health and physiology (Kain et al., 2014).
Long-term stress is known to generate life-long changes in stress susceptibility that are correlated with epigenetic changes (Jensen, 2014). Thus, it is expected that if animals are constantly subjected to stress and systemic hormonal changes, this exposure will imprint the epigenome of, for example, blood cells. Epigenetic changes in blood cells will then serve as markers of past exposure to stress. Research in humans (Malan-Muller et al., 2014) and monkeys (Provencal et al., 2012) has shown that stress affects DNA methylation in blood cells. The epigenome of blood cells can provide a meaningful assessment of biological processes involved in stress because disruptions of the HPA axis have systemic consequences (Zannas and West, 2014). Because different practices in the production environment will generate different levels of stress in animals, it is practical (from the perspective of evaluation of long-term stress) to understand how stress correlates with specific epigenetic profiles in production animals. Chickens provide an excellent model to study epigenetic profiles of long-term stress, as they represent the second-most consumed meat source worldwide (behind pig meat) (OECD-FAO, 2015) and because avian species, unlike mammals, contain nucleated red blood cells (RBCs). This allows for accurate epigenetic profiling because RBCs are a cell type that is simple to purify that can be obtained from live animals.
The present study aimed to investigate whether the rearing conditions of chickens can be identified by looking at epigenetic patterns in their RBCs later in life. We compared RBC DNA methylation in a group of birds reared in cages (a common housing system, with low environmental complexity) with that of birds reared in open aviaries (which represents a complex environment). These two rearing conditions are likely to differ with regard to the amount of stress the birds are exposed to. Previous studies have shown that chickens reared in a complex aviary system are less fearful, use elevated areas of the pen more often as adults (Brantsaeter et al., 2016), and have better spatial working memory (Tahamtani et al., 2015) than laying hens reared in a simpler cage environment. The aviary-reared birds also moved more and spent more time close to a human and a novel object compared with the cage-reared birds, when tested at 19 and 23 weeks of age (Brantsaeter et al., 2016). These results are indicative of aviary rearing reducing the birds’ underlying fearfulness. Additionally, cage rearing was found to reduce the birds’ short-term memory 2 months after transfer from the rearing farm compared with aviary-reared birds, when assessed in a hole-board memory test (Tahamtani et al., 2015). The objective of using this model was to generate a proof-of-concept for future detection of long-term stress in production animals using epigenetic measurements in cell types of easy accessibility in live animals. Our comparison revealed 115 genomic windows with significant differences in RBC DNA methylation between experimental groups, which were located within or in the vicinity of 53 genes and within 22 intronic regions. The identification of a correlation between RBC epigenetic profiles and long-term stress will overcome limitations that exist when evaluating stress through hormonal levels or visual health assessments, which do not provide reliable accounts of long-term stress.
MATERIALS AND METHODS
All experimental protocols employed in the present study that relate to animal experimentation were approved by the Institutional Animal Care and Use Committee at the Norwegian University of Life Sciences under resolution ID number 6190, in order to ensure compliance with international guidelines and regulations for the ethical use of animals in scientific research.
Subjects and rearing treatments
The study was conducted using non-beak-trimmed, female Dekalb white chickens, Gallus gallus domesticus (Linnaeus 1758), aged 0–23 weeks with normal health status. Birds were hatched at a commercial hatchery and immediately brought to a rearing farm. All birds were housed within the same room. Initially, all birds were kept confined inside the aviary row, with access to food and water. When the birds were 4 weeks of age, half of them were given access to the aviary corridors, as this is the normal procedure in aviary-rearing systems. These were named aviary-reared birds (AV). The remaining birds were confined to the aviary row for the entire rearing period; they were named cage-reared birds (CA). These two rearing conditions were maintained until the birds were 16 weeks of age. After the rearing period had ended, a random subset of birds from each treatment was moved to the experimental facilities for blood sampling, which occurred at 24 weeks of age. A schematic representation of the experimental design is shown in Fig. 1A.
Rearing system conditions
The housing system in the single room in which all birds were housed was Natura Primus 1600 (Big Dutchman; www.bigdutchmanusa.com) designed for aviary rearing of laying hen pullets. This system consisted of cages stacked in three tiers placed on either side of a corridor to allow inspection by the caretaker. Cage dimensions were 120×80×60 cm (length×width×height). Each aviary cage contained a 120 cm feed trough, one 120 cm perch and five drinking nipples. All the cages could be opened at the front, allowing the birds to move freely between each tier and the floor of the corridor. Ramps ran from the floor to the second tier to increase ease of access for pullets. When cage doors were in the open position, perches extended from the front of the first and second tiers. The density was 25 birds m−2 for both treatments during the first 4 weeks of life. Chick paper covered 30% of the wire mesh floor of the cages in sufficient amounts to last until the birds were released into the corridors.
During rearing, all birds were exposed to the same light intensity, light schedule and temperatures, as recommended by the General Management Guide for Dekalb White Commercial Layer (Hendrix, 2015). They were provided with ad libitum access to feed using a chain dispersal system and ad libitum access to water. The feed type was conventional pullet feed (Felleskjøpet, Lillestrøm, Norway; Kromat oppdrett 1 for 0- to 6-week-old birds, Kromat avl egg 1 for 6- to 8-week-old birds, and Kromat oppdrett 2 for 8- to 15-week-old birds).
Blood collection and DNA extraction
Blood samples were collected from 21 individuals (9 AV and 12 CA) at 24 weeks of age. Before blood sampling, chickens were sedated using 0.5 ml kg−1 Zoletil mix, which contains 10 ml Rompun (xylazine 20 mg ml−1) and 0.75 ml Butomidor (Butorphanol 10 mg ml−1) mixed with Zoletil powder (tiletamine HCl 125 mg and zolazepam HCl 125 mg). Blood samples were collected as soon as the birds were considered unconscious, which occurred within a maximum time frame of 10 min. Birds were then humanely killed by cervical dislocation. Each blood sample was collected using a 1 ml syringe and a BD Microlance cannula (21 gauge×½ in, 0.80×40 mm). A total of 160 µg of blood was transferred from each sample into two heparinized glass capillaries, which were then centrifuged at 1006 g for 5 min. After centrifugation, the tubes were manually broken into two pieces, one of which contained the haematocrit fraction, which was placed inside 1.5 ml micro-centrifuge tubes and stored in a −80°C freezer until further analyses.
DNA extraction was performed through proteinase K digestion. Initially, 10 µl of the haematocrit fraction was incubated with 200 µl of extraction buffer (1 mol l−1 Tris-HCl, 0.5 mol l−1 EDTA, 10% SDS) and 20 µl of 0.1 mol l−1 DTT at 65°C for 15 min. Then, incubation with 20 µl of proteinase K (20 mg ml−1) was performed overnight at 55°C under rotation. After proteinase K digestion, samples were incubated with Protein Precipitation Solution (Promega) for 15 min on ice and centrifuged for 20 min at 13,000 rpm and 4°C in a benchtop microcentrifuge. The supernatant (1 ml) was transferred to a new tube and DNA was precipitated with an equal amount of 100% isopropanol. In addition, 3 µl of glycogen (5 mg ml−1) was added to improve further visualization of DNA pellets. After 30 min of incubation at 4°C, the samples were centrifuged at 13,000 rpm and 4°C for 30 min. The supernatants were discarded and the DNA pellets were washed with ice-cold 70% ethanol, followed by centrifugation at 13,000 rpm and 4°C for 10 min. The supernatants were discarded again and the pellets in the tubes were dried out in a heating block at 55°C for 5 min. DNA pellets were re-suspended in 200 µl of ultrapure water.
DNA methylation analyses
In order to perform DNA methylation analyses in a cost-effective manner, we combined a genotyping by sequencing method (GBS; Pértille et al., 2016) with the technique of methylated DNA immunoprecipitation (MeDIP; Guerrero-Bosagna and Jensen, 2015). We have recently described the optimization of each of these two methodologies separately for use with chicken DNA (Guerrero-Bosagna and Jensen, 2015; Pértille et al., 2016). This combination of methods was needed because current methods that assess DNA methylation in reduced genomes perform such a reduction through enzymatic digestion targeting restriction sites that contain CpG sites (Gu et al., 2011). Moreover, such an approach is highly biased towards CpG islands (Gu et al., 2011). Our approach, instead, reduces the genome by digestion at restriction sites unrelated to CpGs and is unbiased towards CpG islands.
We first digested the genome with PstI (Thermo Scientific) as previously described (Pértille et al., 2016). After this fragmentation had generated a significantly reduced genome (approximately 2% of its original size) and enrichment of small fragments in a suitable range for Illumina sequencing (200–500 bp) (Pértille et al., 2016), the methylated fraction was captured by an anti-methyl-cytosine antibody (2 µg µl−1; catalogue number C15200006, Diagenode, Denville, NJ, USA) as previously described (Guerrero-Bosagna and Jensen, 2015). The output of the MeDIP was used as the input of GBS. The GBS method uses ligation steps in which a barcode adapter (identifying individual samples) and a common adapter for the Illumina sequencing barcoding system are ligated at each end of the digested DNA fragments (Poland and Rife, 2012). Because of the barcoding system, the GBS technique enables the creation of a sequencing library with DNA pooled from several individuals (Elshire et al., 2011; Poland and Rife, 2012). Once the barcodes and adaptors are ligated, PCR is performed followed by the clean-up of primer dimers and unbound adapters (Elshire et al., 2011; Poland and Rife, 2012). A detailed description of the method for use in chickens has previously been reported (Pértille et al., 2016). The use of the present approach, in which these two methodologies are combined, allowed us to scan the RBC methylome of 21 chickens using only half of an Illumina sequencing lane. Sequencing was performed paired-end with a read length of 125 bp on the Illumina HiSeq2500 platform, at the SNP&SEQ facilities of the SciLifeLab (Solna, Sweden). The dataset supporting the conclusions of this article is available at the European Nucleotide Archive (ENA) repository (EMBL-EBI), under accession number PRJEB21356 (www.ebi.ac.uk/ena/data/view/PRJEB21356).
Single-nucleotide polymorphism (SNP) calling was performed using the Tassel v.3.0 program (Glaubitz et al., 2014) following the default TASSEL-GBS Discovery Pipeline. The alignment of quality-trimmed reads was performed using Bowtie2 tool v.2.2.5 (Langmead and Salzberg, 2012) against the chicken reference sequence (Gallus_gallus 4.0, NCBI). The filtering parameters used were 5% for minimum minor allele frequency and 50% for minimum site coverage.
For the methylated DNA sequencing, data quality trimming was performed in paired-end short reads with the SeqyClean tool v.1.9.10 (https://github.com/ibest/seqyclean) using a Phred quality score ≥24 and a fragment size ≥50. The quality of the reads was checked before and after cleaning by FastQC v.0.11.3 (www.bioinformatics.babraham.ac.uk/projects/fastqc). The Stacks v.1.39 program was used for data de-multiplexing (Catchen et al., 2011) of quality-trimmed reads. In this procedure, each read stored in a FASTQ file has an identification map key file, a barcode containing matching information for the respective sample. The expected reads begin with one of the individual barcodes and are followed by the cut site remnant for PstI, which contains the sequence CTGCA. Fragments are then grouped into individual files, which correspond to individuals identified by their respective barcodes. The option ‘very sensitive-local alignment’ was used in the Bowtie2 tool v.2.2.5 (Langmead and Salzberg, 2012) for the alignment of quality-trimmed reads against the chicken reference sequence Gallus_gallus 4.0 (NCBI). Default parameters for paired- and single-end sequences were used. The coverage depth of each sample was checked using Samtools v.0.1.19 (Li et al., 2009) with the ‘depth’ option.
Because low methylated DNA material is obtained after MeDIP of the PstI-reduced genome, some samples will contribute with very low DNA amounts to be sequenced. These individuals will show a low total number of reads distributed in a few genomic regions, generating a skewed distribution of methylated sites along the genome. This will result in an overestimation of the coverage values in those CpG sites that happened to be covered by reads. To prevent this, we defined a minimum cut-off index in order to select high-quality sequenced samples for further testing of differences between experimental groups. This cut-off index was defined by dividing the percentage of the chicken genome covered (%Cov) by the sequencing coverage average for each sample (S.Cov) and multiplying by 100 (i.e. cut-off index=%Cov/S.Cov×100). It was needed to combine the information provided by each of these variables into one single index because average individual coverage per se could not be used as the sole criteria as some samples have high coverage but only in small factions of the genome. The value of 1.0 was defined as a minimum threshold, which signifies that at least 1% of the genome is well covered in that specific sample. Because previous analyses of the chicken genome by GBS alone showed good coverage of ∼2% of the genome (Pértille et al., 2016), it is reasonable to expect that a smaller percentage than this will be well covered when using GBS in combination with MeDIP. By plotting S.Cov versus %Cov, we obtained an exponential relationship (Fig. S1) in which it can be observed that having a cut-off index of 1 effectively eliminates samples with low %Cov and/or low S.Cov. Based on this, individuals showing a cut-off index below 1.0 were discarded from further analyses.
Following read alignment, all analyses were performed using bioinformatics packages from the ‘R’ Bioconductor repository. The BSgenome.Ggallus.UCSC.galGal4 package was uploaded as the reference genome. The MEDIPS R-package was used for basic data processing, quality controls, normalization and identification of differential coverage. In order to avoid possible artefacts caused by PCR amplification, MEDIPS allows a maximum number of stacked reads per genomic position. This is done by using a Poisson distribution of stacked reads genome-wide. The default parameter of P=0.001 was used as the threshold for the detection of stacked reads. The reads that passed this quality control were then standardized to 100 bp by extending smaller reads to this length, which is the paired-end read size generated by the Illumina HiSeq platform. The genome was divided into adjacent windows of 300 bp length, which is the average length of expected contigs generated by our GBS approach, as well as the program default. MeDIP-seq data were transformed into genome-wide relative methylation scores by a CpG-dependent normalization method (Chavez et al., 2010). This normalization is based on the dependency between short-read coverage and CpG density at genome-wide windows (Down et al., 2008) and can be visualized as a calibration plot. A calibration plot was generated using one of the 10 individuals that passed the cut-off index to generate a coupling set (object that groups information about CpG density genome-wide). Based on this, a threshold for a minimum sum of counts across all samples per window was defined (minRowSum=10). Sequencing data for each individual were then assigned to one of the experimental groups (AV or CA) and differential coverage (i.e. differential methylation) was calculated between the two conditions. Adjacent windows showing significant change were then merged to generate the differentially methylated region (DMR) obtained. For this, the default value of 1 was used within the function MEDIPS.mergeFrames, allowing the neighbouring significant windows to be merged with a 1 base pair gap between them. The merged windows were annotated against the chicken reference genome (Gallus_gallus 4.0, NCBI) using the Variant Effect Predictor (VEP) tool (McLaren et al., 2010). The edgeR and heatmap.2 packages (and extensions) were used for the generation of plots.
The internet-based tool Consensus PathDB (Kamburov et al., 2013) (http://cpdb.molgen.mpg.de) was used to perform an analysis of biological pathways enriched by genes with DMRs found in our study, as well as gene ontology analyses of these genes. Consensus PathDB (Kamburov et al., 2013) integrates interaction networks based on published information in humans. These interaction networks include complex protein–protein, genetic, metabolic, signalling, gene regulatory and drug–target interactions, as well as biochemical pathways (Kamburov et al., 2013). Another internet-based tool used in this study to identify over-represented pathways related to our gene list was Reactome (Croft et al., 2011), which is an open-source curated bioinformatics database of human pathways and reactions (www.reactome.org). The advantage of Consensus PathDB over Reactome is that it is capable of accessing a variety of databases that contain previously described biological pathways (e.g. Kegg, Biocarta, Reactome, Wikipathways). However, in order to use Consensus PathDB, the genes in the chicken genome had to be extrapolated to humans, as it does not accept the ENSEMBL chicken genome annotation. Therefore, Reactome, which did accept the input of chicken genes with the ENSEMBL identifier, was also used. These two tools therefore provided complementary information about our gene list.
The present experiment compared the RBC methylome of chickens reared in open aviaries with that of chickens reared in cages, to detect whether epigenetic profiles in RBCs could be identified as correlating to each of these rearing conditions. The experimental procedure is summarized in Fig. 1B. The RBCs of 21 chickens were extracted in total, nine reared in open aviaries and 12 reared in cages. A combination of the GBS and MeDIP methods was used to identify genome-wide changes in DNA methylation.
The reduced-methylated DNA fraction from RBCs of these animals was sequenced, and bioinformatic analyses were then performed and filter parameters were applied. Our quality control procedure selected sequencing data from four AV and six CA animals for further statistical analyses. In order to account for potential genetic effects generated by the treatments, we tested whether fixed allelic differences could be identified between the experimental groups. A total of 248,170 unique sequence tags were obtained from the 21 individuals submitted to the Tassel pipeline, and 83.1% were aligned against the chicken reference genome (Gallus_gallus 4.0, NCBI). A total of 21,093 SNPs were identified across the 10 individuals (four AV and six CG) that passed the cut-off index. Among the SNPs identified, we did not observe any fixed allelic differences between the CA and AV groups, which could have stochastically appeared as a result of the separation of animals into two groups or emerged as a consequence of the treatment itself.
For the DNA methylation analysis, our method interrogated changes in 810,186 CpG sites per individual, which corresponds to ∼7.6% of all CpGs in the chicken genome. An MA plot showing the log fold-change of AV/CA counts per 300 bp genomic window, which represents changes in DNA methylation, against the normalized window counts is shown in Fig. 2A. Genomic windows with significant changes in DNA methylation between groups (P<0.0005) are indicated by red dots. A principal components analysis (PCA) performed using the windows with significant differences in counts (P<0.0005) between the AV and CA groups confirmed that all individuals in the analysis matched the initial experimental group separation (Fig. 2B). Our comparison revealed that 115 windows showed significant differences in DNA methylation between experimental groups (Table S1). A heat map showing the windows with significant differences is presented in Fig. 3. Adjacent windows showing differential coverage were merged into DMRs between the experimental groups (see Materials and methods for detailed information), which were located within or in the vicinity of (5 kb upstream or downstream of the DMR, based on default criteria in the VEP tool) 53 genes and within 22 intronic regions. Table S2 describes the chromosomal location of all DMRs, the number of CpGs within them, their annotation, as well as their location within or in the vicinity of genes. Fig. 4 summarizes the location of these regions relative to genes (Fig. 4A), as well as their chromosomal location (Fig. 4B). The fold-changes in DNA methylation of the DMR and the direction of these, e.g. hyper- or hypo-methylation of CA versus AV animals, are shown in Fig. 5.
A network analysis was performed with the DMR-associated genes in Consensus PathDB, which connects biological pathways and gene ontology information (Table S3). A simplified pathway showing hypothetical effects of the DMR affected by the treatment is shown in Fig. S2A (redundant information was discarded, e.g. same biological processes shown as being affected by different databases). This analysis showed that DMR-associated genes are mainly enriched in biological processes such as G-protein activation (comprising ∼10% of the genes in that pathway), mitogen-activated protein kinase (MAPK) signalling (involving five genes in our list) and purine ribonucleotide binding (involving 14 genes in our list). P- and q-values of all significantly affected pathways are shown in Table S3. In addition to these main affected pathways, less-enriched pathways are shown in Table S3. Of interest is also the appearance in the network of processes such as visual photo-transduction, opioid signalling, mRNA processing and cytoskeleton organization.
The network analysis performed with Reactome, in turn, shows that genes with altered DNA methylation in our list primarily target pathways in the immune system (Fig. S2B), followed by signal transduction pathways involved in opioid signalling, regulation of the photo-transduction cascade and G-protein activation (Fig. S2C). A less-affected pathway was the metabolic pathway, which showed effects in the sub-pathway inhibition of insulin secretion by adrenaline and noradrenaline. A scheme summarizing the main pathways hypothetically affected is shown in Fig. 6.
Stress has been reported to be associated with DNA methylation alterations in the brain. For example, infant rats exposed to parental maltreatment present long-term DNA methylation and gene expression changes in the brain-derived neurotrophic factor (BDNF) gene in the frontal cortex (Roth et al., 2009). However, from the perspective of using epigenetic tools to determine the history of stress in live animals, it is of interest to determine whether epigenetic changes can also be observed in cell types of easy access such as blood cells.
A few studies have reported epigenetic changes in blood related to stress. For example, adult rats previously exposed to traumatic conditions during early life exhibit an altered microRNA profile in the blood, brain and spermatozoids compared with non-traumatized individuals (Gapp et al., 2014). In humans (Malan-Muller et al., 2014) and monkeys (Provencal et al., 2012), alterations in DNA methylation in peripheral blood cells have been shown to correlate with previous stress. As birds have nucleated RBCs, they represent a model organism in which DNA methylation can be measured in live individuals, and in an easily accessible and simple to purify cell type.
The present study evaluated the effects of early-life conditions on adult DNA methylation patterns in a farm animal. This was performed in RBCs of adult hens that had been reared in groups exposed to different levels of environmental complexity. The aim was to identify epigenetic profiles in RBCs of adult hens associated with different rearing conditions. The rearing conditions to which hens were subjected in the current study cause long-term differences in fearfulness as indicated by differences in inhibition of behaviour and avoidance of a human and a novel object in a novel test arena (Brantsaeter et al., 2016). Although we have not documented stress-related physiological differences between the treatment groups during the rearing phase (first 16 weeks of age), the fact that fear responses are by definition associated with physiological stress suggests that the rearing treatments induce distinct long-term alterations in the stress response. On the one hand, birds in the complex aviary environment are likely to be exposed to a higher degree of mild intermittent stress. On the other hand, confinement in the more barren cage environment may generate a sustained and long-term stress as a result of deprivation. Interestingly, evidence indicates that the aviary environment may be harsher and more challenging than the cage environment, as indicated by the fact that mortality of aviary-housed birds is normally twice as high as that of cage-housed birds (Janczak and Riber, 2015). In addition to fear responses, these rearing conditions are also associated with different levels of cognitive capabilities observed later in life in birds from the same groups as in the present experiment (Tahamtani et al., 2015).
A number of genomic regions presented changes in RBC DNA methylation between the different rearing conditions tested. These DMRs were more common in regulatory regions and less so in intergenic regions (Fig. 4). Such regulatory regions relate to promoters, promoter flanking regions, enhancers, CTCF binding sites, transcription factor binding sites or open chromatin regions, based on the categorization within the VEP functional annotation tool (McLaren et al., 2010). It is well documented that epigenetic mechanisms, particularly DNA methylation, within these regions are major players in the regulation of gene expression (Bogdanović et al., 2016; Wan and Bartolomei, 2008; Weber et al., 2007; Weber and Schübeler, 2007). Interestingly, most of our DMRs are within these regions, and thus have the potential to directly affect gene expression. In addition, a great number of these DMRs are present in intronic regions, which suggests they could be involved in intron retention and splicing, as DNA methylation has recently been reported to have an important role in these processes (Wong et al., 2017). DMRs were absent in chromosomes 3, 9, 14, 18, 23, 32 and W. All the other chromosomes presented a fairly even distribution of DMRs, although chromosome 25 contained the highest number (Fig. 4). The genes associated with these DMRs were tested in pathway network analyses to determine whether they would significantly affect biological processes. It should be emphasized that this network analysis is only a proxy to orient research into the functional relationship of genes in connection with DMRs. Also, it is important to mention that in order to use Consensus PathDB, the genes in the chicken genome had to be extrapolated to humans, as this tool does not accept the ENSEMBL chicken genome annotation. Reactome, however, accepts the input of chicken genes with the ENSEMBL identifier, but it is not as well connected to other databases as Consensus PathDB is. Therefore, these two tools were used to provide complementary information about our gene list.
We used Consensus PathDB and Reactome to search for biological pathways enriched with genes found to be associated with the DMRs reported here. For this, we tested whether at least two of the genes in our list would belong to a single biological pathway previously described in the associated databases. Within Consensus PathDB, we also performed a gene ontology analysis to determine possible common functional roles of these genes. Consensus PathDB analyses demonstrated that differentially methylated genes are involved in pathways related to G-protein activation (in particular, those involved in the opioid response and the photo-transduction cascade), MAPK signalling and purine ribonucleoside binding (related to post-transcriptional processes). MAPKs are known to regulate a wide array of cell functions relating to regulation of gene expression in cellular processes such as proliferation, differentiation, mitosis, apoptosis and survival (Pearson et al., 2001). Interestingly, MAPKs such as p38, MK2 and MK3 are known to mediate the stress response, regulating the transcriptional activation of so-called immediate early genes in mammalian cells (Ronkina et al., 2011). The involvement of purine ribonucleoside (i.e. AMP and GMP) binding has been given minor attention in research investigating stress responses. However, of interest are recent data showing the mediation of purine ribonucleoside binding in the antidepressant side-effects of phosphodiesterase inhibitors (i.e. etazolate, an anxiolytic drug; sildenafil, a drug used in the treatment of erectile dysfunction) in mice (Wang et al., 2014).
Pathway analysis with Reactome gave similar results, as equivalent signal transduction pathways were shown to be affected. The main pathway affected according to Reactome analysis was the immune system. One reason for this effect on the immune system could be that animals living in a confined space exhibit higher levels of stress, which are known to be correlated with an altered immune response. In humans, for example, individuals with a history of post-traumatic stress have compromised immune systems, with a reduced number of lymphocytes and T cells, reduced natural killer cell activity, and reduced production of interferon gamma and interleukin-4 (Kawamura et al., 2001). Also, housing conditions have been correlated with decreased immune response in farm animals. For example, dairy calves housed in smaller stalls present reduced lymphocyte proliferation in comparison to calves in larger stalls (Ferrante et al., 1998). In mice, the bedding type has been shown to influence the intestinal immune system (Sanford et al., 2002).
In addition, many sub-pathways were affected within the immune system. Altered signal transduction pathways include opioid signalling, regulation of the photo-transduction cascade and G-protein activation. Interestingly, opioid signalling has for a long time been related to housing conditions in farm animals. For example, in pigs, opioid receptor density is affected by housing conditions and is inversely correlated with stereotypic behaviour duration (Zanella et al., 1996). Also in pigs, the expression of opioid receptors in the amygdala is substantially different between individuals maintained in an enriched versus conventional housing environment (Kalbe and Puppe, 2010). Although not much research has been done on the role of opioids in chickens, it has been reported that opioid systems modulate social attachment and isolation stress (Sufka et al., 1994; Warnick et al., 2005). This is concordant with findings in rats showing that social isolation increases the responsiveness of the kappa opioid receptor (Karkhanis et al., 2016). An interesting finding in the present paper is that the opioid system could be affected not only in the central nervous system but also in peripheral cells. Further research needs to be done to understand the role of peripheral opioid systems in the modulation of the stress response. Although not many studies have focused on the correlation between photo-transduction and stress, research in chickens has shown that the immune response varies with light cycles in a circadian fashion, controlled in part by the pineal gland, which, among other cell types, contains B-lymphocytes (Bailey et al., 2003). Vasotocin receptors, which belong to the G-protein receptor family, have been reported to mediate the stress response in chickens. The recently characterized neuropeptides in this family (VT2R and VT4R) are known to be involved in the stress response, particularly within the cephalic lobe of the anterior pituitary (Kuenzel et al., 2013). Again, how these neuronal effects translate to peripheral signalling is an interesting matter of future investigation.
In addition, the Reactome metabolic pathway showed some effects on the inhibition of insulin secretion by adrenaline and noradrenaline. Experiments with perfused (canine) pancreas show that insulin secretion is strongly inhibited by adrenaline or noradrenaline (Iversen, 1973). In turn, adrenaline and noradrenaline levels are known to vary not only as a result of stress (Henry, 1992; Ishibashi et al., 2013; Muller et al., 2013) but also in connection with the conditions under which animals are kept in captivity (Muller et al., 2013). For example, in porpoises, free-ranging animals present higher blood levels of both adrenaline and noradrenaline than animals in rehabilitation or under human care (Muller et al., 2013). It is not surprising that high adrenaline or noradrenaline levels in free-ranging animals will lead to the inhibition of insulin and a concomitant rapid increase in circulatory glucose levels, concordant with the high energy demands of animals living in free-ranging conditions. However, it is intriguing that such a mechanisms could be epigenetically regulated.
An interesting suggestion from our data is that a compromised immune system response could be imprinted in the epigenome of RBCs after animals are reared under specific conditions of stress. As DNA methylation patterns are altered, it is suggested that the different rearing conditions leave an epigenetic mark in the RBCs that will in turn affect the functioning of biological processes such as the immune response, possibly in a permanent manner. Further experiments are needed to elucidate whether altered physiological measures of immune responses correlate with developmentally altered epigenetic patterns in farm animals.
The aim of the current study was to identify epigenetic profiles of early developmental stress-related environmental effects in RBCs. We identified distinguishable DNA methylation profiles relating to each treatment. The future goal is that the present results can be used as a proof-of-concept for the identification of epigenetic marks related to past stress conditions that occur in the production environment. Future experiments should evaluate whether sets of DMRs could constitute reliable ‘epigenetic signatures’ of specific and controlled stress conditions in extended populations of animals. The present study reports for the first time DNA methylation changes in RBCs of adult hens when reared in conditions of differing environmental complexity. We describe that these changes in DNA methylation are associated with genes involved in biological functions such as the immune response, and cell signalling related to MAPK, G-protein and opioid pathways. These results prompt interesting questions regarding the role of early-life stimuli in altering epigenetic patterns that could be involved in these mechanisms. Moreover, questions also arise regarding the role RBCs play in G-protein and opioid pathways in the stress response.
Conceptualization: M.B., A.M.J., P.J., C.G.-B.; Methodology: F.P., M.B., J.N., C.G.-B.; Software: F.P., C.G.-B.; Formal analysis: F.P., P.J., C.G.-B.; Resources: M.B., J.N.; Writing - original draft: F.P., P.J., C.G.-B.; Writing - review & editing: F.P., M.B., J.N., L.C., A.M.J., P.J., C.G.-B.; Supervision: L.L.C., A.M.J.; Project administration: A.M.J., P.J.; Funding acquisition: L.L.C., A.M.J., P.J.
F.P. is grateful for funding from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES; Brazilian government) for PhD research internship (PDSE, Programa de Doutourado-sanduíche no exterior) at Linköping University (Sweden), and for post-doc funding from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP, project no. 2016/20440-3). A.M.J., J.N. and M.B. acknowledge that parts of this study involving the use of live animals were funded by the Foundation for Research Levy on Agricultural Products (FFL), the Agricultural Agreement Research Fund (JA), and Animalia (Norwegian Meat and Poultry Research Centre) through the Norges Forskningsråd (grant number 207739). L.L.C. is a recipient of a research productivity scholarship from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq; Brazilian government) and receives funding from Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP; São Paulo State, Brazil; grant 2014/08704-0). C.G.-B. and P.J. appreciates funding from the European Research Council advanced grant 322206 ‘Genewell’ to P.J.
The dataset supporting the conclusions of this article is available from the European Nucleotide Archive (ENA) repository (EMBL-EBI), under accession number PRJEB21356 (www.ebi.ac.uk/ena/data/view/PRJEB21356).
The authors declare no competing or financial interests.