Integrated multi-omics data to identify the key genes affecting egg-laying performance in late-phase hens.

Hongrui Yang: Writing – review & editing, Writing – original draft, Visualization, Validation, Data curation. Lei Wang: Resources, Investigation. Li Rong: Visualization, Formal analysis. Jiuhong Nan: Formal analysis. Jinping Cai: Investigation. Yunxia Zhao: Supervision, Funding acquisition. Guiyu Zhu: Methodology, Conceptualization. Shijun Li: Writing – review & editing, Supervision, Project administration, Methodology, Funding acquisition, Conceptualization.

To characterize egg production patterns during the late laying phase, laying rates from weeks 44 to 49 were recorded for 872 Wenchang chickens. The results revealed a population mean laying rate of 57.0 % (SD = 0.178), with a coefficient of variation of 31.1 %, indicating substantial variation in late-phase productivity. Statistical analysis demonstrated that the egg-laying rates in this population followed an approximately normal distribution ( Fig. 1 A). A total of 151 hens (17.3 % in total) were classified as high egg production (laying rate ≥ 76.0 %), 511 (58.6 % in total) with medium egg production (50.0 % < laying rate < 76.0 %), and 210 (24.1 % in total) with low egg production (laying rate ≤ 50.0 %). Notably, despite the generally low and differential productivity, 12 (1.3 % in total) high-performing individuals maintained laying rates > 85.0 %. To investigate the mechanisms underlying this variation, we established experimental the HEL (n = 3, mean egg-laying rate = 78.4 %) and LEL (n = 3, mean egg-laying rate = 36.3 %) groups. Statistical analysis confirmed that a highly significant difference between the two groups (p = 0.005; Fig. 1 B). Fig. 1 Population egg-laying performance and ovarian histology. (A) Histogram of egg-laying rate for Wenchang chickens during the late laying phase (44 to 49 weeks of age, n = 872). (B) Egg-laying rates in the HEL and LEL groups. Data are presented as mean ± SD. ** p < 0.01. (C) H&E-stained paraffin sections of ovarian tissue from the HEL and LEL groups. H, the HEL group; L, the LEL group; C, cortex; M, medulla; *, white follicles; △, atretic follicles. Some small white follicles and primordial follicles are not labeled. Fig 1 dummy alt text Population egg-laying performance and ovarian histology. (A) Histogram of egg-laying rate for Wenchang chickens during the late laying phase (44 to 49 weeks of age, n = 872). (B) Egg-laying rates in the HEL and LEL groups. Data are presented as mean ± SD. ** p < 0.01. (C) H&E-stained paraffin sections of ovarian tissue from the HEL and LEL groups. H, the HEL group; L, the LEL group; C, cortex; M, medulla; *, white follicles; △, atretic follicles. Some small white follicles and primordial follicles are not labeled. To investigate histological differences in ovaries associated with divergent egg production, we prepared H&E-stained paraffin sections from the ovarian stroma of the HEL and LEL groups. Histological analysis revealed that the HEL ovaries contained numerous follicles of varying sizes with compact cortical and medullary structures. In contrast, the LEL ovaries exhibited depleted follicular reserves and sparse cortical and medullary structures ( Fig. 1 C). To investigate the determinants of egg-laying variation at the transcriptomic level, we performed RNA-seq library construction, sequencing, and analysis on ovarian stroma from the HEL and LEL groups. Data quality control results demonstrated the acquisition of > 39.8 million clean reads per sample and genomic mapping rates exceeding 92 % ( Table S1 ). Visualization of FPKM distributions revealed consistent gene expression patterns across all six libraries ( Fig. S1A and S1B ). PCA further identified significant correlations between egg-laying rates and PC2, PC3, and PC9 through linear regression modeling ( Fig. S1C and 1D ). To further investigate transcriptional differences in the ovarian stroma of the HEL and LEL groups, DEGs were identified with thresholds of p < 0.05 and |log₂(FC)| ≥ 1. Among 23,726 genes detected in transcriptome data, 800 DEGs were identified, comprising 317 upregulated and 483 downregulated genes in the HEL group ( Fig. 2 A, Table S2 ). These DEGs exhibited distinct clustering in hierarchical heatmap analysis ( Fig. 2 B). Nine randomly selected DEGs showed consistent expression trends in both RNA-seq and RT-qPCR assays ( Fig. 2 C). Furthermore, comparative analysis of fold-change measurements also revealed a highly positive correlation (R² = 0.8917, Fig. 2 D). Besides, 440 of the 800 DEGs represent unannotated genes requiring extensive functional characterization. Fig. 2 Identifying of DEGs in the HEL and LEL groups. (A) Volcano plot of DEGs. The screening criteria were p < 0.05 and |log 2 FC| ≥ 1. Red dots: upregulated genes; blue dots: downregulated genes; grey dots: non-significant genes. (B) Clustering heatmap of gene expression in different samples. Color reflects the expression level of a DEG in the sample. (C) RT-qPCR validation. Data are presented as mean ± SD. *p < 0.05, ** p < 0.01. Specificity validation for the primers used is shown in Fig. S2 . (D) Correlation analysis between RNA-seq and RT-qPCR results. Each black dot represents a gene, and the dashed line indicates the linear trend. Fig 2 dummy alt text Identifying of DEGs in the HEL and LEL groups. (A) Volcano plot of DEGs. The screening criteria were p < 0.05 and |log 2 FC| ≥ 1. Red dots: upregulated genes; blue dots: downregulated genes; grey dots: non-significant genes. (B) Clustering heatmap of gene expression in different samples. Color reflects the expression level of a DEG in the sample. (C) RT-qPCR validation. Data are presented as mean ± SD. *p < 0.05, ** p < 0.01. Specificity validation for the primers used is shown in Fig. S2 . (D) Correlation analysis between RNA-seq and RT-qPCR results. Each black dot represents a gene, and the dashed line indicates the linear trend. Significant enrichment (p < 0.05) was detected for 144 and 111 GO terms among upregulated and downregulated DEGs, respectively. After multiple-testing correction, 11 upregulated and 29 downregulated terms remained significant (padj < 0.05; Fig. 3 A and 3B; Table S3 ). Among all DEGs, the term “integral component of membrane” contained the largest number of genes, including 29 upregulated genes, and 40 downregulated genes. This was followed by the term “extracellular space”, which contained 20 upregulated genes, and 23 downregulated genes. Notably, many significantly enriched GO terms were related to immune processes, including “immune response”, “defense response to bacterium”, “killing of cells of other organism”, and “CCR6 chemokine receptor binding”. Fig. 3 Function enrichment analysis of differentially expressed genes. (A to B) GO functional enrichment for upregulated (A) and downregulated (B) DEGs. The top 10 GO terms with an adjusted p-value (padj) < 0.05 in each category are visualized. (C) KEGG pathway enrichment analysis of DEGs. Dot size and color indicating the number of enriched genes and the significance of enrichment, respectively. (D) PPI network of DEGs. Each node represents a DEG. Node size and color represent its degree. The innermost nodes correspond to the top 20 % of DEGs by degree, the middle layer to the 20 % to 53.5 %, and the outermost layer to the 53.5 % to 100 %. Fig 3 dummy alt text Function enrichment analysis of differentially expressed genes. (A to B) GO functional enrichment for upregulated (A) and downregulated (B) DEGs. The top 10 GO terms with an adjusted p-value (padj) < 0.05 in each category are visualized. (C) KEGG pathway enrichment analysis of DEGs. Dot size and color indicating the number of enriched genes and the significance of enrichment, respectively. (D) PPI network of DEGs. Each node represents a DEG. Node size and color represent its degree. The innermost nodes correspond to the top 20 % of DEGs by degree, the middle layer to the 20 % to 53.5 %, and the outermost layer to the 53.5 % to 100 %. A total of 12 KEGG pathways showed significant enrichment in all DEGs (p < 0.05; Fig. 3 C; Table S4 ). Among them, the most notable one was the neuroactive ligand-receptor interaction pathway (p < 0.001). Three enriched pathways are involved in cellular communication with the external environment and the maintenance of tissue homeostasis, including ECM-receptor interaction, focal adhesion and tight junction. It is also enriched in a vascular smooth muscle contraction pathway related to physiological activity. In addition, there were seven related to substance synthesis and metabolism, including tyrosine, linoleic acid, galactose, fructose and mannose, retinol, and arachidonic acid. To identify potential key regulatory genes, we constructed PPI network from the DEGs ( Fig. 3 D). Based on degree ranking, the top 20 % of genes (38 in total) were identified as hub genes ( Table S5 ), suggesting their extensive interactions and potential pivotal contribution in the overall regulatory network. Moreover, preceding functional enrichment analysis revealed that these 38 hub genes were significantly associated with multiple processes involved in follicular development and ovulation regulation, including immune response, negative regulation of endopeptidase activity, response to corticosterone, antimicrobial humoral immune response mediated by antimicrobial peptide, defense response to gram-negative bacterium, extracellular space, integral component of plasma membrane, plasma membrane, motor activity, serine‑type endopeptidase inhibitor activity, tight junction, tyrosine metabolism, neuroactive ligand-receptor interaction, vascular smooth muscle contraction and focal adhesion ( Table 1 ). Table 1 Functional enrichment analysis of hub genes. Table 1 dummy alt text Category Description Upregulated Genes Downregulated Genes GO-BP Immune Response MHCIA6, HLA-F10AL3, MICA, MHCIY - Negative Regulation of Endopeptidase Activity OVAL, OVALY, KNG1, C4A - Response to Corticosterone OVAL, BPIFB2 - Antimicrobial Humoral Immune Response Mediated by Antimicrobial Peptide - CATH1, TF Defense Response to Gram-negative Bacterium - CATH1, TF GO-CC Extracellular Space OVALY, KNG1, MHCIA6, C4A, OVAL, HLA-F10AL3, MICA, MHCIY DBH, GAL, CATH1, TF, CG-16 Integral Component of Plasma Membrane OVALY, MHCIA6, MHCIY, CYP3A5, HLA-F10AL3, MICA DBH, SLC17A8, EGF, CD247, AGTR1 Plasma Membrane - EGF, CD247, AGTR1 GO-MF Motor Activity MYH1B, MYH1D, MYH1F - serine-type endopeptidase inhibitor activity OVAL, OVALY - KEGG Tight Junction MYH1B, MYH1D, MYH1F MYH11 Tyrosine Metabolism - DDC, DBH, TH Neuroactive Ligand-receptor Interaction KNG1 AGTR1, GAL Vascular Smooth Muscle Contraction - ACTG2, AGTR1 Focal Adhesion - VTN, EGF Functional enrichment analysis of hub genes. To investigate the causes of egg-laying variations from an epigenomic perspective, we constructed ATAC-seq libraries using ovarian stroma tissues from the two individuals with the highest and lowest egg-laying rates within each experimental group. Data quality statistics showed that each sample yielded more than 35.8 million clean reads, with mapping rates to the genome exceeding 98 % ( Table S6 ). Spearman correlation analysis among all samples revealed coefficients of 0.80 between the two HEL samples and 0.83 between the two LEL samples. Furthermore, the two groups were clearly separated by cluster analysis ( Fig. 4 A). Fig. 4 Identifying and annotation of chromatin accessibility variation regions. (A) Spearman correlation of chromatin accessibility across samples. Color intensity and number in each cell indicates the correlation coefficient. (B) Venn diagram of candidate DARs between HEL and LEL groups. Numbers indicate candidate DARs unique to or shared between groups. (C) Genome-wide distribution of candidate DARs. The HEL- and LEL-specific candidate DARs are shown in red and blue, respectively. Each row represents a chromosome; the x-axis indicates genomic position and the y-axis shows peak intensity. (D) Genomic annotation of candidate DARs. The x-axis shows the percentage of peaks in each annotation category. Colored bars represent the proportion of candidate DARs associated with specific genomic features. (E) Distribution of TFBL relative to TSS. The x-axis shows the distance of candidate DARs from the TSS (upstream and downstream). Fig 4 dummy alt text Identifying and annotation of chromatin accessibility variation regions. (A) Spearman correlation of chromatin accessibility across samples. Color intensity and number in each cell indicates the correlation coefficient. (B) Venn diagram of candidate DARs between HEL and LEL groups. Numbers indicate candidate DARs unique to or shared between groups. (C) Genome-wide distribution of candidate DARs. The HEL- and LEL-specific candidate DARs are shown in red and blue, respectively. Each row represents a chromosome; the x-axis indicates genomic position and the y-axis shows peak intensity. (D) Genomic annotation of candidate DARs. The x-axis shows the percentage of peaks in each annotation category. Colored bars represent the proportion of candidate DARs associated with specific genomic features. (E) Distribution of TFBL relative to TSS. The x-axis shows the distance of candidate DARs from the TSS (upstream and downstream). To further explore the impact of ovarian epigenome differences in egg production, we identifying candidate DARs of chromatin between the HEL and LEL groups. A total of 17,014 and 16,377 chromatin accessible regions were detected in the HEL and LEL groups, comprising 1,065 and 428 candidate DARs, respectively ( Fig. 4 B). Candidate DARs were evenly distributed across chromosomes without obvious bias ( Fig. 4 C). Genomic annotation of these candidate DARs ( Fig. 4 D) showed a higher proportion in promoter regions (within 2 kb of a transcription start site ( TSS )) in the HEL group (9.86 %) than in the LEL group (5.84 %). Intronic candidate DARs were also more abundant in the HEL group (10.51 %) than in the LEL group (3.04 %). Furthermore, in the HEL group, minor proportions of candidate DARs were found in exons (0.38 %), UTRs (0.38 %), and 300 bp downstream of TSSs (0.09 %). The remaining candidate DARs in both groups were located in intergenic regions. Analysis of transcription factor-binding loci ( TFBL ) distribution relative to TSSs ( Fig. 4 E) revealed a balanced up- and downstream pattern in the HEL group, contrasting with a downstream-skewed distribution in the LEL group. Moreover, the HEL group also had more TFBLs near TSSs (within 1 kb). To investigate the underlying mechanisms of differential gene expression, we identified DEGs with variation in chromatin accessibility within their promoter regions. This revealed 18 DEGs containing candidate DARs in their promoters, ( Fig. 5 A; Table S7 ). Functional enrichment analysis of the DEGs showed that RFT1, MICA, MHCIA6 and MHCIY were primarily associated with terms “immune response”, “extracellular space” and “integral component of membrane”, and all of them were upregulated in the HEL group. Furthermore, preceding PPI network analysis identified MICA, MHCIA6 and MHCIY among these genes as top 20 % hub genes, so these three were defined as key candidates. Their promoter candidate DARs were shown in Fig. 5 B and Table 2 . Fig. 5 Definition of key genes and prediction of potential TFs. (A) Venn diagram for candidate gene screening. Intersecting between DEG promoters (RNA-seq) and promoter-annotated candidate DARs (ATAC-seq). (B) Chromatin accessibility profiles at key gene promoter regions. Gray-shaded areas highlight the candidate DARs. (C) The TFs predicted to bind promoters of key genes. Fig 5 dummy alt text Table 2 List of the key genes. Table 2 dummy alt text Gene log2FC (RNA-seq) p-value (RNA-seq) Chromosome Peak Start Peak End MICA 4.67 1.05E-04 16 2,068,675 2,068,985 MHCIA6 5.39 1.21E-04 16 2,044,792 2,045,249 MHCIY 3.27 2.13E-03 16 1,915,878 1,916,593 Definition of key genes and prediction of potential TFs. (A) Venn diagram for candidate gene screening. Intersecting between DEG promoters (RNA-seq) and promoter-annotated candidate DARs (ATAC-seq). (B) Chromatin accessibility profiles at key gene promoter regions. Gray-shaded areas highlight the candidate DARs. (C) The TFs predicted to bind promoters of key genes. List of the key genes. To elucidate the mechanistic basis for elevated expression of the three key genes associated with differential chromatin accessibility, we performed GO enrichment analysis on the genes with p < 0.05 and |log₂FC| ≥ 0. This result revealed genes which enriched in term “DNA-binding transcription factor activity, RNA polymerase II-specific (GO:0000981, padj = 0.040)” were considered differentially expressed TFs ( Table S8 ). The candidate DARs within the promoters of the three key genes were screened for potential TF binding sites, and seven TFs with high binding affinity were predicted (relative score > 90 %): HOXB6, E2F1, GATA3, TFAP2A, WT1, NKX3-2 , and KLF10 ( Table S9 ). Given that the TF binding models in JASPAR are primarily based on human or mouse data, protein sequence homology between these TFs and their chicken orthologs was assessed. HOXB6 (79.2 %), GATA3 (92.0 %), TFAP2A (90.8 %), and WT1 (86.4 %) exhibited high conservation, whereas the remaining TFs showed low identity (< 70.0 %) ( Fig. S3 ). Therefore, we proposed that HOXB6, GATA3, TFAP2A , and WT1 are potential regulators upregulating the expression of key genes ( Fig. 5 C). These TFs may enhance the transcriptional initiation of key genes by binding to specific motifs within their promoter regions.

All animal procedures in this work were approved by the Animal Welfare Committee of Huazhong Agricultural University (Hubei, China) (Approval No.: HZAUCH-2024-0044). A total of 872 Wenchang chickens were obtained from China Hainan (Tanniu) Wenchang Chicken Co., Ltd. for subsequent experiments. After 12 weeks of age, the chickens were transferred to single cages within the same coop. They were provided a corn-soybean meal diet formulated according to the Chinese national standard “Nutritional Requirements of Yellow Feathered Broilers” (NY/T 3645 2020) and had free access to water. Once the laying phase began, light was increased weekly by 1 h until reaching 16 h, at which point it was maintained. Daily egg production was recorded manually for each chicken from 44 to 49 weeks of age. Based on these data, three hens with laying rates > 76.0 % were randomly selected to form the high egg-laying ( HEL ) group, and three hens with laying rates < 50.0 % were randomly selected to form the low egg-laying ( LEL ) group. The HEL group consisted of H-1, H-2, and H-3, and the LEL group consisted of L-1, L-2, and L-3. The chickens were euthanized and blood samples were collected. Following abdominal dissection, ovaries were harvested. The egg yolk of the follicles was removed, and remanent tissue was washed with cold PBS. Two parts of ovarian stroma tissue were taken for the subsequent library construction and data analyses. A portion of the tissue was minced, placed in TRIzol reagent (15-596-018, Invitrogen, Carlsbad, CA, USA), frozen in liquid nitrogen and stored at −80°C for total RNA extraction. Another portion of the tissue was frozen in liquid nitrogen and stored at −80°C for ATAC-seq library construction. The remaining tissue was placed in 4 % paraformaldehyde ( PFA ) fix solution (BL539A, Biosharp, Anhui, China) for morphological examination. Ovarian tissue samples from the HEL and LEL groups were fixed in 4 % PFA fix solution for 24 h, followed by overnight rinsing under running water. The samples were then processed through a series of steps including graded ethanol dehydration, clearing, paraffin infiltration, and embedding. Subsequently, using a microtome, 4 μm-thick sections were cut from the paraffin blocks and mounted them onto glass slides. The sections were stained with hematoxylin and eosin ( H&E ), and sealed with neutral balsam. The prepared ovarian histological sections were examined under an upright microscope for observation. Total RNA was extracted from ovarian tissues in TRIzol reagent using the phenol chloroform method. Using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), RNA concentration and quality were assessed. Whereafter, RNA integrity was examined using an Agilent 2100 Bioanalyzer Instrument (Agilent Technologies, Palo Alto, CA, USA). Samples with RIN values > 8.0 were used for mRNA library construction according to the TruSeq mRNA Sample Preparation Kit (Illumina, San Diego, CA, USA), and paired-end sequencing was performed on the Illumina NovaSeq 6000 Sequencing System (Illumina, San Diego, CA, USA). Samples H-1, H-2, and H-3 were from the HEL group, and L-1, L-2, and L-3 were from the LEL group. The raw sequence data have been deposited in the Genome Sequence Archive (GSA) ( Zhang, et al., 2025 ) at the China National Center for Bioinformation (CNCB) web server, with the BioProject ID: PRJC054005. Quality control of the raw sequencing data was performed using FastQC (v0.12.1) ( Wingett and Andrews, 2018 ). Libraries with a Q30 > 80 % and a balanced GC content ranging from 40 % to 60 % were retained for subsequent analyses. The cleaned reads were aligned against the chicken reference genome (galGal6, GCF_000002315.5), and transcript-level quantification was performed using RSEM (v1.3.1) ( Li and Dewey, 2011 ) with the splice-aware STAR aligner (v2.7.10a) ( Dobin, et al., 2012 ). Quantification was run with the parameters: –paired-end –star –calc-ci –strandedness reverse, and gene expression levels were calculated and presented as fragments per kilobase million ( FPKM ). The resulting alignment files were sorted and converted to BAM format using SAMtools (v1.17) ( Danecek, et al., 2021 ) for downstream visualization. Principal component analysis ( PCA ) of the samples was performed using the prcomp function in R (v4.5.1). Using the edgeR package (v4.0.16) ( Robinson, et al., 2010 ), differential expression analysis was conducted. Gene expression was assessed based on fold change ( FC ). Differentially expressed genes ( DEGs ) were identified with thresholds of p < 0.05 and |log₂(FC)| ≥ 1. KOBAS ( Bu, et al., 2021 ) was employed to perform GO and KEGG pathway enrichment analyses on the DEGs with the chicken ( Gallus gallus ) as the reference species. The protein-protein interaction ( PPI ) network of DEGs with an interaction score ≥ 0.40 (moderate confidence) was filtered from the STRING database ( Szklarczyk, et al., 2022 ). The PPI network was imported into Cytoscape software (v3.9.1) ( Shannon, et al., 2003 ) for visualization analysis. For RT-qPCR validation, nine genes were randomly selected from the identified DEGs, with GAPDH serving as the endogenous reference. The specificity of primers was verified through NCBI Primer-BLAST analysis and agarose gel electrophoresis following conventional PCR ( Fig. S2 , Table S10 ). Total RNA was reverse-transcribed into cDNA using the HiScript II Q RT SuperMix for qPCR (+gDNA wiper) kit (R223-01, Vazyme, Jiangsu, China). RT-qPCR was then performed using the ChamQ universal SYBR qPCR master mix (Q711-02, Vazyme, Jiangsu, China) on the ABI real-time PCR system (Applied Biosystems, Carlsbad, CA, USA). Statistical analysis was performed on the ΔCt values (ΔCt = Ct target gene - Ct GAPDH ). Statistical significance between the two groups was determined using an unpaired two-tailed Student's t-test. As only a single comparison was conducted between the two groups, no multiple comparison correction was applied. Frozen tissues were pulverized in liquid nitrogen, and approximately 0.05 g of the powder was processed. Unless otherwise specified, all subsequent steps were performed at 4°C. Samples were resuspended in 1 mL of ice-cold 1 × DPBS solution, then centrifuged at 2500 rcf for 5 min to remove impurities. Nuclei were released after incubation with 1 mL of ice-cold lysis buffer for 3 min. The nuclei were counted under a microscope using a hemocytometer and trypan blue stain. A total of 50,000 nuclei per sample were fragmented with Tn5 transposase (A0215, YINGZI GENE, Hubei, China) at 37 °C for 1 h. The products were purified with a DNA cleanup kit (TD413, GENSTONE BIOTECH, Beijing, China) and quantified by Qubit fluorometry (Thermo Fisher Scientific, Waltham, MA, USA). After PCR amplification, the products were purified, size-selected, and recovered using 0.5 ×, 0.8 ×, and 1.1 × DNA Clean Beads, respectively, to construct the final ATAC-seq libraries. Paired-end sequencing was performed on an Illumina NovaSeq 6000. Samples H-1 and H-2 were from the HEL group, L-1 and L-2 were from the LEL group. RNA-seq was also performed on all of them. The raw sequence data have been deposited in the Genome Sequence Archive (GSA) ( Zhang, et al., 2025 ) at the China National Center f or Bioinformation (CNCB) web server, with the BioProject ID: PRJC054005. Following initial filtering, data quality was assessed using FastQC (v0.12.1) ( Wingett and Andrews, 2018 ). Adapters and low-quality bases were trimmed using Trimmomatic (v0.39) ( Bolger, et al., 2014 ) with the parameters: ILLUMINACLIP:NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36. After trimming, a second round of FastQC was performed to confirm that the retained reads had a Q30 > 80 %, a GC content of 40 % to 60 %, and no remaining adapter contamination. Cleaned paired-end ATAC-seq reads were aligned to the chicken reference genome (galGal6, GCF_000002315.5) using BWA-MEM ( Li, 2013 ). SAM files were converted to coordinate-sorted, indexed BAM files using SAMtools. Reads were then filtered to keep properly paired alignments with MAPQ ≥ 30 and to remove mitochondrial (chrM) reads before downstream analyses. These BAM files were then converted to BigWig format using bamCoverage from deepTools (v3.5.1) ( Ramírez, et al., 2014 ) and normalized by reads per kilobase per million mapped reads ( RPKM ) for visualization in IGV (v2.10.0) ( Thorvaldsdóttir, et al., 2013 ). Genome-wide read distribution was summarized using multiBamSummary, and inter-sample correlations were calculated and visualized with plotCorrelation. To generate representative signal tracks for visualization and aggregated-signal comparison, biological replicates within each condition were merged and indexed using SAMtools ( Danecek, et al., 2021 ). Peak calling was performed with MACS2 (v2.2.8) ( Zhang, et al., 2008 ) in two stages: an initial round using the parameters: –nomodel –shift −37 –extsize 73 -g 1065370000, followed by a second round using the estimated average fragment size with parameters: -q 0.05 -g mm -B –nomodel –extsize 260. Putative differential accessibility regions ( candidate DARs ) were obtained using MACS2 bdgdiff with parameters -g 60 -l 260. We note that this aggregated-signal approach does not explicitly model biological variability across replicates and is therefore used here for exploratory identification of candidate regions. Finally, candidate DARs were annotated using the ChIPseeker (v1.38.0, parameters: tssRegion = c(−2000, 2000), addFlankGeneInfo=TRUE, flankDistance=5000) ( Yu, et al., 2015 ). Genomic coordinates of the DEGs identified were retrieved in batch using the UCSC Genome Browser ( Perez, et al., 2024 ), and a BED file containing their promoter regions was generated. Candidate differential chromatin accessibility peaks overlapping with these promoter regions were then identified using BEDTools (v2.31.0) ( Quinlan and Hall, 2010 ). DEGs whose promoters showed differential accessibility were considered candidate genes. Key genes were subsequently selected by integrating the results of functional enrichment and PPI network analyses. GO enrichment analysis of DEGs (p < 0.05, |log₂(FC)| ≥ 0) was performed using KOBAS ( Bu, et al., 2021 ). The DEGs enriched in the term “DNA-binding transcription factor activity, RNA polymerase II-specific” (GO:0000981, padj = 0.0424; Table S11 ) were defined as differentially expressed transcription factors ( TFs ) ( Table S8 ). TF-binding motifs within the previously identified peak regions ( MICA , chr16:2,068,675-2,068,985; MHCIA6 , chr16:2,044,792-2,045,249; MHCIY , chr16:1,915,878-1,916,593) were analyzed for enrichment using the JASPAR database ( Rauluseviciute, et al., 2023 ). The resulting set of candidate TFs were then aligned to their chicken homologs via UniProt ( The UniProt Consortium, 2025 ), and those with high sequence homology (> 79 %) were selected as potential regulators of the key genes.

In summary, this study reveals that ovarian immune homeostasis plays a critical role in sustaining egg production during the late laying phase in chickens. By integrating RNA-seq and ATAC-seq, we identified key immune-related genes whose expression is associated with changes in chromatin accessibility, and uncovered TFs that may coordinate their regulation. This multi-omics approach enabled the systematic delineation of regulatory networks that single-omics analyses alone could not resolve. Our findings provide potential molecular targets for genetic improvement of laying performance in Wenchang chickens and contribute to a broader understanding of the regulatory mechanisms underlying ovarian function and reproductive aging in poultry.

The late laying phase is a critical period influencing overall egg production in chickens. In this study, we integrated RNA-seq and ATAC-seq to investigate differences in gene expression and chromatin accessibility in the ovaries of high- and low-laying Wenchang chickens during this phase ( Fig. 6 ). We identified three key genes, including MICA, MHCIA6 and MHCIY . All of them belong to the non-classical MHC-Ⅰ (MHC-Ⅰb) family and located within the MHC-Y region of chromosome 16 in chickens. This class of genes is characterized by high polymorphism and allogeneic immunogenicity. They are widely expressed on cell surfaces, and may mediate specialized immune responses ( Miller and Taylor, 2016 ; Zhang, et al., 2020 ). In particular, MHC-Ⅰb molecules are immunoregulators critical for T cell development, differentiation, and overall T cell-mediated immunity ( Kim, et al., 2024 ; Won, et al., 2023 ). Fig. 6 Graphical abstract. The graph was drawn using BioRender ( https://www.biorender.com ). Fig 6 dummy alt text Graphical abstract. The graph was drawn using BioRender ( https://www.biorender.com ). Consistent with previous studies, an association between MHC genes and egg production in chickens has been indicated ( Ouyang, et al., 2000 ). A recent study further suggested that upregulated expression of MHC I gene may enhance both the egg-laying rate and immune function in 55-week-old hens ( Geng, et al., 2026 ). Collectively, our results support the growing consensus that ovarian immune homeostasis plays a critical role in reproductive performance. The ovary contains a variety of immune cells, including monocytes, macrophages, T cells, B cells, and natural killer cells ( Shen, et al., 2023 ). Although these cells constitute only a small fraction of the total ovarian cell population, they establish the necessary immune microenvironment to support cyclic ovulation through cytokine-mediated interactions with the follicular pool ( Cacciottola, et al., 2025 ). With advancing age, immune infiltration in the ovary increases and correlates with upregulation of the MHC‑I pathway. This process may enhance interactions between T cells and neighboring cells, thereby reshaping the immune microenvironment ( Wei, et al., 2026 ). During the late laying phase, sustained ovulation results in the gradual accumulation of inflammatory responses within the ovary ( Zhang, et al., 2022 ). This chronic inflammation is related to reduced ovarian T cell numbers, increased stromal fibrosis, depleted follicular reserve, and elevated follicular atresia rate, collectively affecting ovarian immune function ( Barua and Yoshimura, 1999 ; Wang, et al., 2020 ; Zeng, et al., 2024 ). Our findings suggest that MHC-Ib molecules may help alleviate such chronic inflammation and maintain ovarian immune homeostasis, thereby supporting sustained egg-laying capacity. This interpretation is further supported by our ATAC-seq data, which reveal differences in chromatin accessibility in the promoter regions of these key MHC-Ib genes between high and low egg-laying hens, indicating that their transcriptional regulation is influenced by local chromatin structural changes. In addition to these immune-related genes, our integrated analysis identified four TFs as potential regulators of the key genes, including HOXB6, GATA3, TFAP2A , and WT1 . HOXB6, an ANTP-class homeobox protein, participates in cell differentiation by binding to enhancers ( Yang, et al., 2024 ), and has been implicated in female reproductive system pathologies, including ovarian tumors ( Stone, et al., 2003 ), ovarian cancer ( Kar, et al., 2015 ), and endometriosis ( Geng, et al., 2022 ). GATA3, a zinc-finger transcription factor, is essential for T cell development and inflammatory responses ( Bacha, et al., 2025 ; Popp, et al., 2017 ), and also regulates granulosa cell apoptosis and autophagy in the ovary ( Luo, et al., 2023 ). Additionally, GATA3 can cooperate with TFAP2A to co-regulate chromatin accessibility and target gene expression ( Liu, et al., 2023 ). TFAP2A is a TFAP2 family transcription factor that plays an important role in oocyte maturation regulation ( Lin, et al., 2022 ). WT1, another zinc-finger transcription factor, is involved in various processes such as proliferation, apoptosis, adhesion, and intercellular signaling ( Zhu, et al., 2024 ) and is crucial for granulosa cell development ( Cen, et al., 2020 ). In chickens, WT1 expression is specific to the stage of follicular development, with high expression in immature follicles and sharp downregulation before ovulation ( Chun, et al., 1999 ; Tang, et al., 2022 ). Therefore, our observation of upregulated WT1 expression in the HEL group may indicate the presence of more immature follicles, reflecting greater egg-laying potential. Importantly, our ATAC-seq analysis identified differentially accessible regions in the vicinity of these transcription factor binding sites, suggesting that chromatin remodeling may influence the regulatory networks controlling the expression of the three key immune-related genes. Overall, this multi-omics approach offers a multidimensional perspective for deciphering the regulatory mechanisms underlying variation in egg production during the late laying phase. Nevertheless, several limitations should be acknowledged. Our identification of DEGs was based on a nominal p‑value threshold (p < 0.05) without false discovery rate ( FDR ) correction. Although we used the edgeR package to obtain exact p-values returned from the derived probabilities ( Rapaport, et al., 2013 ), and our identification method was strongly supported by RT-qPCR validation and has been widely adopted in relevant studies ( Chen, et al., 2024a ; Duan, et al., 2026 ; Huang, et al., 2025 ; Liu, et al., 2025 ; Li, et al., 2026 ; Sun, et al., 2024 ; Zhou, et al., 2024 ), we recognize that it may increase the risk of false‑positive findings. Differences in chromatin accessibility among ovaries from different laying groups were identified only on an exploratory basis, without considering the biological variability among replicate samples. In addition, the key genes and potential TFs identified in this study are supported only by correlative evidence and lack functional validation. Future studies involving more samples should prioritize applying stricter criteria (e.g., FDR) for DEG identification to obtain more conservative estimates, and should employ replicate‑aware statistical frameworks (e.g., DESeq2 or edgeR) to validate these candidate DARs and estimate true biological variance. For the key genes and potential TFs, chromatin immunoprecipitation and electrophoretic mobility shift assays can be employed to confirm their interaction, and gain- or loss-of-function experiments can be conducted to elucidate the roles of the key genes in downstream regulatory pathways. Furthermore, integrating other epigenomic and three‑dimensional genomic technologies will enable a more comprehensive understanding of ovarian regulatory mechanisms during the late laying phase, thereby providing reliable targets for precision breeding of Wenchang chickens.

The egg-laying cycle of hens comprises three distinct stages: the early, peak and late laying phases. Notably, the late laying phase constitutes over 50 % of the total cycle duration. The egg-laying rate in this phase constrains total egg production and represents a critical focus for enhancing overall production. However, during the late laying phase, the egg-laying rate exhibits a progressive decline with aging in chickens, particularly in broiler breeds ( Zhang et al., 2021 ). A previous study showed that the egg-laying rate in this phase has strong genetic characteristics and population-level variation. In a broiler hybrid line, the heritability and CV of egg-laying rate at 52 weeks of age are 0.426 and 33.3 %, respectively ( Chomchuen, et al., 2022 ). GWAS has also identified several SNPs and candidate genes that influence egg-laying traits in chickens, such as SHROOM2, SYNE3 , and CNNM2 ( Tan, et al., 2025 ). Therefore, genetic investigations focusing on this phase are essential for promoting the improvement of egg production performance in Wenchang chickens. The ovary is a critical component of the avian reproductive system, with its physiological status playing an important role in modulating reproductive traits. The chicken ovary comprises follicles and ovarian stroma. The latter is essential for maintaining normal follicle development, and is stratified into cortex and medulla. The cortex, constituting the outer layer, is the site of follicular growth and maturation. The medulla, located on the inner side, consists of blood vessels, nerves, lymphatic vessels, immune cells, the extracellular matrix, and other ovarian-specific components. Together, these elements generate a complex microenvironment, and perform critical functions in follicle development, ovulation regulation, and maintenance of ovarian immune homeostasis ( Kinnear, et al., 2020 ). In the late laying phase, ovarian aging exacerbates oxidative stress and local inflammation, leading to follicular atresia, granulosa cell apoptosis, and ultimately impairing egg production ( Hao, et al., 2020 ; Zeng, et al., 2024 ). Multiple studies have indicated that genes such as LECT2, BMP4 , and COL1A1 can facilitate granulosa cell proliferation and ensure normal follicular development by modulating inflammatory responses within the ovarian stroma, alleviating endoplasmic reticulum stress, and promoting cell cycle progression ( Chen, et al., 2024b ; Li, et al., 2025 ; Yao, et al., 2020 ). Although this evidence suggests that ovarian decline is a key factor contributing to reduced reproductive capacity, the molecular mechanisms underlying the decline in egg production in chickens during this process remain to be fully elucidated. At present, RNA-seq is widely used in research on egg production traits in poultry, and has been used to identify many key genes and pathways, including GRIA1, PRLR, THBS2 , and IGF1 , as well as pathways such as neuroactive ligand-receptor interaction, ECM-receptor interaction, immune function, MAPK signaling, and calcium signaling ( Bhavana, et al., 2022 ; He, et al., 2022 ; Mu, et al., 2021 ; Zhang, et al., 2019 ). However, a single-omics approach cannot elucidate the regulatory mechanisms underlying ovarian decline during the late laying phase. Differences in gene expression often arise from alterations in the activity of upstream regulatory elements, and chromatin accessibility is a key determinant of their functional capacity. ATAC-seq, which maps genome-wide chromatin accessibility, identifies active regulatory elements such as promoters and enhancers ( Buenrostro, et al., 2013 ). Integrating ATAC-seq with RNA-seq provides a powerful strategy to clarify gene regulatory mechanisms by correlating chromatin accessibility with gene transcriptional levels. This multi-omics approach has been increasingly used to elucidate the genetic regulatory mechanisms of complex traits such as reproduction and organ development. Zhang et al. (2023) identified genes such as FGFR3 related to gonadal development in sex-reversed chickens. Li et al. (2022) demonstrated that transcriptional activity in chicken granulosa cells is regulated by dynamic changes in chromatin structure. Dinh et al. (2023) demonstrated that ovulation stimulation significantly influences cellular chromatin accessibility in mice. This multi-omics integration thus provides an effective approach for elucidating regulatory networks within the ovaries of late-laying hens. The Wenchang chicken, a prominent yellow-feathered breed in China, is valued for its desirable meat quality. However, this breed exhibits poor egg-laying performance, which adversely affects the economic interests of breeding enterprises. Consequently, investigating how to enhance the egg production of Wenchang chickens holds significant importance for the poultry industry. Furthermore, its relatively short domestication history has retained substantial genetic diversity. This is manifested as significant phenotypic variation in egg production, and reflects considerable selective breeding potential ( Ren, et al., 2024 ; Xing, et al., 2020 ). Therefore, Wenchang chicken is a valuable model for elucidating mechanisms by which ovarian differences regulate laying rates. In summary, the molecular mechanisms by which the ovarian stroma contributes to the decline in egg production during the late laying phase remain poorly understood, especially in indigenous breeds such as the Wenchang chicken. In this study, we integrated RNA-seq and ATAC-seq analyses to investigate the mechanisms underlying variation in egg production in Wenchang chickens during this phase. We characterized differences in gene expression and chromatin accessibility in the ovarian stroma between chickens with distinct egg-laying rates, identified three key regulatory genes, and explored the associated cis- and trans-regulatory mechanisms. By elucidating these molecular networks, this study provides a theoretical foundation and identifies potential molecular targets for the genetic improvement of egg-laying performance in Wenchang chickens.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.