Transcriptome and proteome profile analysis of the regulation of chicken ovarian development.

In this study, the egg production rate of TBSF increased most rapidly from week 20 to week 30, rising quickly from approximately 2 % to 70 %, then stabilizing ( Fig. 1 A). The follicles in the ovaries were isolated and arranged one by one. The ovarian morphology and follicle characteristics at week 20 and week 30 showed that the ovaries were intact, containing follicles at various developmental stages ( Fig. 1 B). Statistical analysis of follicle numbers revealed that SWF, LWF, and SYF counts at 30 weeks were significantly higher than those at 20 weeks ( Fig. 1 C). HE staining results indicated that ovarian tissue at 20 weeks was more compact, whereas at 30 weeks it appeared looser ( Fig. 1 D). The number of follicles within the ovary was significantly higher at 20 weeks compared to 30 weeks ( Fig. 1 E). Serum hormone analysis showed that FSH, LH, E2, and P4 levels at 30 weeks were higher than those at 20 weeks ( Fig. 1 F–I). Fig. 1 Analysis of egg production performance, ovarian morphology, follicle counts, serum hormone levels, and gene expression. A. Laying performance statistics from week 20 to week 50. B. Ovarian morphology and follicles at 20 and 30 weeks. C. The number of follicles at different stages in the ovary ( n = 4). D. HE section of the ovary from TBSF. E. Statistical analysis of follicle numbers within the ovary. F. Follicle-stimulating hormone (FSH). G. Luteinizing hormone (LH); H. Estradiol (E2); I. Progesterone (P4). * P  < 0.05 and ** P  < 0.01. Fig. 1 Analysis of egg production performance, ovarian morphology, follicle counts, serum hormone levels, and gene expression. A. Laying performance statistics from week 20 to week 50. B. Ovarian morphology and follicles at 20 and 30 weeks. C. The number of follicles at different stages in the ovary ( n = 4). D. HE section of the ovary from TBSF. E. Statistical analysis of follicle numbers within the ovary. F. Follicle-stimulating hormone (FSH). G. Luteinizing hormone (LH); H. Estradiol (E2); I. Progesterone (P4). * P  < 0.05 and ** P  < 0.01. In this study, we sequenced 8 cDNA libraries on the Illumina NovaSeq PE150 platform to investigate transcript differences between ovarian tissues at two different time points, with 4 RNA-Seq samples in each group. Each library generated between 40,485,546 and 49,564,064 raw reads, and the number of filtered clean reads ranged from 42,273,598 to 42,045,430. These clean reads were then aligned to the reference genome GRCg7b (GCA_016699485.1), resulting in alignment rates between 85.37 % and 90.91 %. The Q20 content exceeded 95 %, and the Q30 content exceeded 94 % ( Table 1 ). The sample correlations and gene expression distributions met standard quality control criteria. ( Supplementary figure 1A-B ). PCA analysis of the samples revealed that the two groups could be distinctly separated into two clusters ( Fig. 1 J), and the clustering heatmap showed significant differences in gene expression patterns between the two groups ( Fig. 1 K). Based on these metrics, the transcriptome data were used for subsequent functional enrichment and differential expression analyses. Table 1 RNA-seq raw sequencing data reads. Table 1 Sample Raw_reads Clean_reads Total_map Unique_map Q20 Q30 20W 44301628 40072716 36884933(92.05 %) 36428884(90.91 %) 98.08 95.16 48624928 46527004 42125379(90.54 %) 41468301(89.13 %) 97.85 94.7 44732052 42779438 37513331(87.69 %) 36518778(85.37 %) 97.57 94.15 41742154 39715234 35839385(90.24 %) 35262971(88.79 %) 97.8 94.59 30W 46274190 44198162 40138374(90.81 %) 39560164(89.51 %) 97.87 94.68 40485546 38110844 35026237(91.91 %) 34507349(90.54 %) 98.01 94.97 42432940 38806540 35600341(91.74 %) 35105267(90.46 %) 98.06 95.16 49564064 47066174 43287874(91.97 %) 42587850(90.49 %) 98.1 95.21 RNA-seq raw sequencing data reads. The transcriptome results of this study detected a total of 4563 differentially expressed genes, of which 1677 were upregulated and 2886 were downregulated ( Fig. 1 L, Supplementary figure 1C and Table S3 ). Subsequently, we performed GO and KEGG enrichment analyses on these differentially expressed genes. KEGG showed that the differentially expressed genes were significantly enriched in pathways such as Ribosome, Cytokine-cytokine receptor interaction, ECM-receptor interaction, Cell cycle, Cell adhesion molecules, Steroid biosynthesis, Oocyte meiosis, and Progesterone-mediated oocyte maturation. ( Fig. 2 A and Table S4 ). GO enrichment analysis Biological Process ( BP ) mainly indicated significant enrichment in pathways related to cell development, extracellular matrix structure, inflammatory response, and angiogenesis. The main enrichment of molecular function ( MF ) includes binding and catalytic activities, extracellular matrix structural constituents, cytokine receptor activity, and collagen binding, among others. Cellular components ( CC ) are mainly enriched in ribosomal subunits, extracellular matrix components, membrane regions, cell junctions, and other cellular structures ( Fig. 1 F-H and Table S5 ). GSEA enrichment analysis identified significant upregulation of pathways such as ECM-receptor interaction, focal adhesion, and cell adhesion molecules in 20-week-old ovarian tissues. ( Fig. 2 B-D). The steroid biosynthesis pathway was observed to be significantly upregulated in ovarian tissues at 30 weeks, corresponding to the peak laying period ( Fig. 2 E). Transcriptome sequencing analysis identified significant enrichment of pathways related to extracellular matrix synthesis and cell proliferation in early ovarian development stages, suggesting potential roles in the developmental processes leading up to the onset of the peak laying period. Moreover, the high expression levels of genes involved in the steroid biosynthesis pathway were observed during the peak laying period, indicating a potential association with hormonal regulation during this stage. Fig. 2 Functional annotation analysis of DEGs. A. Dot plot of pathway enrichment analysis, where the size of each dot represents the number of genes involved in a specific pathway, and the color intensity indicates the P -value, with darker colors representing higher statistical significance. B. Cell adhesion molecules. C. ECM-receptor interaction. D. Focal adhesion. E. Steroid biosynthesis. Normalized enrichment score ( NES ) is shown for each pathway, indicating the direction and magnitude of enrichment (NES > 0: upregulated in 30 W group; NES < 0: downregulated in 30 W group). F-H. Bar plot of GO enrichment analysis, where the length of the bars represents the significance of enrichment, and the curve shows the number of enriched genes. Fig. 2 Functional annotation analysis of DEGs. A. Dot plot of pathway enrichment analysis, where the size of each dot represents the number of genes involved in a specific pathway, and the color intensity indicates the P -value, with darker colors representing higher statistical significance. B. Cell adhesion molecules. C. ECM-receptor interaction. D. Focal adhesion. E. Steroid biosynthesis. Normalized enrichment score ( NES ) is shown for each pathway, indicating the direction and magnitude of enrichment (NES > 0: upregulated in 30 W group; NES < 0: downregulated in 30 W group). F-H. Bar plot of GO enrichment analysis, where the length of the bars represents the significance of enrichment, and the curve shows the number of enriched genes. We employed qRT-PCR to validate the results of RNA-seq, and 11 selected DEGs were used for verification. The results of qRT-PCR are similar to the expression differences obtained by RNA-seq ( Fig. 3 A). The linear regression between the RNA-seq and qRT-PCR results shows a positive correlation, with a correlation coefficient (R) of 0.87, which supports the reliability of the RNA-seq results ( Fig. 3 B). Fig. 3 qRT-PCR validation of RNA-seq results. A. Pearson correlation analysis of qRT-PCR and RNA-seq DEGs. B. Comparison of fold changes of 11 DEGs between qRT-PCR and RNA-seq. Fig. 3 qRT-PCR validation of RNA-seq results. A. Pearson correlation analysis of qRT-PCR and RNA-seq DEGs. B. Comparison of fold changes of 11 DEGs between qRT-PCR and RNA-seq. Transcriptome data provide insights into gene expression in ovarian tissues, but they represent only an intermediate stage in understanding protein-level regulation and functional outcomes. It cannot show the specific expression at the protein level. To better understand the differences in protein abundance between 30 W and 20 W, this study conducted a proteomics study based on Tandem Mass Tag ( TMT ). Quality control showed the cumulative distribution of the coefficient of variation ( CV ) of protein expression in the proteome, with samples from both groups exhibiting good expression stability ( Fig. 4 A). PCA analysis results indicated that the proteins in the two groups clustered well together ( Fig. 4 B). Overall, the protein results were suitable for subsequent analysis. The volcano plot displayed differential expressed proteins between the two groups, and the clustering heat map showed significant differences in protein expression patterns between the two groups. In this study, a total of 6212 proteins were identified in the proteome. Comparative analysis revealed 154 differentially expressed proteins, among which 29 were up-regulated and 125 were down-regulated ( Fig. 4 C-D and Table S6 ). To further understand the functions of the differentially expressed proteins, GO and KEGG enrichment analyses were performed on all differentially expressed proteins. KEGG showed that differentially expressed proteins were significantly enriched in nucleotide metabolism, ABC transporters, arginine and proline metabolism, and focal adhesion signal pathways ( Fig. 4 E Table S7 ). GO enrichment showed that signal pathways related to biological processes were mainly concentrated in extracellular structure, tissue development, and signal transduction ( Fig. 4 F Table S8 ). In the GO enriched BP terms, proteins closely related to tissue development included the annexin family proteins and matrix-producing proteins ( Fig. 4 G). The annexin protein family members ANXA1, ANXA2, ANXA5, ANXA7, along with the matrix development proteins ITGB3, COL5A2, COL1A2, and POSTN, were highly expressed at 20 weeks. Interestingly, at 30 weeks, these proteins are associated with higher expression levels of the matrix production protein YES1 and the antioxidant protein ALDH7A1 ( Fig. 5 A). To understand the interaction relationships between proteins, this study used the STRING database for protein interaction analysis and found that the proteins in the interaction network were mainly concentrated in matrix development and cell proliferation ( Fig. 5 B). Further analysis of the proteomics data revealed that pathways related to cell matrix and tissue development in the ovary may be crucial in the early development of ovarian tissues. To validate the proteomic results, we selected three proteins, Annexin A5 ( ANXA5 ), Cytochrome P450 Family 17 Subfamily A Member 1 ( CYP17A1 ), and Aldolase and Yes Proto-Oncogene 1 ( YES1 ) for verification. The results of Western blotting were consistent with the trends observed as those of the quantitative proteomic analysis ( Fig. 5 C-D). Fig. 4 30 W vs. 20 W Proteome quality control and functional notes. A. Cumulative coefficient of variation (CV) curve, with the X-axis representing the coefficient of variation (CV) and the Y-axis representing the cumulative frequency. The curve compares the stability and variability of protein expression between the two sample groups; B. Principal component analysis (PCA) scatter plot, based on protein expression data, showing the overall differences in proteome expression between the two sample groups; C. Volcano plot, with the X-axis representing log2 fold change and the Y-axis representing -log10 P-value. The color represents the significance levels of the P-values, ranging from blue (not significant) to yellow and red (most significant); D. Heatmap and clustering analysis. D. Shows the expression level differences of specific genes through a heatmap, where each row represents a gene and each column represents a sample. The color intensity indicates gene expression levels, from blue (low expression) to red (high expression). Fig. 4 Fig. 5 Protein interaction networks and Western bloting validation for proteomics. A. Heat map of DAPs involved in ovarian development. B. Protein interaction network diagram. C, D: The protein expression levels of ANXA5, YES1, and CYP17A1 were compared between the W30 and W20 groups. Statistically significant differences are indicated by * P  < 0.05 and ** P  < 0.01. Fig. 5 30 W vs. 20 W Proteome quality control and functional notes. A. Cumulative coefficient of variation (CV) curve, with the X-axis representing the coefficient of variation (CV) and the Y-axis representing the cumulative frequency. The curve compares the stability and variability of protein expression between the two sample groups; B. Principal component analysis (PCA) scatter plot, based on protein expression data, showing the overall differences in proteome expression between the two sample groups; C. Volcano plot, with the X-axis representing log2 fold change and the Y-axis representing -log10 P-value. The color represents the significance levels of the P-values, ranging from blue (not significant) to yellow and red (most significant); D. Heatmap and clustering analysis. D. Shows the expression level differences of specific genes through a heatmap, where each row represents a gene and each column represents a sample. The color intensity indicates gene expression levels, from blue (low expression) to red (high expression). Protein interaction networks and Western bloting validation for proteomics. A. Heat map of DAPs involved in ovarian development. B. Protein interaction network diagram. C, D: The protein expression levels of ANXA5, YES1, and CYP17A1 were compared between the W30 and W20 groups. Statistically significant differences are indicated by * P  < 0.05 and ** P  < 0.01. To further understand the key DEGs and DAPs in the early development of chicken ovaries, an integrated analysis of transcriptome and proteome data was performed. The nine-quadrant plot illustrated the distribution profile of genes and proteins, showing that few genes and proteins had the same trend ( Fig. 6 A). Correlation analysis also revealed a low correlation between genes and proteins ( Fig. 6 B). We screened for co-expressed genes and proteins with the same expression trend and significant differences, identifying 22 overlapping genes ( Fig. 6 C). Identified significant differentially expressed genes from both transcriptome and proteome were subjected to protein interaction analysis, resulting in the discovery of interaction networks associated with tissue development ( Fig. 6 D). Among these, 21 genes were upregulated in the 20 W group ( STAB1, AEBP1, COL12A1, CPXM1, TNFAIP8L3, PALMD, CPN1, LGMN, SERPINF2, OSBP2, FBLN1, SELENBP1, COL1A2, ANXA5, ANXA2, S100A6, CRISPLD2, NEU3, TAL1, IRAK4, EDNRA ), and the YES1 gene was downregulated in the 20 W group ( Fig. 6 E). Integrated transcriptome and proteome analyses identified 22 co-expressed genes and proteins that exhibited significant differential expression between 20 and 30 weeks, suggesting potential involvement in ovarian development. Fig. 6 Integrated analysis of transcriptome and proteome. A. Nine-quadrant plot of transcriptome and proteome, using log2FC to compare the same expression trends of genes and proteins between 30 weeks and 20 weeks. B. Correlation plot between transcriptome and proteome, showing the correlation between transcriptome and proteome data. The x-axis represents log2FC from RNA-seq, and the y-axis represents log2FC from the proteome. The linear regression line and correlation coefficient (r) along with its significance (P-value) in the plot indicate the statistical correlation between the two. C. Venn diagram of transcriptome and proteome, displaying the number of genes with increased and decreased expression in the transcriptome and proteome and their intersections. D. Protein interaction network plot of overlapping genes in the transcriptome and proteome. E. Heatmap of co-expressed genes in transcriptome and proteome. The figure shows the heatmap of co-expressed genes in the transcriptome and proteome at 30 weeks and 20 weeks. Fig. 6 Integrated analysis of transcriptome and proteome. A. Nine-quadrant plot of transcriptome and proteome, using log2FC to compare the same expression trends of genes and proteins between 30 weeks and 20 weeks. B. Correlation plot between transcriptome and proteome, showing the correlation between transcriptome and proteome data. The x-axis represents log2FC from RNA-seq, and the y-axis represents log2FC from the proteome. The linear regression line and correlation coefficient (r) along with its significance (P-value) in the plot indicate the statistical correlation between the two. C. Venn diagram of transcriptome and proteome, displaying the number of genes with increased and decreased expression in the transcriptome and proteome and their intersections. D. Protein interaction network plot of overlapping genes in the transcriptome and proteome. E. Heatmap of co-expressed genes in transcriptome and proteome. The figure shows the heatmap of co-expressed genes in the transcriptome and proteome at 30 weeks and 20 weeks. Correlation analysis of all genes and proteins revealed a positive correlation (Spearman coefficient > 0, P < 0.05) in 181 mRNA-protein pairs, and a negative correlation (Spearman coefficient < 0, P < 0.05) in 80 mRNA-protein pairs ( Fig. 7 A). The positively correlated mRNA-protein pairs were enriched in ECM-receptor interaction, focal adhesion, and steroid biosynthesis ( Fig. 7 B). In the early stages of egg production, ITGA8, ITGA7, and ITGB1, which are associated with extracellular matrix development, showed significant differences in both transcriptomic and proteomic analyses ( Fig. 7 C-E). During the peak period of egg production, CYP17A1 and HSD17B7, which are involved in steroid biosynthesis, and YES1, which is associated with cell proliferation and development, also exhibited significant differences ( Fig. 7 F-H). Fig. 7 Genes and proteins correlation analyses. A. Histogram of correlation analysis for all mRNA-protein. B. KEGG enrichment for mRNA-protein with significant correlations. C—H. Analysis of gene and protein expression profiles and correlation coefficients at 30 weeks and 20 weeks, C(ITGA8), D(ITGA7), E(ITGB1), F(CYP17A1), G(HSD17B7), H(YES1). Fig. 7 Genes and proteins correlation analyses. A. Histogram of correlation analysis for all mRNA-protein. B. KEGG enrichment for mRNA-protein with significant correlations. C—H. Analysis of gene and protein expression profiles and correlation coefficients at 30 weeks and 20 weeks, C(ITGA8), D(ITGA7), E(ITGB1), F(CYP17A1), G(HSD17B7), H(YES1).

All chickens in this study were sourced from Xichang Fengxiang Poultry Industry Co., Ltd. (Taihe, China). The processes of hatching, brooding, rearing, and egg-laying of TBSFs took place under the auspices of this company. Operations were conducted in strict adherence to the disease purification and feeding program developed by the company, ensuring scientific and uniform management. For the purposes of this study, 1200 TBSFs were relocated to the same laying hen house after reaching the age of 100 days, where they were housed in individual cages. The statistics of egg-laying performance were carried out in accordance with what was done in our previous research ( Huang et al., 2025 ). This study was approved by the Animal Welfare Committee of Zhejiang University (Approval No.: ZJU20190149). At 20 weeks of age, four laying hens with uniform body size were selected. Similarly, at 30 weeks of age, four laying hens with uniform body size were randomly selected. All samples were fasted for 12 hours before sampling. Bloodletting was performed using cervical dislocation to minimize the hens' suffering. The abdominal cavity was carefully opened to extract the ovaries, removing atretic follicles and post-ovulatory follicles, and the follicles were then counted and classified into SWF, LWF, SYF, LYF, and preovulatory follicles. After removing the preovulatory follicles, the remaining ovarian tissue was washed with pre-cooled PBS, minced, placed into test tubes, mixed, and then quickly immersed in liquid nitrogen for subsequent experiments. The levels of follicle-stimulating hormone ( FSH ), luteinizing hormone ( LH ), estradiol ( E2 ), and progesterone ( P4 ) in serum were measured using an ELISA kit (Jiangsu Meimian Industrial Co., Ltd., Jiangsu, China) and a microplate reader. Catalog numbers: FSH: 627; LH: 12062; E2: 627; P4: 9556. According to the method of Zhang et al. (2025) , four ovaries from each of the two groups were fixed in 4 % paraformaldehyde (G1101-500ML, Servicebio, Wuhan, China), embedded in paraffin blocks, then sectioned, stained with hematoxylin and eosin, and examined microscopically. Samples quality control and the preparation of libraries for transcriptome sequencing were initiated with an RNA integrity evaluation utilizing the RNA Nano 6000 Assay Kit via the Bioanalyzer 2100 system (Agilent Technologies, CA, USA). Isolation of mRNA from total RNA was accomplished with poly-T oligo-attached magnetic beads, which was then followed by mRNA fragmentation using divalent cations at elevated temperatures and the synthesis of first strand cDNA via random hexamer primers and M-MuLV Reverse Transcriptase, lacking RNase H activity. The second strand cDNA was synthesized using DNA Polymerase I and RNase H, with blunt ends achieved through exonuclease/polymerase activity. Following adenylation of the DNA fragment ends, adaptors featuring hairpin loop structures were ligated, and cDNA fragments ranging from 370∼420 bp were selected via the AMPure XP system. The amplification process utilized Phusion High-Fidelity DNA polymerase, universal PCR primers, and an Index Primer, which was followed by the purification and library quality assessment using the Agilent Bioanalyzer 2100 system. The clustering of index-coded samples was executed on a cBot Cluster Generation System employing the TruSeq PE Cluster Kit v3-cBot-HS (Illumina, San Diego, CA, USA), setting the stage for sequencing on the Illumina NovaSeq platform (Illumina, San Diego, CA, USA), which generated 150 bp paired-end reads. Quality control of raw fastq data was conducted using fastp software, which generated clean reads by removing adapters, reads containing poly-N, and low-quality reads. Additionally, metrics such as Q20, Q30, and GC content of the clean data were calculated, forming the basis for all subsequent analyses. For reads mapping, the reference genome and gene model annotation files were obtained directly from the genome database. The index of the reference genome was constructed with Hisat2 v2.0.5, and paired-end clean reads were aligned to the reference genome using the same tool. Gene expression levels were quantified using featureCounts v1.5.0-p3, with Fragments Per Kilobase of exon per Million fragments mapped ( FPKM ) calculated to consider both sequencing depth and gene length, making it a widely adopted method for estimating gene expression. Differential expression analysis between two groups was carried out using the DESeq2 R package, which uses a negative binomial distribution model to identify differentially expressed genes, adjusting P-values with the Benjamini and Hochberg method ( Benjamini and Hochberg, 1995 ), controling the false discovery rate. Genes with an adjusted P-value of ≤ 0.01 and |log2 fold change| ≥ 1 were identified as differentially expressed. We conducted Gene Ontology ( GO ) and Kyoto Encyclopedia of Genes and Genomes ( KEGG ) enrichment analyses of differentially expressed genes using the R software package( Yu et al., 2012 ), where GO terms and KEGG pathways with corrected P-values less than 0.05 were considered significantly enriched. Additionally, Gene Set Enrichment Analysis ( GSEA ) enrichment analysis was performed for the entire gene set using the same software package. The chicken ovary tissue samples were retrieved from a −80°C freezer, ground into powder at low temperatures, and fully lysed in a centrifuge tube. The lysate was centrifuged at 12,000 g for 15 minutes, and the supernatant was collected. Dithiothreitol ( DTT ) was added, and the mixture was incubated at 56°C for 1 hour. After the reaction was complete, Iodoacetamide ( IAM ) was added and incubated in the dark for another hour. Four volumes of pre-cooled acetone were added, and the mixture was precipitated at −20°C for 3 hours. The precipitate was collected by centrifuging at 12,000 g for 15 minutes, washed with 1 mL of pre-cooled acetone, centrifuged again to collect the precipitate, air-dried, and then dissolved in a protein dissolving solution. A 120 µg aliquot of each protein sample was taken, and protein dissolving solution was added to bring the volume to 100 µL. After adding trypsin and Triethylammonium Bicarbonate ( TEAB ) buffer, the mixture was well-mixed and digested overnight at 37°C. An equal volume of 1 % formic acid was added, mixed well, and centrifuged at 12,000 g for 5 minutes. The supernatant was collected and desalted using a Octadecyl silane ( C18 ) desalting column. The column was washed three times with 1 mL washing solution, then eluted three times with 0.4 mL elution solution. The eluted samples were combined and lyophilized. They were then reconstituted in 100 µL of 0.1 M TEAB buffer, 41 µL of TMT labeling reagent dissolved in acetonitrile was added, and the mixture was incubated at room temperature for 2 hours. Then, 8 % ammonia solution was added to terminate the reaction. Equal volumes of the labeled samples were combined, desalted, and lyophilized. Fractionation was performed using an HPLC system (Thermo Fisher). Proteome Discoverer 2.2 software was used to obtain the relative quantitative values of peptide spectrum matches ( PSMs ) for each sample based on the peak areas of the original spectra. The relative quantitative values of unique peptides were then corrected based on the quantitative information of all PSMs contained in the identified unique peptides. Finally, the relative quantitative values of each protein were corrected based on the quantitative information of all unique peptides contained in each protein. A T-test was conducted on the relative quantitative values of each protein in the two comparative samples. Proteins were considered upregulated when the fold change ( FC ) was ≥1.2 and P-value ≤ 0.05, and downregulated when the FC was ≤0.83 and P-value ≤ 0.05. Annotating proteins using the GO ( www.geneontology.org ) and KEGG ( http://www.genome.jp/kegg/ ) databases helps to further understand the biological processes and functions they are involved in. Protein interactions were analyzed using the STRING database, and network diagrams were drawn using Cytoscape 3.10.1. To validate the results of RNA-seq, 11 differentially expressed genes were selected for analysis. The validation was carried out by real-time quantitative PCR ( qPCR ) using the SYBR Premix PCR Kit (Sangon Biotech, Shanghai, China) on the CFX96 Touch Real-Time PCR Detection System (Bio-Rad, USA). Primers were designed using the NCBI primer design tool according to the chicken mRNA sequences in GenBank, the primer information is shown in the Table S1 . Each reaction was performed in triplicate. The 2 -ΔΔCt method was used for the statistical analysis of qPCR data. The relative expression levels of proteins in ovarian tissue were detected using Western blotting. Total protein from ovarian tissue were extracted using RIPA lysis buffer mixed with protease inhibitors, and the supernatant concentration was determined using a BCA protein assay kit (Beyotime Biotechnology, Jiangsu, China). Proteins were separated by Bis–Tris SDS-PAGE. A semi-dry transfer system (Bio-Rad, USA) was used to transfer the proteins from the gel to a polyvinylidene fluoride ( PVDF ) membrane, ensuring efficient and uniform transfer. The PVDF membrane was blocked with 5 % skim milk at room temperature for 1 hour, followed by overnight incubation with the primary antibody at 4°C. Finally, the membrane was incubated with horseradish peroxidase-conjugated secondary antibodies (HuaBio, Hangzhou, China) at room temperature for 1 hour. To ensure accurate protein quantification, internal reference proteins were selected based on their molecular weight similarity to the target proteins. Specifically, GAPDH was used as the loading control for YES1 and CYP17A1 due to their comparable molecular weights to α-Tubulin, whereas α-Tubulin served as the reference for ANXA5, which has a molecular weight closer to that of GAPDH. This approach ensured clear separation of bands and minimized potential variability caused by molecular weight discrepancies. Protein bands were quantified using Image Lab 6.1 software (Bio-Rad, USA). Information on primary antibodies and dilution ratios is provided in Table S2 . All experiments were performed in triplicate, and data are expressed as the mean ± standard error of the mean ( SEM ). Differences between groups were analyzed using an independent t-test. P < 0.05 was considered statistically significant.

In this study, we investigated ovarian tissues at the early egg-laying stage and peak egg-laying period through morphological observation, serum biochemical analysis, transcriptomic profiling, and proteomic characterization. Histological examination revealed that ovarian tissue was more compact at the early stage but exhibited a looser structure during the peak laying period. Furthermore, reproductive hormone levels were significantly elevated during the peak egg-laying phase. We identified potential pathways such as cell adhesion molecules, ECM-receptor interaction, focal adhesion, and steroid biosynthesis, along with candidate genes including COL12A1, COL1A2, ANXA2, ANXA5, OSBP2, LGMN, EDNRA, CRISPLD2, SERPINF2, CYP17A1, YES1 and HSD17B1 . These may be related to early ovarian development after sexual maturation and the maintenance of high egg production rates. Further research is recommended to investigate the mechanism of action of the above pathways and candidate genes on ovarian development and egg production.

Evaluating egg production performance is a crucial economic trait indicator for poultry. By studying the ovarian development process in chickens, important genes related to ovarian development and key pathways involved in egg production can be identified ( Zhou et al., 2020 ). This study integrates transcriptome and proteome data to reveal the molecular mechanisms during the early ovarian development in chickens, identifying key candidate genes and biological processes associated with early ovarian development in hens. The results of this study elucidate the multi-omics expression profiles and potential regulatory networks of early ovarian development in TBSF. Follicle number is an important indicator of egg production performance, and follicle selection is a critical process during laying ( Zhang et al., 2023 ). In this study, the total follicle count at 30 W was higher than at 20 W, whereas the number of SYF was lower at 30 W compared to 20 W. This suggests that follicle selection intensifies during peak laying, accelerating the transition of SYF to LYF to support sustained ovulation. Ovarian morphology showed clear structural differences between 20 W and 30 W. HE staining revealed that ovarian tissue at 20 W was more compact, while that at 30 W appeared looser. More follicles were observed within the ovarian tissue at 20 W than at 30 W. It is speculated that many follicles had not yet protruded from the ovarian surface during the early laying stage, whereas most had matured and emerged by 30 W, leading to fewer internal follicles and a looser tissue structure ( Zhang et al., 2025 ). Hormones are key regulators of ovulation in poultry, with FSH, LH, E2, and P4 playing coordinated roles in maintaining ovarian development and ovulatory dynamics ( Jayaraman and Kumar, 2017 , Lin et al., 2011 ). In this study, levels of FSH, LH, E2, and P4 were higher at 30 W than at 20 W, indicating increased hormone synthesis and ovarian demand at this stage. These hormones likely contribute to the precise regulation required to maintain hormonal homeostasis during the peak laying period. In this study, signaling pathways and biological processes associated with collagen synthesis and cell adhesion were found to be significantly enriched and upregulated during the early stages. Collagen is a major component of the extracellular matrix ( ECM ) and plays a crucial role in regulating various cellular functions. The ECM not only provides structural support for the normal physiological activities of tissue cells but also plays an indispensable role in immune regulation under both homeostatic and pathological conditions ( Conway and Jacquemet, 2019 ). Cell adhesion molecules ( CAMs ) are proteins that facilitate cell adhesion to the ECM, contributing to cell signaling, immune responses, and developmental processes. CAMs are essential for promoting tight junctions, gap junctions, and other forms of cell-cell interactions, forming the basis of tissue integrity and communication. Focal adhesions are large dynamic protein complexes that anchor the cell cytoskeleton to the ECM ( Hu et al., 2023 ). This interaction helps transmit ECM signals to the cells, thereby influencing cell morphology, mobility, and cell cycle progression. It is closely related to the ECM-receptor interaction pathway, which bridges the extracellular environment with intracellular signaling mechanisms ( Honig and Shapiro, 2020 ). The genes related to collagen synthesis, such as COL4A1, COL4A2, COL12A1, COL1A2 , and COL6A1 , were significantly upregulated at 20 weeks of age (the early egg-laying stage). Additionally, the proteins COL12A1 and COL1A2 were also significantly upregulated in the proteome during the early egg-laying stage. However, we also observed a low correlation between gene and protein expression, as illustrated by the nine‑quadrant plot. This phenomenon highlights the complexity of gene expression regulation, since mRNA levels do not always directly correspond to protein abundance. Factors contributing to this discrepancy may include post‑transcriptional regulation, variation in translation efficiency, protein stability, or differences in subcellular localization of gene products. Such inconsistency is commonly observed in biological systems and underscores the importance of integrating transcriptomic and proteomic data to more comprehensively understand the molecular mechanisms of ovarian development. Similar discrepancies have been reported in other studies and may reflect the dynamic nature of protein synthesis and degradation processes, which do not necessarily align with gene expression levels ( Liu et al., 2016 , Naesens and Sarwal, 2010 ). COL12A1 encodes type XII collagen, a non-fibrillar collagen belonging to the collagen superfamily. Type XII collagen enhances tissue structural stability by binding to fibrillar collagen and is crucial for maintaining tissue structural integrity ( Izu and Birk, 2023 ). COL1A2 encodes the α2 chain of type I collagen, which, together with the α1 chain (encoded by COL1A1 ), constitutes the major component of type I collagen. This type of collagen is widely present in skin, bones, tendons, ligaments, and other connective tissues, providing essential mechanical strength and elasticity ( Gelse et al., 2003 ). Similarly, the integrin family genes ITGA8 ( Kulus et al., 2023 ), ITGA7 ( Li et al., 2023 ), and ITGB1 ( He et al., 2012 ), which play pivotal roles in cell adhesion and extracellular matrix interactions, were also significantly upregulated during this stage. These genes provide structural and ovarian tissue stability for the rapid proliferation of cells in early ovarian tissue. At the onset of egg production, poultry follicles and ovarian medulla require a large amount of collagen and stromal cells to support cell proliferation and development. Enhanced intercellular communication in the ovarian stroma of hens facilitates the rapid provision of collagen content for increased follicle attachment area. Genes related to collagen synthesis and cell adhesion play an important role in the rapid development of the ovarian stroma during the early egg-laying stage. During the peak egg-laying period (30 weeks), the genes involved in the steroid biosynthesis pathway ( MSMO1, FDFT1, DHCR24, LSS, CYP51A1, NSDHL, SQLE , and SC5D ) exhibited significant high expression. These genes are closely associated with ovarian development ( Pan et al., 2024 , Yu et al., 2022 ). Li et al. found that MSMO1, FDFT1, DHCR24, LSS, CYP51A1 , and SQLE are closely linked to steroid biosynthesis in the hypothalamic-pituitary-ovarian ( HPO ) axis of chickens, and are potential candidate genes influencing follicular development, luteinization, and reproductive system function ( Li et al., 2020 ). Enhanced expression of FDFT1 may promote steroid synthesis ( Zhang et al., 2024 ). Additionally, the significantly correlated gene-protein pairs CYP17A1 and HSD17B7 showed prominent differences during the peak period. This phenomenon reveals the central role of cholesterol and its derivative steroid hormones in the development of ovarian follicles and maintenance of the egg-laying cycle. Cholesterol, as an essential component of cell membrane structure, not only ensures the physical stability of cells and normal signal transduction but is also the fundamental substance for synthesizing steroid hormones such as estrogen ( Ershov et al., 2021 ; Ellinger and Chatuphonprasert, 2022 ). Estrogen particularly influences follicle maturation, promotes corpus luteum formation, and eggshell calcification, being a key factor in maintaining normal egg production in hens ( Mishra et al., 2019 ). During the peak egg-laying period, the demand for estrogen and other steroid hormones significantly increases, and the high expression of these genes may be a physiological adaptation to meet the high demand for steroid and hormones during the peak egg-laying period ( Çiftci, 2012 ). By enhancing steroid biosynthesis, these genes ensure the continuous production of estrogen and other related hormones, thereby supporting the ongoing development and maturation of ovarian follicles and maintaining high levels of egg production. In this study, we found that the ANXA family genes ANXA2 and ANXA5 , which are associated with cell proliferation and tissue development, were differentially highly expressed in both the transcriptome and proteome during the early egg-laying stage. Annexins are a multifunctional family of Ca 2+ -regulated membrane phospholipid-binding proteins that are highly conserved across different organisms. The N-terminal domain of Annexin A2 ( ANXA2 ) contains binding sites for P11 protein ( S100A10 ), plasminogen, and tissue plasminogen activator (t-PA) ( Kwon et al., 2005 ; Flood and Hajjar, 2011a ). The activity of plasmin is crucial in many physiological processes, such as ECM protein hydrolysis, cell migration, tissue repair/remodeling, and angiogenesis ( Pepper, 2001 ; Syrovets and Simmet, 2004 ; Flood and Hajjar, 2011b). Researchers have noted that the acetylation of ANXA2 enhances its binding to the epidermal growth factor receptor ( EGFR ), thereby activating the EGFR signaling pathway, which significantly affects the proliferation and apoptosis of ovarian granulosa cells( Zhou et al., 2024 ). ANXA2 plays an important role during the process of embryo implantation, promoting embryonic attachment and endometrial receptivity ( Wang and Shao, 2020 ). In a study on ovarian tissues of hens with different egg production rates, ANXA2 was identified as a candidate functional gene associated with high egg-laying performance ( Yang et al., 2008 ). Similarly, ANXA5 is a Ca 2+ -dependent phospholipid-binding protein that belongs to the annexin family, and its gene is induced by Gonadotropin-releasing hormone ( GnRH ) ( Crompton et al., 1988 ; Kawaminami et al., 2002 ). Studies have shown that GnRH agonists enhance the expression of ANXA5 and LHβ in the pituitary of hypogonadal female mice, while a decrease in ANXA5 mRNA levels inhibits LH secretion in pituitary cells ( Yonezawa et al., 2015 ). Additionally, in mouse ovaries, the expression of ANXA5 increases under hCG stimulation, particularly during the transition of follicular granulosa cells into luteal cells ( Tungmahasuk et al., 2018 ), indicating that ANXA5 plays an important role in regulating granulosa cell function under hormonal influence. During the early development of chicken ovaries, ANXA5 and ANXA2 may maintain ovarian tissue homeostasis and maturation by promoting cell migration, angiogenesis, and hormonal regulation. The protein encoded by the Oxysterol Binding Protein 2 ( OSBP2 ) gene belongs to the oxysterol-binding protein family and is mainly expressed in the retina, testes, and fetal liver ( Moreira et al., 2001 ). The OSBP2 protein plays a critical role in regulating intracellular lipid metabolism, particularly in the transport and distribution of cholesterol ( Charman et al., 2014 ). In reproduction, OSBP2 is especially important for the differentiation of germ cells, and its absence can lead to male infertility, characterized by abnormal sperm morphology and decreased motility ( Udagawa et al., 2014 ). These findings suggest that OSBP2 plays a key role in maintaining reproductive health and function. Although this gene has not yet been reported in the ovary, based on its high expression in early ovaries and its role in reproductive function, we hypothesize that the OSBP2 protein plays an important role in the maturation and development of the early ovaries of egg-laying hens by regulating the transport and metabolism of key lipids such as cholesterol. Legumain ( LGMN ) is a cysteine protease first isolated from pig kidney and is expressed in various tissues( Chen et al., 1997 ). LGMN can activate other enzymes or proteins, including CTSL, matrix metalloproteinase-2, progelatinase A, and Toll-like receptor 9 ( Sepulveda et al., 2009 ); it is involved in the remodeling of intracellular and extracellular matrices and mediates apoptotic pathways ( Zeeuwen et al., 2009 ). After the ovary begins egg production, the ovary continuously enlarges, undergoing extensive tissue remodeling. The protease LGMN may participate in the tissue remodeling process and maintain the proliferation of cells within the ovarian tissue. Endothelin receptor type A ( EDNRA ) plays a crucial role in ovarian maturation. EDNRA is a G protein-coupled receptor whose activation can induce the contraction of smooth muscle cells in the follicle wall, a key mechanism for follicle rupture and oocyte release ( Meidan and Levy, 2002 ). In the ovary, EDNRA is co-expressed with endothelin receptor type B ( EDNRB ) and promotes folliculogenesis, steroidogenesis, oocyte maturation, ovulation, and luteal function ( Cho et al., 2012 ). EDNRA and its ligand endothelin-1 promote oocyte maturation, mainly expressed in granulosa cells and cumulus cells surrounding the oocyte, and are critical for intercellular communication and oocyte maturation( Kawamura et al., 2009 ). In laying hens, the increased expression of EDNRA in the ovary may indicate that the follicles are preparing for maturation and ovulation. Cysteine-rich secretory protein LCCL domain-containing 2 ( CRISPLD2 ) is a member of the cysteine-rich secretory proteins, antigen 5, and pathogenesis-related 1 proteins ( CAP ) superfamily. It is associated with various functions such as cell differentiation, migration, inflammation, and immunity ( Gibbs et al., 2008 ). This gene has been widely reported in organ morphogenesis ( Quinlan et al., 2007 ). Studies have shown that this gene regulates the proliferation, apoptosis, and migration of fetal lung fibroblasts, aids in mesenchymal-epithelial signaling, enhances wound repair, and is associated with the abnormal expression of multiple ECM (extracellular matrix) genes that regulate lung development and repair ( Yoo et al., 2014 ). Interestingly, CRISPLD2 is regulated by P4 and progesterone receptor ( PGR ) in the uterus, with high expression during decidualization and sustained expression throughout pregnancy, and its expression is dysregulated in patients with endometriosis. Similar P4/PGR-dependent regulation of CRISPLD2 gene expression has been observed in rat granulosa cells ( Sriraman et al., 2010 ). CRISPLD2 may be related to hormonal response, tissue repair, and maintenance of the extracellular matrix in early ovarian development and egg production in chickens, ensuring normal ovarian function and continuity of egg production.

The authors declare that they have no competing financial interests. Supplementary figure 1. Quality control of transcriptomic data and analysis of differentially expressed genes. A. Pearson correlation matrix between samples; B. Gene expression distribution, with each box plot representing the distribution of all gene expression values within a sample, including median, quartiles, and outliers; C. Bar chart of DEGs, with blue bars representing the number of downregulated genes and red bars representing the number of upregulated genes.

Eggs are a crucial component of the human diet, and the ovaries are the primary organs for egg production in poultry. Proper ovarian development is essential for the reproductive performance of poultry. Unlike mammals, chickens have only one functional ovary on the left side, which contains numerous primordial follicles ( Johnson, 2015 ). As hens reach sexual maturity, their ovaries gradually develop and enter the laying phase. At this stage, follicles at various developmental stages are attached to the ovarian surface and can be classified into small white follicles ( SWF ), large white follicles ( LWF ), small yellow follicles ( SYF ), and large yellow follicles ( LYF ) based on their size and function. These follicles subsequently undergo selection and further development into preovulatory follicles ( PRF ), which are ranked by size as F1, F2, F3, F4, F5, and F6. After ovulation, 2–4 post-ovulatory follicles remain, which lack oocytes and are interconnected within the ovarian stroma ( Onagbesan et al., 2009 ; Hlokoe et al., 2022 ). In addition, the ovary regulates follicular functions through autocrine and paracrine signaling pathways, involving a series of gene transcription and protein expression events under complex regulatory mechanisms ( Zhao et al., 2023 ). Egg production is also influenced by multiple factors, including nutritional status ( Liu et al., 2020 ) and environmental conditions ( Konkol et al., 2020 ). Although there is increasing research on the mechanisms of ovarian development, studies focusing on the initial egg-laying stage to the peak period in hens are still limited. In recent years, with the rapid development of high-throughput sequencing and mass spectrometry technologies, these techniques have been widely used. Transcriptomics can provide regulatory information at the gene level( Ghosh et al., 2018 ), while proteomics reveals regulatory mechanisms at the protein level( Jaswal et al., 2021 ). Data from a single omics approach can only offer information on a specific molecular layer within a biological sample and cannot comprehensively reveal the complex biological processes and interactions of cells or tissues. Currently, multi-omics integrative analysis has become an important and routine method for analyzing the molecular mechanisms underlying complex traits in animals. Researchers predict the potential genes regulating chicken feather color by combining transcriptomics and proteomics( Wang et al., 2019 ). Li et al. revealed the intrinsic mechanisms underlying the earlobe color differences in Jiangshan Black-bone chickens through transcriptomic and proteomic analysis( Li et al., 2024 ). A recent study revealed the critical role of lactate dehydrogenase A ( LDHA ) in the ovulation process of domestic chickens through combined proteomic and transcriptomic analysis( Nie et al., 2024 ). The Taihe black-boned silky fowl is a high-quality indigenous chicken breed from the Taihe County, Jiangxi Province, China. This breed is renowned for its ten unique characteristics and its dual value in both medicine and food( Tu et al., 2009 ). Its use as a dietary therapy product dates back over 2200 years( Mi et al., 2018 ). However, the TBSF currently exhibits a delayed onset of egg production, and there is a lack of systematic understanding of ovarian function during the egg-laying process. These factors limit its economic benefits and developmental potential. Comprehensive analysis of the transcriptome and protein abundance in ovarian tissues is performed. By comparing the transcriptome and quantitative proteome differences between the ovaries of TBSF at 30 weeks and 20 weeks of age, the study identifies key differentially expressed genes ( DEGs ) and differentially abundant proteins ( DAPs ) involved in the development from the onset of egg production to the peak production period. These findings contribute to a deeper understanding of ovarian development in this local chicken breed, offering valuable insights for breeding strategies.

The current work was funded by the 10.13039/501100013076 Major Scientific and Technological cooperation between Zhejiang University and Taihe County Government, grant number 2021-KYY-517102-0023.