Thermal Manipulation During Embryogenesis Impacts H3K4me3 and H3K27me3 Histone Marks in Chicken Hypothalamus.

Changes in gene activity through epigenetic alterations induced by early environmental challenges during embryogenesis are known to impact the phenotype, health, and disease risk of animals. Learning how environmental cues translate into persisting epigenetic memory may open new doors to improve robustness and resilience of developing animals. It has previously been shown that the heat tolerance of male broiler chickens was improved by cyclically elevating egg incubation temperature. The embryonic thermal manipulation enhanced gene expression response in muscle (P. major) when animals were heat challenged at slaughter age, 35 days post-hatch. However, the molecular mechanisms underlying this phenomenon remain unknown. Here, we investigated the genome-wide distribution, in hypothalamus and muscle tissues, of two histone post-translational modifications, H3K4me3 and H3K27me3, known to contribute to environmental memory in eukaryotes. We found 785 H3K4me3 and 148 H3K27me3 differential peaks in the hypothalamus, encompassing genes involved in neurodevelopmental, metabolic, and gene regulation functions. Interestingly, few differences were identified in the muscle tissue for which differential gene expression was previously described. These results demonstrate that the response to embryonic thermal manipulation (TM) in chicken is mediated, at least in part, by epigenetic changes in the hypothalamus that may contribute to the later-life thermal acclimation.

Introduction

During the last several decades, it became clear that the epigenome of eukaryotes dynamically responds to the environment. Stress, diet, behavior, toxins, and other factors modulate gene expression through epigenetic modifications such as chemical additions to DNA and histones that can be stably maintained during mitotic divisions and do not change the primary DNA sequence (Fresard et al., 2013; David et al., 2017a). Early embryogenesis is a critical window of sensitivity to the environment, and environmentally induced changes can be transmitted during development by subsequent cell divisions. Among the best-studied examples are the persistent epigenetic differences in humans associated with prenatal exposure to famine during the Dutch hunger winter (Heijmans et al., 2008) and the maternal diet supplemented with folate that induces pseudo-agouti phenotypes in mice (Feil and Fraga, 2012). Early temperature exposure is also known to affect the phenotypes of plants and animals through epigenetic reprogramming. In yeast, histone post-translational modifications, and in particular H3K4 methylation, have emerged as primary actors in the mitotic memory of the environmental stress response (Fabrizio et al., 2019). In plants, the molecular basis of vernalization (the acceleration of flowering by prolonged cold) has been extensively studied and notably involves the quantitative accumulation of trimethylation of histone H3 at lysine 27 (H3K27me3) during cold periods (Coustham et al., 2012). Heat acclimation-mediated cross-tolerance resulting from the enhancement of innate cytoprotective pathways in rats was shown to involve epigenetic mechanisms such as post-translational histone modification and altered levels of chromatin modifiers during the acclimation phase (Horowitz, 2017). In chickens, post-natal heat acclimation during the 3rdor 5thday of life was shown to improve temperature tolerance at 10 days of age through epigenetic changes such as DNA methylation at BDNF and histone post-translational modifications at BDNF, Eif2b5 and HSP70 in the hypothalamus (Yossifoff et al., 2008; Kisliouk et al., 2010; Kisliouk et al., 2011; Kisliouk et al., 2017).

Broiler chickens have long been selected for growth performance traits, with significant progress in genetic selection during the last several decades. However, the huge improvement in chicken body weight and muscle growth has not been associated with a similar increase in visceral organs such as heart and lungs necessary for hyperventilation and heat loss (Havenstein et al., 2003). As a consequence, chickens have became less able to cope with extreme environmental temperatures, as illustrated by the 2003 hot spells in France that resulted in the death of several millions of meat-type chickens (Yahav and McMurtry, 2001; St-Pierre et al., 2003). In order to improve bird thermal tolerance and welfare without impairing growth in later life, a protocol of early exposure to heat during egg embryonic development has been proposed (Piestun et al., 2008). The TM treatment during embryogenesis consists in increasing egg incubation temperature to 39.5°C for 12 h/d, in contrast to a continuous incubation at 37.8°C under standard condition, during a time-window that coincides with the development of the hypothalamic-pituitary-thyroid axis (between E7 and E16; chicken eggs hatch around E21) (Piestun et al., 2008; Loyau et al., 2013). The TM treatment, that can easily be applied by hatcheries, was shown to modify several physiological parameters such as lowering the internal temperature and to improve heat tolerance of broilers during acute heat stress at slaughter age (5 weeks old) with little effect on animal performance (Piestun et al., 2008; Loyau et al., 2013).

To investigate the impact of TM at the molecular level, a muscle transcriptome analysis was performed on TM and control (C) chickens that were either heat challenged or not challenged before sampling at 5 weeks of age (Loyau et al., 2016). Only a few genes (28) were found to be differentially expressed (DE) between C and TM. Interestingly, 759 genes were found to be DE between the heat-challenged TM (32°C for 5h at 35 days of age) and TM individuals that were left at thermally neutral condition (21°C) while only 128 were found to be DE in the comparison between heat-challenged controls and controls left in thermally neutral conditions (Loyau et al., 2016). Hence, a number of genes may have been conditioned to respond to further heat exposure by a mechanism that may involve epigenetic changes. In addition, the analysis showed that 26 genes involved in chromatin organization, remodeling, and gene silencing were DE when TM chickens were exposed to a heat challenge (vs. TM in thermal neutral conditions) (Loyau et al., 2016). We therefore hypothesized that an epigenetic reprogramming induced by the embryonic heat exposure may have modified part of the chromatin landscape thus affecting gene expression in favor of heat acclimation. To test this hypothesis, we performed a whole-genome analysis of H3K4me3 and H3K27me3 histone marks. Both P. major muscle and hypothalamus tissues of TM and C chickens were investigated because of the transcriptome changes previously shown and the involvement of brain in postnatal acclimation respectively (Yossifoff et al., 2008). Here, we show that embryonic TM affected H3K4me3 and, to a lesser extent, H3K27me3 in the hypothalamus of 35-day-old chickens whereas the effects were very limited in the muscle. Most differential peaks occurred within the 5’ part of genes that were involved in biological functions relevant to neurodevelopmental, metabolic, and gene regulation functions. This study showed that several genes that code for components of the corticotropin-releasing hormone (CRH) signaling pathway were affected at the histone level, suggesting that TM may affect biological processes similar to those affected by the postnatal acclimation model.

Results

TM Has a Long-Lasting Impact on the Epigenome of Developing Chickens

TM was applied (or not, for controls) on Cob500 eggs and the reproducibility of the embryonic treatment was compared to previous experiments by measuring cloacal temperature at hatching, as this has been shown to be decreased in TM chickens (Loyau et al., 2015). As expected, we found that cloacal temperature of male TM animals was significantly reduced at hatching (38.7°C in average, n = 54) compared to male controls (38.9°C, n = 68; p-value < 0.05). We then sampled hypothalamus and muscle samples of 35-day-old male TM and C chickens. ChIP followed by Illumina sequencing was performed using antibodies specific to H3K4me3 and H3K27me3 modifications. For hypothalamus samples, ∼23 to 75 thousands of peaks were obtained ( and ). As expected, H3K27me3 peaks were broader than H3K4me3 peaks ( and ). 785 and 148 differential peaks (DP) were identified for H3K4me3 and H3K27me3 respectively in this tissue ( ). By contrast, only 45 and 20 peaks were found to be DE for H3K4me3 and H3K27me3 respectively in the muscle, representing less than 0.1% of detected peaks. A majority of DP showed reduced signal in TM animals for both marks and tissues ( ).

Table 1

Summary of ChIP-seq results.Three biological replicates were sequenced per tissue and mark.

	Hypothalamus		Muscle P. major
	H3K4me3	H3K27me3	H3K4me3	H3K27me3
Analyzed peaks	30597	23098	75943	26562
Number of DP	785	148	45	20
% DP	2.57	0.64	0.06	0.08
Number of DP up	147	37	7	4
Number of DP down	638	111	38	16
Genes containing DP	975	69	44	17

Reads were mapped against the chicken genome Galgal5, and peaks were detected using PePr for H3K4me3 and epic for H3K27me3. DP, Differential Peaks. Comparisons (DP up, DP down) were calculated as TM vs. control.

We next investigated the distribution of DP. Hypothalamic H3K4me3 DP were found preferentially in genes, near the transcription start sites (TSS), in the first exon and in the first intron ( ). This result is consistent with genome-wide distribution of this mark in all the H3K4me3 peaks detected in the study ( ). Only DP outside of genes appeared underrepresented suggesting that H3K4me3 changes in the hypothalamus of TM animals occurred preferentially within genes. However, H3K27me3 hypothalamus DP were evenly distributed within the body of genes ( ), likely reflecting the broader size of peaks. However, one quarter of hypothalamus H3K27me3 DP were found outside of genes, which is similar to the genome-wide distribution of this mark in our experiment ( ). Muscle distribution of DP was roughly similar to the hypothalamus, with the exception of the H3K4me3 distribution that appeared less skewed toward the 5’ end of genes ( ).

Figure 1

Characteristics of differential peaks. (A) Percentage of ChIP-seq H3K4me3 differential peaks across genomic features in hypothalamus relative to all differential peaks, obtained using GenomeFeatures. Promoter and downstream regions were defined upstream of Transcription Start and End Sites (TSS and TES, respectively), respectively. (B) Barplot showing the ratio for each feature between the percentage of differential peaks (DP) per feature relative to all DP and the percentage of all peaks per feature vs. all peaks for H3K4me3 in hypothalamus. The color legend is the same as in . (C) Venn diagram showing the number of genes containing a DP that are specific or common for H3K4me3 (K4) and H3K27me3 (K27) marks in hypothalamus (HT) and muscle (M) tissues. (D) Box plot of H3K4me3 relative enrichment to input at CRY2, CDKL5, OCLN, and RGS4 loci in Control (blue) and TM (red) replicates. Box boundaries represent the first and third quartiles. The median is indicated by the bold horizontal line dividing the interquartile range. Upper and lower ticks indicate the 10th and 90th percentiles. Corresponding p-values are indicated at the top of each graph

We next sought to identify DP that were localized within genes (DP genes). For the hypothalamus, we found 975 genes containing at least one H3K4me3 DP (95% unique) and 69 genes containing at least one H3K27me3 DP. Muscle DP were found only in 44 and 17 genes for H3K4me3 and H3K27me3, respectively ( ). Overall, 25 DP genes (2.5%) were found in at least two conditions (shared between marks and/or tissues; and ). Notably, eight genes including the Protocadherin alpha-2(PCDHA2) subunit involved in the establishment and function of specific cell to cell connections in the brain and seven nonannotated genes contained DP in three conditions ( ). Moreover, 17 genes contained DP in two conditions, either in a tissue specific or a mark specific fashion with the exception of another protocadherin (PCDHA11) that contained a H3K4me3 DP in the hypothalamus and a H3K27me3 DP in muscle ( ). The complete list of DP genes is provided in .

To validate the ChIP-seq data, we performed ChIP-qPCR on biological replicates from the same experiment on a set of H3K4me3 DP in the hypothalamus. Primers were designed within or encompassing DP regions at CCNF (cyclin F), CRY2 (Cryptochrome 2), CDKL5 (Cyclin Dependent Kinase Like 5), OCLN (Occludin), and RGS4 (Regulator Of G Protein Signaling 4) with an increase in TM condition. A nondifferential control at the GAPDH locus was included in the analysis. We found that H3K4me3 levels were significantly increased in all regions in TM hypothalamus samples from biological replicates thus confirming ChIP-seq data except for CCNF that only showed a nonsignificant trend toward an increase ( and data not shown). To investigate the relationship between DP and gene expression, we also checked the expression of CCNF, CDKL5, CRY2, and OCLN by RT-qPCR. However, none of these appeared to be altered by the presence of DP as no transcript level was found significantly altered ( ).

Differential Peaks Were Found Preferentially in Neurodevelopmental, Metabolic and Transcriptional Regulator Genes

To gain insights into the biological functions that may be affected by the presence of DP, we performed a functional analysis based on gene ontology (GO) analysis using the ViSEAGO tool (Brionne et al., 2019). Enrichment tests of GO terms associated with biological processes led to the identification of 154 enriched terms among the four list of genes ( , , and and ). These terms were grouped into 26 clusters that could be further categorized into four major processes: neurodevelopment (36 GO terms), metabolism (39 GO terms), gene regulation (30 GO terms) and immunity, cellular and other processes (49 GO terms; ).

Figure 2

Functional analysis of hypothalamus and muscle genes bearing H3K4me3 and H3K27me3 differential peaks.The clustering heatmap plot of functional sets of gene ontology (GO) terms was obtained using ViSEAGO. From left to right are shown the major processes, the conditions defining the clusters (histone marks, tissues) defining the clusters, the cluster number, a heat map with the number of GO terms in each set and a dendrogram based on BMA semantic similarity distance and Ward’s clustering criterion.

Table 2

Description of gene ontology (GO) clusters. For each cluster, the number of GO terms and the common ancestor term are shown

Cluster	GO Terms (total: 154)	Common ancestor term
1	10	Cellular metabolic process
2	9	Regulation of metabolic process
3	5	Regulation of macromolecule metabolic process
4	4	Regulation of macromolecule biosynthetic process
5	4	Regulation of nucleobase-containing compound metabolic process
6	9	Regulation of RNA biosynthetic process
7	6	Organic substance metabolic process
8	4	Nucleobase-containing compound metabolic process
9	7	Cellular biosynthetic process
10	4	Neuron projection morphogenesis
11	6	Anatomical structure morphogenesis
12	7	Gliogenesis
13	7	Cellular developmental process
14	6	System development
15	2	Cerebral cortex cell migration
16	4	Ensheathment of neurons
17	9	Multicellular organismal process
18	2	Positive regulation of muscle hypertrophy
19	5	Negative regulation of biological process
20	6	Microtubule cytosekeleton organization
21	10	Cellular process
22	6	Process
23	11	Biological regulation
24	5	Metal ion transport
25	3	Cytokine secretion
26	3	Protein localization to cell cortex

Six functional clusters comprised GO terms from different conditions. The 7thcluster, associated with gene expression regulation at the protein modification level, was defined by five GO terms corresponding to H3K4me3 DP (four for hypothalamus, one for muscle) and one GO term corresponding to hypothalamic H3K27me3 DP. However, none of the 144 genes defining this cluster were in common between conditions. Moreover, two clusters related to metabolic processes (2 and 23) and one related to cellular processes (21) resulted from GO terms enriched for H3K4me3 DP from both tissues. Finally, two clusters related to cellular processes (17 and 20) resulted from GO terms enriched for both H3K4me3 and H3K27me3 DP in the hypothalamus.

Interestingly, all clusters defining the neurodevelopmental function (10–16) were specific to the H3K4me3 mark and the hypothalamus, suggesting that this mark may contribute to the neurogenic processes of acclimation. This finding was confirmed by an analysis using the Ingenuity Pathway Analysis (IPA) software demonstrating that biological functions associated with hypothalamic H3K4m3 DP were mostly related to the nervous system development, including cellular development, growth and proliferation, as well as the cell-to-cell signaling ( ).

IPA Canonical Pathway Analysis Revealed an Impact of TM on H3K4me3 in Several Genes Related to the CRH Pathway and Previously Identified in the Postnatal Model

To gain further insights into the molecular mechanisms affected by TM, we investigated the canonical pathways using the IPA software and Fisher’s exact test. We focused on the hypothalamic H3K4me3 DP as it was the condition for which most DP were found. The most significantly enriched canonical pathways were mainly related to cellular growth and development such as AMPK signaling and Oct4 (Octamer-binding transcription factor 4) function in embryonic cell, but were also related to the nervous system and behavior with calcium signaling, opioid signaling pathway, and synaptic long-term potential ( ). Interestingly, the CRH signaling pathway was also significantly enriched (p-value < 0.05; and ). Molecules identified in this pathway included ATF2 (Activating transcription factor 2), a transcriptional factor shown to be involved in the maintenance of the chromatin structure following embryonic stress in Drosophila (Seong et al., 2011), and CRHR2 (Corticotropin Releasing Hormone Receptor 2), a receptor of the CRH hormone shown to be involved in determining later-life stress response of thermally-manipulated chicks (Cramer et al., 2015). In addition, the downstream target BDNF (Brain Derived Neurotrophic Factor; ) was previously reported to be a major epigenetic target of postnatal thermal treatment in chicks (Yossifoff et al., 2008). To explore the impact of H3K4me3 changes of several members of the CRH pathway in TM animals, we measured the gene expression by RT-qPCR of several members of the pathway that contained a DP (CRHR2, ATF2 and MAPK11 – Mitogen-Activated Protein Kinase 11; ). We also investigated the expression of other key members of this pathway or downstream targets that did not contain a DP (CRHR1, BDNF, and CREB1–CAMP Responsive Element Binding Protein 1). Our results showed that none of these genes showed expression changes at a significance level of 0.05 even though BDNF and CRHR1 showed a trend toward a decrease (both at P = 0.066; ). Therefore, while TM appeared to significantly affect H3K4me3 marks for a subset of CRH pathway members, we could not demonstrate a clear link with gene expression in the biological replicates assayed.

Table 3

List of the canonical pathways identified by Ingenuity Pathway Analysis. The 25 significant pathways, with p-value < 0.05, are listed. The number and name of molecules participating in each pathway are shown in the “n =“ column

Ingenuity Canonical Pathways	-log(p-value)	n=	Molecules
AMPK Signaling	3.81	21	CHRNA5, ATF2, PPM1B, CAB39, PRKAR2B, PIK3CB, SMARCC1, FOXO6, RAB9B, RAB3A, PPARGC1A, CCND1, SMARCD3, PPP2R5D, CHRNB2, PRKAR1B, PPP2R2B, ELAVL1, PHF10, MAPK11, CHRM4
RAN Signaling	3.57	5	CSE1L, KPNA4, KPNA6, RAN, RANBP1
Role of Oct4 in Mammalian Embryonic Stem Cell Pluripotency	2.79	7	RARA, JARID2, BMI1, NR5A2, WWP2, CCNF, KDM5B
Geranylgeranyldiphosphate Biosynthesis	2.69	2	FDPS, GGPS1
Calcium Signaling	2.68	17	CHRNA5, ATF2, NFATC2, PRKAR2B, CAMK2D, MICU1, TRPC3, CACNA1B, ATP2B1, PPP3CA, ATP2B2, MCU, CHRNB2, PRKAR1B, CACNA2D2, TPM1, RYR2
CDK5 Signaling	2.47	11	PPP1R14C, KRAS, PRKAR2B, ADCY7, ITGB1, PPP2R5D, PRKAR1B, PPP2R2B, MAPK11, GNAL, MAPK6
Trans, Trans-farnesyl Diphosphate Biosynthesis	2.22	2	FDPS, GGPS1
BMP Signaling Pathway	1.83	8	KRAS, ATF2, PRKAR2B, PITX2, TLX2, PRKAR1B, MAPK11, SMURF1
Netrin Signaling	1.81	7	NFATC2, PRKAR2B, PRKAR1B, CACNA1B, PPP3CA, CACNA2D2, RYR2
GP6 Signaling Pathway	1.71	10	SCHIP1, COL4A6, LAMC3, PRKD3, KLF12, COL4A5, COL24A1, COL25A1, PIK3CB, PTK2
Cardiac Beta-Adrenergic Signaling	1.65	11	PPP1R14C, PDE4A, PRKAR2B, ADCY7, PDE4D, PPP2R5D, PRKAR1B, PPP2R2B, AKAP1, PDE10A, RYR2
Opioid Signaling Pathway	1.62	17	KRAS, ATF2, PRKAR2B, PRKD3, ADCY7, CLTC, CAMK2D, CACNA1B, PPP3CA, GNAL, RGS3, MAPK6, PRKAR1B, CACNA2D2, RGS4, GRK6, RYR2
Sonic Hedgehog Signaling	1.61	4	GLIS2, PRKAR2B, PRKAR1B, HHIP
Corticotropin Releasing Hormone Signaling	1.6	11	ATF2, GUCY1A2, PRKAR2B, PRKD3, UCN, ADCY7, PRKAR1B, CACNA1B, MAPK11, CACNA2D2, CRHR2
Synaptic Long Term Potentiation	1.59	9	PPP1R14C, KRAS, GRM2, ATF2, PRKAR2B, PRKD3, CAMK2D, PRKAR1B, PLCL2, PPP3CA
Dopamine-DARPP32 Feedback in cAMP Signaling	1.58	12	PPP1R14C, ATF2, GUCY1A2, PRKAR2B, PRKD3, ADCY7, PPP2R5D, PRKAR1B, PPP2R2B, PLCL2, PPP3CA, KCNJ10
cAMP-Mediated Signaling	1.55	16	GRM2, ATF2, PDE4A, CNGA4, PRKAR2B, ADCY7, CAMK2D, PPP3CA, GNAL, PDE10A, PDE4D, PRKAR1B, CHRM4, RGS4, AKAP1, DUSP4
RAR Activation	1.52	13	RARA, PRKAR2B, PRKD3, ADCY7, PIK3CB, SMARCC1, ALDH1A2, PPARGC1A, SMARCD3, PRKAR1B, PHF10, MAPK11, SRA1
Protein Kinase A Signaling	1.49	21	ATF2, TCF3, PDE4A, NFATC2, CNGA4, PRKAR2B, PRKD3, ADCY7, CAMK2D, PPP3CA, PTK2, PDE10A, ANAPC10, PPP1R14C, PTPRG, PDE4D, PRKAR1B, PLCL2, AKAP1, DUSP4, RYR2
Sperm Motility	1.48	9	PDE4A, GUCY1A2, CNGA4, PRKAR2B, PRKD3, PDE4D, PRKAR1B, PLCL2, PTK2
Regulation of Cellular Mechanics by Calpain Protease	1.46	6	KRAS, CNGA4, CCND1, ITGB1, TLN2, PTK2
Synaptic Long Term Depression	1.44	12	IGF1R, KRAS, GRM2, GUCY1A2, PRKD3, PPP2R5D, PPP2R2B, PLCL2, CACNA1B, CACNA2D2, GNAL, RYR2
ERK/MAPK Signaling	1.41	13	KRAS, ATF2, PRKAR2B, TLN2, PIK3CB, PTK2, HSPB2, PPP1R14C, ITGB1, PPP2R5D, PRKAR1B, PPP2R2B, DUSP4
Renin-Angiotensin Signaling	1.4	10	KRAS, ATF2, PRKAR2B, PRKD3, ADCY7, SHC2, PRKAR1B, MAPK11, PIK3CB, PTK2
Mitotic Roles of Polo-Like Kinase	1.32	5	STAG2, PPP2R5D, PPP2R2B, SMC3, ANAPC10

Muscle Differential Peaks Did Not Localize Within De Genes Previously Identified

Prior to this study, a chicken transcriptome microarray (accession number GSE70756) was performed to compare P. major gene expression profiles of TM and C chickens reared at 21°C and under heat challenge condition at 32°C for 5 h (TM-HC and C-HC, respectively) (Loyau et al., 2016). Among the DE genes in the muscle, 28 were found to be DE between C and TM, 128 between C-HC and C, and 759 between TM-HC and TM. We investigated whether the DP gene list overlapped to these DE gene lists. To that end, the microarray data was reannotated with the Galgal5 annotation used in this study, leading to a similar list of 24, 127, and 724 DE genes included in the analysis for C vs. TM, C-HC vs. C, and TM-HC vs. TM comparisons, respectively ( ). None of the 44 H3K4me3 and 17 H3K27me3 DP genes corresponded to DE genes from all comparisons ( ). We also investigated the presence of DP in candidate genes that were shown to be DE by targeted RT-qPCR approach in the pectoral muscle of TM animals compared to C, including the DIO3 (Iodothyronine Deiodinase 3) and PGC1-α (PPARG Coactivator 1 alpha) metabolic genes (Loyau et al., 2014). However, none of these DE genes contained a DP in the present analysis (data not shown).

Discussion

Many studies investigated the impact of embryonic heat treatment in birds (for a review see Loyau et al., 2015) at the phenotypic, metabolic, and gene expression levels, but the molecular mechanisms underlying these changes were not fully understood. Recently, the hypothesis of an epigenetic contribution was raised by the identification of long-term gene expression effects including the alteration of transcript abundance of key epigenetic players (Loyau et al., 2016). Here, we provided the first demonstration of an epigenetic effect induced during embryogenesis in broiler chickens by the genome-wide mapping of two histone marks, H3K4me3 and H3K27me3, selected for their known function in epigenetic memory.

With only 61 DP for both marks, muscle chromatin of TM chickens displayed little changes compared to control animals at 35 days of age. This surprising result, given the large number of DE genes resulting from TM in muscle of heat-challenged broilers identified in a previous study (Loyau et al., 2016), is unlikely to be due to technical limitations given the very high sequencing depth employed in this study, above 74 M uniquely mapped reads for H3K27me3 instead of 40 M as recommended by ENCODE (Landt and Marinov, 2012) and 41 to 61 M uniquely mapped reads for H3K4me3 whereas ENCODE recommends 20 M for this mark in humans ( ). The number of biological replicates (3) was also in agreement with ENCODE recommendations (> 2). Finally, the ChIP muscle preparation was optimized as shown in a previous study (David et al., 2017b) and led to the identification of a large number of peaks ( ) that were of expected shape and distribution ( and , , and ). Therefore, we are confident in concluding that TM had a limited impact on muscle H3K4me3 and H3K27me3 marks in our study. We previously showed in a transcriptome study based on the same P. major muscle tissue from the same Cobb 500 line of the same age and raised in the same experimental unit, that a large number of genes were DE in heat-challenged animals when they were TM. We hypothesized that some of these DE genes may bear epigenetic marks that may condition their response when animals were further exposed to a heat challenge in later life. Among the 61 genes that contained a DP identified in our study, none were found in the transcriptome data reannotated to match the current annotation. This suggests that both H3K4me3 and H3K27me3 marks are unlikely to be directly responsible for the DE genes resulting from TM previously observed in muscle. Whether or not other epigenetic marks or nonepigenetic mechanisms might be involved in long-term muscle TM response remains to be shown.

TM had a much stronger impact on both marks in the hypothalamus of 35-day-old chickens, with about 150 DP for H3K27me3 and 785 for H3K4me3. About 75% of H3K4me3 DP occurred at the 5’ end of genes, which is mainly explained by the sharp distribution of this mark around gene TSS. Nonetheless, H3K4me3 DP were underrepresented in intergenic regions, suggesting that TM is likely to preferentially impact H3K4me3 peaks possibly controlling gene expression. We confirmed the differential enrichment identified at four candidate loci by ChIP-qPCR in four biological replicates per condition, suggesting that the epigenetic impact of TM was reproducible within the flock. However, among the five DP genes tested for differential expression by qRT-PCR, none seemed to be altered in expression ( and ). Even though H3K4me3 is globally associated with transcriptional activity, many reports also showed that epigenetic marks, in particular when analyzed individually, are not always instructive of transcriptional status (Dong and Weng, 2013; Howe et al., 2017). Instead, TM-induced DP may also promote changes in expression in response to a stimulus (such as heat) rather than a constitutive change in gene expression that may be deleterious for development under normal conditions. For instance, a post-natal heat conditioning of 3-day-old chicks affected the expression of key factors such as BDNF an DNMT3A (Yossifoff et al., 2008), EZH2 (Kisliouk et al., 2011), HSP70 (Kisliouk et al., 2017), CRH (Cramer et al., 2015), and Eif2b5 (Kisliouk et al., 2010) only during a 24-hour heat exposure at 10 days of age but not in control rearing conditions. To test this hypothesis, a new experiment would be required to measure the expression of DP genes in the hypothalamus during a heat challenge at 35 days of age as performed previously (Loyau et al., 2016).

The HSP70 family of chaperone proteins was shown to be involved in the cytoprotection of heat acclimated rats (Horowitz, 2017). In addition, the methylation level of a distal part of the HSP70 promoter was suggested to reflect heat-stress-related epigenetic memory in postnatally heat-conditioned chicks (Kisliouk et al., 2017). All HSP and HSF were checked for DP in all conditions. None appeared to display a differential peak with the exception of HSPB2 for which a significantly higher peak of H3K4me3 signal was found in the hypothalamus. Despite not being a member of the HSP70 family, the small HSPB2 chaperone was recently shown to be associated with neuropathies when disregulated (Yu et al., 2018) and thus may play a role in neuron proliferation, a process shown to be important for acclimation (Matsuzaki et al., 2017).

To explore which biological functions may have been impacted by histone mark changes induced by TM, we performed a comprehensive GO analysis including all conditions (mark, tissue). No significantly enriched GO term was found for muscle H3K27me3, likely due to the very low number of DP genes (17) in that condition. Thirty-nine GO terms identified were related to metabolic processes. This is consistent with findings showing that TM was shown to affect the expression of several metabolic regulator genes in muscle (Loyau et al., 2014) and to the fact that heat acclimation is known to impact metabolic rate in several animal models (Collin et al., 2001; Shein et al., 2007; Piestun et al., 2009). Thirty-six GO terms were associated with neuronal development. IPA biological function analysis confirmed an impact of TM on H3K4me3 peaks localized at genes controlling neurodevelopmental functions. This finding is consistent with previous reports in rats showing that neurogenesis plays an important role in the establishment of heat acclimation in this species (Matsuzaki et al., 2009; Shido and Matsuzaki, 2015). Heat exposure of 5-week-old rats was notably shown to induce GABAergic and/or glutamatergic heat-responsive neurons in the preoptic area of the hypothalamus, influencing autonomic thermoregulation in the long term (Matsuzaki et al., 2017). Thirty GO terms were associated with the regulation of gene expression, which is consistent with a previous hypothesis suggesting that heat exposure during embryogenesis may impact the dynamics of chromatin architecture to allow access of the regulatory transcription machinery to favor efficient response of TM birds to heat (Loyau et al., 2016). Interestingly, one functional cluster was shared by both marks and tissues (excluding muscle H3K27me3, for which no GO term was found as stated above). This cluster is associated with gene regulation at the protein modification level (sialylation, peptidyl-methionine modification or ubiquitination) reinforcing the idea that chromatin structure is a key player for long-term acclimation (Loyau et al., 2016). Finally, 49 GO terms were associated with various cellular processes that notably involve immune responsein agreement with previous findings showing that, in the muscle, genes controlling immune or inflammatory responses were found be DE in TM broilers during heat stress in later life (Loyau et al., 2016). Additionally, in rats, heat acclimation was also shown to alter immune functions (Schneider and Zuhl, 2016). Altogether, TM affected the enrichment level of histone marks at a number of genes involved in several functions previously shown to be impacted by the TM or other heat acclimation models.

One interesting finding highlighted by the IPA H3K4me3 DP analysis was the identification, among others, of the CRH signaling pathway as significantly enriched ( ). The CRH pathway is involved in the neuroendocrine, autonomic, and behavioral stress responses (Kovács, 2013). Several key players of the CRH signaling pathway, including the ones cited below, were identified as bearing a H3K4me3 DP. Among them, CRHR2 encodes one of the two downstream G protein-coupled receptors of CRH involved in the increase of intracellular cAMP levels (Inda et al., 2017). A downstream target in the pathways identified in our study is MAPK11, encoding one of the four p38 MAPKs which was shown to participate to heat acclimation in rats (Horowitz, 2014). A substrate for this enzyme is the transcriptional factor ATF2, and the ATF2 gene also contained a DP. Remarkably, ATF2 was shown to be involved in heterochromatin formation and heterochromatin disruption following heat stress-induced activation (Seong et al., 2011). ATF2 belongs to the cAMP-responsive element (CRE) family that comprises the transcription factor CREB which target genes include BDNF. BDNF was reported more than a decade ago to be a key player in the thermal-experience-dependent plasticity of the hypothalamus in the postnatal acclimation model in chicks (Katz and Meiri, 2006). BDNF changes in expression were shown to be under the control of epigenetic alterations in the postnatal model (Yossifoff et al., 2008). BDNF was not part of the list of DP-containing genes but we tested to see if changes in the CRH pathway may influence its expression in the hypothalamus of TM animals. We therefore tested the gene expression of CRHR2, ATF2, MAPK11 (DP genes) alongside the expression of CRHR1, CREB1, and BDNF. None of these genes appeared significantly misregulated, but BDNF showed a decrease in expression at the significance threshold of 0.1 ( ) that was replicated twice (data not shown). CRHR1 also showed a decrease at the same significance threshold. More biological replicates might be necessary to verify this hypothesis. Nonetheless, it is worth mentioning that BDNF expression was shown to be increased by heat during postnatal conditioning (Katz and Meiri, 2006), and during a heat challenge following conditioning at 10 days of age (Yossifoff et al., 2008), but not after conditioning in absence of heat challenge. It would therefore be interesting to analyze BDNF expression, among other components of the CRH pathway, during or right after a heat exposure as was previously performed (Loyau et al., 2016). Finally, it is worth noting that we previously showed that the corticosterone (stress hormone) response seemed to be reduced in heat stressed TM chickens (vs. TM) compared to heat-stressed control chickens (vs. C) (Loyau et al., 2013). Given the interplay between corticosterone levels and CRH expression (Makino et al., 1994), it remains to be shown whether the changes in CRH signaling were induced directly during TM treatment and/or are a consequence of different stress response of TM animals during their lifetime until sampling.

In conclusion, this study provides the first line of evidence that embryonic TM in chicken involves lasting epigenetic changes in the hypothalamus that may contribute to the thermal acclimation of the animals.

Methods

Animals

Experiments were performed in accordance with the legislation governing the ethical treatment of birds and were approved by the French Ministry of Higher Education and the Val De Loire Animal Ethics Committee (Authorization N°APAFIS#4608-201603211212171v2). Cobb 500 eggs were supplied by Hendrix Genetics (Saint-Laurent-de-la-Plaine, France). The 500 eggs were incubated and hatched and animals were raised at the INRA UE1295 PEAT Poultry Experimental Facility (2018, https://doi.org/10.15454/1.5572326250887292E12). Half of the eggs were incubated in TM conditions (39.5°C for 12h/day and 65% RH during E7-16, the rest of the time at 37.8°C and 56% RH) and the other half was incubated in standard condition (37.8°C, 56% RH). The 120 TM chicks and 143 control chicks were hatched, monitored for cloacal temperature using a KIMO KTT-310 probe thermometer (KIMO, Montpon Ménestérol, France) and sexed by venting. Male chickens were raised in four pens (two for TM and two for C) and slaughtered for Pectoralis Major muscle and hypothalamus sampling at 35 days of age. Statistical analysis of temperature at hatching was performed using a Student test. Six animals were used for muscle ChIP-seq experiments (three TM and three C) and six animals were used for hypothalamus ChIP-seq experiments (three TM and three C). For ChIP-qPCR validations, hypothalamus samples from six male animals (three TM and three C; 35 days old) from the same experimental batch were sampled. For RT-qPCR analysis, hypothalamus samples from 10 other male animals (six TM and four C; 35 days old) from the same experimental batch were sampled.

Chromatin Immunoprecipitation, Library Construction, and Sequencing

Muscle native chromatin immunoprecipitation (ChIP) was performed as previously described (David et al., 2017b). For hypothalamus samples, tissues were ground in a mortar using a pestle both cooled using liquid nitrogen. Ground powder was transferred in 1.5 ml TPX tube containing fixation buffer (1% FA in PBS with 1 X Complete™ protease inhibitors, PBS-C, Roche) for 5 min at room temperature under agitation. The reaction was stopped by adding Glycine to 0,125M final. The samples were centrifuged for 5 min at 3,000 g and washed twice using 1-ml ice-cold PBS-C. Cells were then centrifuged 5 min at 5,000 g and cellular pellets were resuspended in 300 µl sodium Dodecyl Sulphate (SDS) Lysis Buffer (50 mM Tris-HCl pH 8, 10 mM EDTA pH 8, 1% SDS, 10% glycerol, 1 X Complete™ protease inhibitors). Tubes were vortexed 20 s, incubated 10 min on ice and vortexed 20 s again. Lysed nuclei were sonicated using a Bioruptor (Diagenode, Denville, USA) set to high setting (30 s ON, 30 s OFF; 4 × 5 min). Chromatin immunoprecipitation was then performed as previously described (David et al., 2017b) using 5 µl of H3K27me3 antibody (07–449, batch #2506493, Merck-Millipore, Billerica, USA) or 5 µl of H3K4me3 antibody (07–473, batch #JBC188194 for hypothalamus and #2717639 for muscle, Merck-Millipore, Billerica, USA). The Sequencing libraries were prepared using the NEB Next Ultra II DNA Library Prep kit for Illumina (New England Biolabs, Ipswich, USA). Ten ng of H3K27me3 IP DNA, 10 ng of H3K4me3 IP DNA and 100 ng of input DNA were used. For hypothalamus, 10, 7, and 4 cycles of PCR amplification were performed, respectively; for muscle, 7, 8, and 4 cycles of PCR amplification were performed, respectively. Agencourt AMPure XP beads (Beckman Coulter, Brea, USA) were used for the 200-bp size selection of DNA fragments. Library concentration was measured by Qubit™ dsDNA HS Assay Kit (Invitrogen, Carlsbad, USA).

All samples were sequenced by the GenomEast genomic platform (IGBMC, Illkirch, France). Sequencing was realized in single-end read with 50-bp fragments by an Illumina Hi-Seq 4,000 sequencer.

Bioinformatic and Biostatistic Analyses of Peaks

Sequencing quality was verified by a FASTQC analysis version 0.11.8 (Babraham bioinformatics). All reads were mapped to the version 5 of the chicken reference genome (Gallus gallus 5.0, GenBank assembly accession GCA_000002315.3) using Bowtie 2 version 2.2.6.2 (Langmead and Salzberg, 2012) with default options. Read duplicates were removed using SAMtools version 1.7 (Li et al., 2009). H3K27me3 broad peak detection for both tissues was realized using epic version 0.1.25 (https://github.com/biocore-ntnu/epic), a reimplementation of SICER (Xu et al., 2014) (options: –fragment-size 50–gaps-allowed 2–false-discovery-rate-cutoff 0.05). For hypothalamus samples, H3K4me3 peak detection was realized using PePr (Zhang et al., 2014) (version 1.1.14) and default settings (option: –peaktype = SHARP). For H3K4me3 peaks in muscle samples, PePr window size of peak detection was fixed to 240 bp (options: –peaktype = SHARP -w 240) as peaks were too broad using default settings in this tissue. A union between controls and TM peaks was performed independently for each tissue to include in the analysis peaks specific from each treatment and peaks common between treatments. Counts at union peak level were computed using featureCounts (Liao et al., 2014) version 1.6.2 (options -F SAF -s 0).

Differential peaks (DP) between control and TM samples were identified using edgeR package version 3.20.8 (McCarthy et al., 2012) for the four data sets (H3K27me3 and H3K4me3 from muscle and hypothalamus samples). Raw counts were normalized for library size and peak composition using trimmed mean of M-values method as recommended by Dillies et al. (2013). Descriptive and diagnostic graphs, like MDS-plot, permitted to detect a technological batch effect corresponding to library preparation and sequencing. A generalized linear model was fitted on normalized counts for each peak including batch and challenging factors (TM vs. C). A likelihood ratio test was performed to test the effect of the challenging factor. Raw p-value were adjusted for multiple testing using the Benjamini-Hochberg procedure to control False Discovery Rate (Benjamini and Hochberg, 1995). Peaks with an adjusted p-value less than 0.05 were considered as differential.

Annotation of DP and analysis of their distribution was made with the GenomeFeatures R Package version 0.99 (https://forgemia.inra.fr/aurelien.brionne/GenomeFeatures) for promoter (3 Kbp upstream and 500 bp downstream of TSS), flanking exons (5’UTR, 3’UTR), exons, introns, and downstream (1 Kbp from TES) features on both strands. A DP was considered within a gene when at least 1 bp of the peak overlapped to the gene, from the TSS until the transcriptional end site (TES). Peaks were visualized with Integrative Genome Viewer (Robinson et al., 2011). BigWig files were obtained from bam files using bamCoverage from deepTools suite (Ramírez et al., 2016) version 3.0.0. Normalization factors were calculated by performing a ratio between the bam file which have the smallest number of unique reads aligned per condition and the bam file of interest of the corresponding condition. Then, mean of BigWig files in each condition and each tissue was realized with bigWigMerge and bedGraphToBigWig (UCSC bioinformatic tools (Kent et al., 2010).

Functional Analysis

Biological interpretations were carried out using the GO public database with the use of Biological Process (BP) category using ViSEAGO R package (Brionne et al., 2019). Associated Gene terms were retrieved from EntrezGene database (Sept. 2018) for Gallus gallus and experimental GO terms from orthologous genes were added. Enrichment tests were performed using exact Fisher’s test for each mark and tissue (four lists of genes). Enriched GO terms (p < 0.01) were grouped into functional clusters using hierarchical clustering based on Wang’s semantic similarity between GO terms respecting GO graph topology and Ward’s criterion. To go further in the interpretation, these functional clusters were grouped using hierarchical clustering based on BMA distance between sets of GO terms and Ward’s criterion.

Ingenuity Pathway Analysis (IPA, Ingenuity Systems, Mountain View, CA; version 377306M, content version 27216297) was used to explore canonical pathways. Core analysis, originally developed to analyze the impacted biological functions from expression data set, was used on identified DP genes with gene symbol and the associated adjusted p-value to discover the impacted molecular and cellular functions.

ChiP-qPCR

Hypothalamus were mechanically homogenized in 1 ml of 0,1X fixation buffer in PBS with 1 X Complete™ protease inhibitors (Roche) using 5-mm stainless steel beads (Retsch) and a Mixer Mill MM 400 (Retsch, 1 min at 10 Hz). The rest of the experiment was performed as the ChIP-seq experiments previously described using H3K4me3 antibody (07–473, batch #2839113, Merck-Millipore, Billerica, USA). Immunoprecipitated DNA was eluated in 60 µl of TE buffer. A selection of putative candidates was obtained from ChIP-seq DP data. DP selection was based on the adjusted p-value (< 0.05) and a fold change >2 or <0.5. qPCR primers were designed in or around the DP coordinates with NCBI Primer blast (https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK_LOC=BlastHome) ( ). The qPCR was performed as previously described (David et al., 2017b) on 2 µl of ChIP DNA or 1 µl of input DNA in triplicates with 5 µM of each primer and Takyon No ROX SYBR 2X MasterMix blue dTTP (Eurogentec, Liege, Belgium) following manufacturer’s instructions. Reactions were performed on a LightCycler®480 Instrument (Roche Diagnostics, Basel, Switzerland) using the following program: denaturation 5 min at 95°C, 50 amplification cycles (10 s at 95°C, 15 s at 60°C, and 15 s at 72°C), melting curve (5 s at 95°C, 1 min at 65°C, continuous at 95°C) and cooling. Enrichments were determined with the percent input method [100*2^ (adjusted input - Ct (IP)]. Statistical analyses of ChIP-qPCR were performed using R (version 3.5.1) by comparing TM (n = 3) and C (n = 3) enrichment using a Wilcoxon signed rank test with one-sided (alternative hypothesis as less) at the 0.05 level of significance.

RNA Isolation, Reverse Transcription and Quantitative PCR for Gene Expression

Gene expression was evaluated by reverse transcription followed by real-time PCR as previously described. Total RNA was extracted from hypothalamus of biological replicates using AllPrep RNA/DNA Mini kit (Qiagen), according to the manufacturer’s instruction. RNA integrity was verified by electrophoresis on agarose gel. RNA was quantified on a Nanodrop ND-1000 UV-Vis Spectrophotometer.

cDNAs were synthesized from 2 µg of total RNA using the Superscript II enzyme (Invitrogen) and hexamer random primers (Promega) in 20 µl final volume, following the manufacturer’s instructions. Quantitative real-time PCR (qPCR) was carried out with Takyon qPCR Kits (Eurogentec) and using a LightCycler®480 Instrument II system (Roche) with 384-well plates (4TI-0382, 4titude).

Three technical replicates were performed for each sample and a standard curve protocol was used to evaluate gene expression. Primer sequences ( ) were designed with NCBI Primer blast or Primer3Plus (http://www.bioinformatics.nl/cgi-bin/primer3plus/primer3plus.cgi) software. To assess the amplification of the correct cDNA fragments, every PCR product size was checked on a 2% agarose gel and sequenced by Sanger sequencing (Genewiz). Relative expression was normalized to the expression of three reference genes (selected among four reference genes tested): GAPDH, POLR2E, and CytB using geNorm software (Vandesompele et al., 2002). Statistical analyses of gene expression were performed using R (version 3.5.1) by comparing TM (n = 6) and C (n = 4) enrichment using a Wilcoxon signed rank test with two-sided at the 0.05 level of significance.

Data Availability Statement

The datasets generated for this study can be found in the ChIP-seq data have been deposited in the ArrayExpress database at EMBL-EBI (www.ebi.ac.uk/arrayexpress) under accession number E-MTAB-7297 (“ChIP-seq of H3K4me3 and H3K27me3 histone marks in Pectoralis major muscle of 35 days-old control and thermally-manipulated chickens”) and E-MTAB-7300 (“ChIP-seq of H3K4me3 and H3K27me3 histone marks in the hypothalamus of 35 days-old control and thermally manipulated chickens”).

Ethics Statement

The animal study was reviewed and approved by the French Ministry of Higher Education and the Val De Loire Animal Ethics Committee (Authorization N°APAFIS#4608-201603211212171v2).

Author Contributions

VC, AC, and S-AD designed the experiments. VC, S-AD, AVC, TB, SC, and NC performed the experiments. VC, CG, BP, CH-A, AB, AVC, YB and S-AD analyzed the results. VC, AVC, CG, AB, and CH-A wrote the manuscript. All authors read and approved the final manuscript.

Funding

Studies were made possible by the financial support of INRA Department Animal Physiology and Livestock Systems (PHASE), project STRESSEPIMARK. S-AD PhD work was cofunded by INRA-PHASE and Region Centre-Val de Loire. AVC and CG were funded by the ANR JCJC “QuailHeatE” Research Program (ANR-15-CE02-0009-01).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.