Sex-Specific Methylomic and Transcriptomic Responses of the Avian Pineal Gland to Unpredictable Illumination Patterns

Abstract

In the production environment of chickens, exposure to unpredictable light patterns is a common painless stressor, widely used to influence growth rate and egg production efficiency. The pineal gland, a key regulator of circadian rhythms through melatonin secretion, responds to environmental light cues, and its function is modulated by epigenetic mechanisms. In this study, we investigated how the pineal gland methylome and transcriptome (including micro-RNAs) interact to respond to a rearing exposure to unpredictable illumination patterns, with a particular focus on sex differences. We conducted an integrative multi-omic analysis-including methylomic (MeDIP-seq), transcriptomic (RNA-seq), and miRNA expression profiling-on the pineal gland of Hy-Line White chickens (n = 34, 18 females, 16 males) exposed to either a standard 12:12 light-dark cycle (control) or a randomized, unpredictable light schedule from Days 3 to 24 post-hatch. Our findings reveal that unpredictable light exposure alters the pineal gland methylome and transcriptome in a sex-specific manner. However, while transcriptomic differences between sexes increased due to the stress, methylomic differences decreased, particularly on the Z chromosome. These changes were driven by females (the heterogametic sex in birds), which became more male-like in their pineal methylome after exposure to the illumination stress, leading to reduced epigenetic sexual dimorphism while maintaining differences at the gene expression level. Further, we implemented a fixed sex effect model as a biological proof of concept, identifying a network of 12 key core genes interacting with 102 other genes, all linked to circadian regulation and stress adaptation. This network of genes comprises a core regulatory framework for circadian response. Additionally, tissue-specific expression analysis and cell-type specific expression analysis revealed enrichment in brain regions critical for circadian function, including neuronal populations involved in circadian regulation and the hypothalamic-pituitary-thyroid axis. Together, these findings provide strong evidence of sex-specific epigenetic transcriptomic responses of the pineal gland upon illumination stress and offer valuable insights into the interplay of different omic levels in relation to circadian response.

Keywords: circadian; epigenetics; gene expression; illumination; light; methylome; miRNA; pineal gland; sex differences; transcriptomic.

Figure 1

Treatment process of randomly divided newly hatched chickens. The control group was kept at a standard 12:12 light–dark cycle for their entire life. The chronic light stress group was exposed to unpredictable illumination patterns on days 7–24 of age.

Figure 2

Representation of multidimensional scaling (MDS) analysis highlighting distinctions based on sex (a; left), treatment (b; center), and their combination (c; right), through the integration of DEG (top), DMR (center), and DMiR (bottom) data. As object in the analysis, we used the genes from DEG (), windows from DMR (), and miRNA from DMiR () analysis.

Figure 3

Number of differentially expressed genes (DEG) (a), differentially methylated regions (DMRs) (b), and differentially expressed microRNAs (DMiR) (c) across the five contrasts analyzed in the study (SvC, MSvFS, MCvFC, MSvMC, FSvFC), with a significance threshold of p ≤ 0.05.

Figure 4

Barplot illustrating the distribution of omic features across the 5 contrasts (MSvFS, MCvFC, MCvMC, FSvFC, and SvC) and their overlaps, featuring (a) differentially expressed genes (DEG), (b) genes related to differentially methylated regions (DMRs), and (c) target genes identified from differentially expressed microRNAs (DMiR) analysis.

Figure 5

Annotation of DMRs showcases their connection to gene regions, such as intronic, exonic, promoter, and distal areas (a), and their distribution across chromosomes, with a focus on the Z chromosome (b), across the various contrasts explored in this research.

Figure 6

Venn diagram comparing DMR counts between sexes exposure to long‐term stress and annotation against the chicken reference genome of lost (left), retained (center), and gained (right) DMRs between sex of the stressed group in relation to the control group (a). Annotation of control and stressed DMRs relative to chicken genes, considering chromosomal location (b).

Figure 7

Heatmaps showing DMRs categorized as loss (a), gain (b), and remnants (c) between control (blue) and stress (yellow) conditions, stratified by sex (black: female; gray: male), with red intensity indicating methylation levels.

Figure 8

Venn diagram illustrating the distribution of omic features, differentiating DMiRs, DEGs, and DMRs, across the sex fixed model (SvC). The accompanying bar plot categorizes genomic features based on the position of DMRs relative to their nearest annotated gene, highlighting their distribution alongside DMiRs and DEGs.

Figure 9

Relative levels and direction of significant omic changes in genes containing modifications in two omic levels: (a) DEG‐DMR, (b) DMR‐DMiR, and (c) DEG‐DMiR. Positive fold changes indicate that the transcriptome (DEG), miRNA (DMiR), or DMR‐related gene is higher expressed or hypermethylated in the stress group compared to the control group.

Figure 10

Venn diagram showing (a) unique and intersected genes among different omic level features (p < 0.05) obtained using SvC after the module approach,(b) gene enrichment using KEGG pathways, and (c) GO terms.

Figure 11

Gene–protein interaction network obtained from the module inference approach from DEG (green), DMR (blue), and DMiR (red) of all the inputted genes from the module approach. The inner circles represent the original genes from the three genomic levels, and the outer circles represent the genes that emerged from the module interaction inference. Red lines represent the interactions between the 12 core genes and other genes.

Figure 12

Tissue specific expression analysis (TSEA) (a) and cell‐type specific expression analysis (CSEA) (b) of the 114 intersected genes among different omic levels (p < 0.05), obtained using SvC after the module approach.