Identification of Key Genes Regulating Body Weight in Qingyuan Partridge Chickens During Development Using RNA-Sequence Analysis

TL;DR

This study identifies key genes involved in body weight regulation in Qingyuan partridge chickens during development using RNA sequencing.

Contribution

The study systematically identifies key genes and pathways associated with body weight regulation in Qingyuan partridge chickens.

Findings

Transcriptomic analysis identified 3521 genes specifically expressed at one day of age, enriched in ribosomal biosynthesis and cell proliferation pathways.

WGCNA identified two gene modules containing 1563 hub genes significantly correlated with body weight.

26 key genes, including CALM2, HSP90AA1, and CHRND, were closely associated with muscle growth and body weight regulation.

Abstract

Background: The Qingyuan partridge chicken is a high-quality local chicken breed in China. Its weight gain directly affects breeding efficiency. This study used RNA sequencing to analyze gene expression dynamics in the breast muscle tissue of Qingyuan partridge chickens at 1, 35, 70, and 105 days of age. Methods: This study employed RNA-sequencing, integrated with differential expression analysis, weighted gene co-expression network analysis (WGCNA), and short time-series expression miner (STEM) analysis, to systematically investigate the transcriptomic dynamics in breast muscle tissue across four developmental stages. Results: Phenotypic analysis revealed a significant increase in both body weight (BW) and breast muscle weight with age (p < 0.05). Transcriptomic analysis identified 3521 genes specifically expressed at the age of one day compared with the other 3 ages. These were significantly enriched in pathways related to ribosomal biosynthesis, cytoskeletal regulation, and cell proliferation (p < 0.05). Turquoise and black modules were identified by WGCNA, containing 1563 hub genes, which significantly correlated with BW. Integration of differentially expressed genes and STEM analysis selected 26 BW-related key genes closely associated with muscle growth, including calmodulin 2 (CALM2), heat shock protein 90 alpha family class A member 1 (HSP90AA1), and cholinergic receptor nicotinic delta subunit (CHRND). Protein–protein interaction analysis revealed two functional networks centered around these genes. Enrichment analysis of the STEM profiles indicated that upregulated genes were significantly enriched in autophagy and the ErbB, FoxO, mTOR, and insulin signaling pathways, while genes related to the ribosome, cell cycle, and PPAR signaling pathways were downregulated. Conclusions: This study identified BW-related key genes and pathways, enriching our knowledge of the functional maintenance of chicken BW.

Figures

Tables

Table 3. Sequencing Data Quality Assessment in Qingyuan Partridge Chicken.

Sample	Raw Reads	Clean Reads	Aligned Read Count	Q20 %	Q30 %	GC %	Quality Control Rate %	Mapping Rate %
D1A	48,131,334	46,608,636	40,773,916	99.19	96.85	50.04	96.84	87.48
D1B	47,704,262	46,292,478	40,261,096	99.06	96.53	51.09	97.04	86.97
D1C	51,243,234	49,686,352	43,321,517	99.10	96.67	50.75	96.96	87.19
D1D	43,895,822	42,458,668	37,138,147	99.14	96.72	50.94	96.73	87.47
D35A	50,289,340	48,672,028	40,186,747	99.16	96.92	53.10	96.78	82.57
D35B	50,494,292	48,555,720	39,865,927	99.13	96.63	53.18	96.16	82.10
D35C	46,279,926	44,640,068	37,439,475	99.15	96.78	52.54	96.46	83.87
D35D	48,220,244	46,459,386	38,646,409	99.13	96.66	53.24	96.35	83.18
D70A	46,393,034	44,949,420	37,243,300	99.04	96.49	53.07	96.89	82.86
D70B	49,560,150	47,562,514	40,093,168	99.18	96.84	52.25	95.97	84.30
D70C	50,313,524	48,594,060	41,901,512	99.17	96.8	51.94	96.58	86.23
D70D	47,830,538	46,435,794	41,016,043	99.07	96.61	51.34	97.08	88.33
D105A	52,252,684	50,752,254	43,153,095	99.15	96.73	52.18	97.13	85.03
D105B	45,424,256	44,000,958	36,056,874	99.16	96.7	51.94	96.87	81.95
D105C	55,872,666	54,028,506	45,718,646	99.17	96.63	52.08	96.70	84.62
D105D	48,963,352	47,198,036	39,645,539	99.19	96.9	51.78	96.39	84.00

Full text

1. Introduction

The Qingyuan partridge chicken is considered a quintessential representative of China’s high-quality local chicken breeds, with a distinguished reputation for its succulent meat and unique flavor [1]. Its meat possesses significant economic value and market potential due to its high quality and unique characteristics. However, its primary drawback is its relatively low body weight (BW), which directly impacts meat production rates and economic returns. Body weight is widely regarded as a pivotal metric for evaluating individual growth, development, and economic performance in poultry production [2]. This parameter exerts a direct influence on the profitability of poultry farming. Weight is an extremely complex quantitative trait, governed dynamically by the synergistic interaction of major genes and numerous minor ones [3].

Genomic research has extensively explored the genetic basis of economic traits in domestic chickens. This work has led to the identification of numerous genes implicated in growth and carcass characteristics. Key among these are insulin-like growth factors (IGFs) [4,5], growth hormone secretagogue receptor (GHSR) [6], and adhesion G protein-coupled receptor G6 (ADGRG6) [7]. Additionally, studies across multiple species like chicken and duck have established a positive correlation between insulin-like growth factor 2 mRNA binding protein 1 (IGF2BP1) and traits such as breast muscle weight (BMW) and overall body size [8,9,10]. In native Korean and indigenous Iranian chickens, body weight was correlated with variants in WDR37, KCNIP4, SLIT2, PPARGC1A, MYOCD, ADGRA3, DCN, MEOX2, and CACNB1 [11,12]. Studies on Asian Game and Bantam breeds have proposed MTM1 and SFPR2 as potential candidate genes for body size [13]. Furthermore, associations with growth traits have been reported for SOX5, IGF1, and several MYH1 genes (MYH1A, MYH1B, MYH1D) in commercial broilers and Jinghong layer chickens [14,15].

Body growth involves a series of biological events, including the balance between cell proliferation and differentiation, the dynamic balance between protein synthesis and degradation, reprogramming of energy metabolism, and multi-level regulation of endocrine and paracrine signaling [16]. Research indicates that growth traits in vertebrates are tightly regulated by multiple classical signaling pathways, including the insulin/IGF axis, mTOR, MAPK, and TGF-β pathways [17,18,19,20]. Key genes within these pathways and their interaction networks collectively determine the final growth phenotype. Growth hormone (GH) and insulin-like growth factor 1 (IGF1) can act synergistically or independently through the GH-GH receptor (GHR)-IGF1 signaling pathway to promote skeletal muscle growth and enhance muscle mass [21].

Research on the genetic mechanisms underlying poultry BW has focused on imported commercial breeds, while research on gene expression during body development and the key regulatory genes in local breeds like the Qingyuan partridge chicken is lacking. RNA sequencing (RNA-Seq) technology can be used to systematically decipher this aspect [22]. It can map dynamic gene expression changes across developmental stages and identify differentially expressed genes (DEGs) associated with traits like BW. Combining this technique with bioinformatics methods like weighted gene co-expression network analysis (WGCNA) and short time-series expression miner (STEM) can effectively screen key genes with core regulatory functions from large datasets.

This study was based on the hypothesis that BW gain in Qingyuan partridge chickens is regulated by the temporal dynamics of specific gene expression and key signaling pathways during development. Specifically, we hypothesized that in breast muscle tissue across distinct developmental stages (1, 35, 70, and 105 days of age), there exist core sets of genes and functional modules significantly correlated with BW traits. The expression patterns of these genes and pathways—particularly those related to muscle growth, metabolism, and cell proliferation—were expected to change systematically with age, collectively governing BW increase. By integrating transcriptomic profiling, WGCNA and STEM, we aimed to systematically identify these key BW-regulating genes and pathways, thereby elucidating the molecular mechanisms underlying BW development in Qingyuan partridge chickens.

2. Materials and Methods

2.1. Animal Genetic Background, Phenotypes, and Sample Collection

Qingyuan partridge roosters were obtained from Guangxi Jinling Agriculture and Animal Husbandry Group Co., Ltd., Nanning, China, to ensure a uniform genetic background. All chickens were reared under standard, uniform conditions with ad libitum access to water and the same diet. A total of 16 left breast muscle samples were collected from four critical developmental stages: day 1 (D1, neonatal stage), day 35 (D35, early rapid growth phase), day 70 (D70, pre-pubertal stage), and day 105 (D105, mature stage), with four biological replicates per time point. These stages were selected to comprehensively cover the key physiological transitions in muscle development. Body weight (BW) was recorded at all four time points. Breast muscle weight (BMW) was recorded at D35, D70, and D105. BMW was not measured at D1 because the breast muscle tissue was too small for accurate dissection and weighing. All collected tissue samples were immediately snap-frozen in liquid nitrogen and stored at −80 °C until further analysis.

2.2. RNA Library Construction and RNA-Seq Analysis

Total RNA was extracted from the breast muscle tissue samples using Trizol reagent (Invitrogen, Carlsbad, CA, USA). The quality of the RNA was assessed after separation by electrophoresis on a 1% agarose gel, and the RNA concentration was determined by spectrophotometry using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA, USA). RNA-Seq was performed by Beijing Novogene Co., Ltd. Total RNA was enriched for mRNA using oligo(dT) magnetic beads and then fragmented under high-temperature conditions. First-strand cDNA was synthesized using random hexamers as primers and M-MuLV reverse transcriptase, followed by second-strand synthesis with DNA polymerase I and RNase H [11]. After end repair, A-tailing, and adapter ligation, cDNA fragments of 370–420 bp were selected for PCR amplification for library construction. After quality control, 150 bp paired-end sequencing was performed on an Illumina Novaseq platform (Illumina, San Diego, CA, USA). Raw data underwent fastp quality control to remove adapters and low-quality reads, where bases with Qphred values ≤5 constitute over 50% of the total read length, yielding high-quality cleaned data. Q20, Q30, and GC content were calculated. Clean reads were aligned to Gallus gallus genome assembly GRCg6a using HISAT2, and the mapped reads were then assembled and their abundance estimated as FPKM (Fragments Per Kilobase of transcript per Million mapped reads) using String Tie [12]. These data indicate satisfactory sample quality, highly accurate sequencing, and strong usability for subsequent analyses (Table 3).

2.3. Screening for DEGs

To obtain count data for differential expression analysis, gene-level raw read counts were quantified from the alignment files using the featureCounts tool (v2.0.3) from the Subread package. Gene alignment statistics were obtained using Feature Counts from Subread, and gene expression levels were calculated using the FPKM. In DESeq2, the raw counts are normalized using the median-of-ratios method, and a negative binomial generalized linear model is fitted. The Wald test was used to assess statistical significance. DEGs between two groups (comparison couples: D35 vs. D1, D70 vs. D1, D126 vs. D1, D70 vs. D35, D126 vs. D70) were identified using DESeq2 with the thresholds set to p < 0.05 and |Fold Change| > 1.5 [23]. Bidirectional clustering analysis of DEGs was conducted based on FPKM values.

2.4. Construction of the WGCNA and Screening of Hub Genes

Weighted gene co-expression network analysis was performed using the WGCNA package in R 4.5.0. Abnormal samples in WGCNA will affect the analysis of module results, so before WGCNA, abnormal outlier samples were removed, 15 samples were clustered, and single samples were found. The 13,867 genes identified were used to build a co-expression module, by first calculating the correlation coefficient of genes. To determine whether genes have similar expression patterns, it is necessary to set an appropriate soft threshold [24], so that the connection between genes in the network obeys the scale-free network. When the scale-free fitting index is 0.85, the soft threshold is 3, indicating good network connectivity. Therefore, 3 was chosen as the most appropriate soft threshold for building a co-expression module. Then, a co-expression network was built through an appropriate soft threshold, and genes were classified according to their expression pattern [25]. Genes with similar expression pattern were classified into one module, and different colors represent different modules. To associate modules with traits, principal component analysis was performed on all genes within the modules [26]. Modules significantly correlated with BW (p < 0.05) were selected for further analysis, and hub genes were screened using the criteria |GS| > 0.6 and |MM| > 0.8.

2.5. Short Time-Series Expression Miner Analysis

The non-parametric short time-series expression miner (STEM) clustering algorithm, version 1.3.13 [27], was used to cluster and visualize possible profiles and expression changes in DEGs screened from three comparison groups (D1 vs. D35, D35 vs. D70, D70 vs. D105). STEM was run using the log-normalized data option. The chosen parameters (a maximum unit change of 2 between time points and a maximum of 20 model profiles) represent the software’s default settings, aimed at identifying major biological expression trends. Crucially, the key gene clusters derived from this analysis were independently validated by the results of co-expression network analysis (WGCNA), ensuring the robustness of the clustering outcomes. The maximum unit change in model profiles between time points was adjusted to 2, and the maximum number of model profiles was set to 20.

2.6. Functional Module Enrichment Analysis

Enrichment analysis was performed on genes specifically expressed at D1 and in key gene profiles (Profile 4 and 17) filtered by STEM. Gene Ontology (GO; http://www.geneontology.org/) was used to identify relevant GO terms divided into three categories: biological process (BP), cellular component (CC), and molecular function (MF). Pathway enrichment analysis using the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/) identified metabolic and signal transduction pathways associated with the regulated genes. Key modules of the hub genes were investigated using the KOBAS (version 3.0) software (http://kobas.cbi.pku.edu.cn/kobas3/?t=1), which identified enriched KEGG pathways and GO terms associated with significant phenotypic traits (p  <  0.05).

2.7. Protein–Protein Interaction Analysis

Protein–protein interactions (PPIs) associated with DEGs involved in cellular differentiation regulatory networks were explored using the STRING (version 12.0) gene/protein interaction database [28,29]. The PPI network was constructed by setting a minimum threshold for interaction reliability (interaction score > 0.15) [28]. The PPI network diagram was constructed using STRING (https://string-db.org/) and visualized with Cytoscape 3.7.1 software.

2.8. Quantitative Real-Time PCR (qPCR)

Six key candidate genes (CALM2, CHRND, HSP90AA1, SLC8A3, CD9, PENK) were strategically selected for qPCR validation based on the following criteria: (1) they were identified as top hub genes in the WGCNA co-expression modules most significantly correlated with muscle growth traits; (2) they exhibited highly significant differential expression across developmental stages in RNA-seq analysis; and (3) they represent core biological pathways implicated in muscle development and metabolism. Primers were designed by NCBI/Primer-Blast software based on the gene sequences, and sequence alignment was performed by BLASTN (Table 1). Total RNA was the same RNA used for sequencing and was reverse transcribed from 2.0 µg total RNA per sample using a FastQuant RT Kit (Tiangen, Beijing, China) to obtain the cDNA. The qPCR reaction was performed in a final volume of 10 µL, comprising 5 µL of 2× iQTM SYBR Green Supermix, 0.5 µL each of upstream and downstream primers, 3 µL of cDNA and 1 μL ddH_2_O. The reaction condition involved an initial denaturation at 95 °C for 3 min, followed by 40 amplification cycles of 95 °C for 3 s and 60 °C for 34 s, on a QuantStudio 7 Flex system (Applied Biosystems, Foster City, CA, USA). All sample amplifications were performed in triplicate. β-actin served as an internal reference gene, and the relative gene expression was calculated using the 2^−ΔΔCt^ method.

2.9. Statistical Analysis

Statistical analysis of differences between means was performed using IBM SPSS Statistics for Windows, Version 26.0 (IBM Corporation, Armonk, NY, USA). A p-value of less than 0.05 was considered statistically significant. Data are presented as means  ±  SEM.

3. Results

3.1. Phenotypic Analysis

We measured the BW and BMW of 16 Qingyuan partridge chickens at various developmental stages (Table 2). BW increased significantly over time (p < 0.05). The BW was 24.83 ± 0.44 g at 1 day old (D1), 430.75 ± 26.18 g at 35 days old (D35), 1190.00 ± 34.45 g at 70 days old (D70), and 2292.5 ± 60.41 g at 105 days old (D105). The BMW was 29.95 ± 3.48 g at D35, 106.75 ± 5.94 g at D70, and 262.9 ± 19.24 g at D105. Pairwise comparisons showed significant differences between all age groups (p < 0.05).

3.2. Transcriptome Sequencing Quality Assessment and Alignment

We obtained 43 million raw reads (Table 3). After quality control, 42 million valid sequences were retained for further analysis, with a quality control rate of 95.97–97.13%. The Q20 rate exceeded 99.04% for all samples, Q30 coverage exceeded 96.49%, and GC content remained above 50.04%. When aligned to the junglefowl reference genome, the effective sequence alignment rate exceeded 81.95%.

3.3. Identification and Functional Enrichment Analysis of Genes Specifically Expressed at 1 Day Ago

The number of DEGs identified in comparisons involving D1 breast tissue was high, with 7238, 7349, and 8533 DEGs detected for the contrasts of D35 vs. D1, D70 vs. D1, and D105 vs. D1, respectively. The lowest number of DEGs (1284) was observed in D70 vs. D35 (Figure 1A). After excluding duplicates from pairwise comparisons in D35 vs. D1, D70 vs. D1, and D105 vs. D1, we retained 3521 DEGs unique to D1 (Figure 1B).

KEGG and GO enrichment analyses were performed using 2512 known genes among the 3521 retained. KEGG pathway analysis revealed significant enrichment in key pathways governing cell proliferation, differentiation, and basal metabolism (Figure 1C, Table S1). Most prominent among these was the ribosome pathway, enriched for 62 ribosomal protein genes, including RPS13, RPS11, and RPL37A, indicative of highly active protein synthesis. DEGs including PDGFRA, PDGFRB, ITGA7, VCL, and ACTN1 were enriched in pathways regulating the actin cytoskeleton and focal adhesion. This suggests robust cytoskeletal remodeling and enhanced cell–matrix adhesion—processes essential for myoblast migration, fusion, and myotube formation. Key signaling pathways regulating muscle development, including apelin (e.g., ACTA2, SLC9A1, PIK3R4, SLC8A1), adrenergic signaling in cardiomyocytes (e.g., SLC9A1, PIK3R5, SLC8A1, MYH7), MAPK (e.g., TGFB2, PDGFRA, MAP3K4, GADD45B, MAPK9, PDGFRB, FGF16, PDGFB, TGFBR1), calcium signaling (e.g., PDGFRB, PDGFRA, SLC8A1), and cell cycle (e.g., TGFB2, GADD45B, PCNA), were also significantly enriched.

GO analysis corroborated the KEGG findings, showing significant enrichment in processes and components essential for rapid growth (Figure 1D, Table S2). Translation (RPL family) and DNA replication (e.g., PCNA, SLC19A1) were the terms in BP, underscoring hyperactive proliferation and anabolism. Regulation of cell shape (e.g., CSF1R, PLXND1, MYH9) and sarcomere organization (e.g., ACTN1, TNNT3, CAPN3) were also enriched, highlighting ongoing differentiation. For CC, genes were predominantly localized to the ribosome (RPL family), actin cytoskeleton (e.g., MYH9, MYO1F, MYO1G), focal adhesion (e.g., NRP1, ACTN1, MYH9), and adherens junction (e.g., CRB1, MYH9, FAT3), emphasizing the machinery for biosynthesis and structural morphogenesis. Accordingly, “structural constituent of ribosome” (RPL and RPS families) was the MF term.

3.4. Identification of Key Genes Related to Body Weight

We identified 60 common DEGs after performing pairwise comparisons across the four age groups (Figure 2A). Among these, 52 DEGs were selected based on consistent growth trends, 20 upregulated and 32 downregulated. Similar expression patterns were observed in breast muscle DEGs between D70 and D105, whereas completely opposite expression patterns were observed in DEGs between D1 and D105 (Figure 2B). Following standard preprocessing, we retained 15 samples and 13,867 genes for downstream analysis. Turquoise and black modules were identified as significantly correlated with BW, with 1563 hub genes identified across both modules (p < 0.05; Figure 2C–E).

By integrating the 52 DEGs, we found 39 genes common to both modules (Figure 3B). We performed STEM analysis on the 5125 genes obtained by taking the intersection of DEGs across the three comparison groups (D1 vs. D35, D35 vs. D70, D70 vs. D105). STEM analysis of expression profiles across adjacent age groups clustered 5125 DEGs into 20 distinct groups based on expression profile similarity, with significant enrichment in seven of them (Figure 3A). Seven profiles represent different trends. Due to the similarity between the expression patterns of Profile 17 and 4 and the growth trend, which suggests they may play a key role in growth regulation, we have therefore prioritized them for further analysis. Expression profile analysis revealed that 578 genes in Profile 17 were downregulated, while 532 genes in Profile 4 were downregulated. Intersecting with the 39 genes yielded 26 DEGs consistent with both the expression profiles and BW (Figure 3B).

PPI network analysis showed that these 26 key genes formed two protein interaction networks: 14 genes, including CALM2 and HSP90AA1, formed one network, while the others, including MATN2 and PENK, formed the second (Figure 3C). Among the 26 genes, SLC8A3, CHRND, CALM2, and ISCU are directly associated with muscle growth; GNMT, AHCY, and NPR2 are linked to body weight and metabolic regulation; HIF1AN, MATN2, CD9, and HSP90AA1 indirectly participate in regulating muscle growth and metabolism; and others, including FOXE3, JMJD6, LAPTM4B, and PENK, are primarily associated with non-muscle biological processes.

The relative mRNA expression levels of six DEGs (CHRND, CALM2, HSP90AA1, SLC8A3, CD9, and PENK) aligned with the observed BW trends (Figure 3D), and showed significant correlations with BW (p < 0.05). Specifically, CALM2 and HSP90AA1 correlated positively with BW (R^2^ > 0.9), whereas CHRND, SLC8A3, CD9, and PENK correlated negatively with BW (0.4 > R^2^ > 0.6).

3.5. Enrichment Analysis of Genes Consistent with the Key Gene Expression Profiles

Enrichment analysis of KEGG and GO pathways was performed for genes from Profiles 17 and 4 (Figure 4A–D, Tables S3–S6). Genes in Profile 17 were significantly enriched in multiple signaling pathways associated with muscle growth and development. Key pathways, such as autophagy–animal, ErbB, FoxO, mTOR, and insulin signaling (e.g., MAPK3, PIK3R3), are directly involved in muscle protein synthesis, cell proliferation, and metabolic regulation. Key genes, including CALM2, MAPK3, PIK3R3, and RAF1, were recurrently identified across multiple pathways, suggesting roles as central regulators coordinating muscle growth and energy homeostasis.

Genes in Profile 4 were prominently enriched in pathways related to the ribosome (RPL and RPS families), RNA transport (EIF3 series), ribosome biogenesis in eukaryotes (e.g., DKC1 and UTP4), and the cell cycle (e.g., CDK1 and E2F1), underscoring a central role in protein synthesis, a fundamental process in muscle growth. Additional enrichment in the cellular senescence and PPAR signaling pathways implies broader involvement in cell proliferation, differentiation, and lipid metabolism. Key cell cycle regulators, including CDK1, E2F1, CCNB2, and CCNA2, further support the importance of this profile in coordinating cellular proliferation and muscle development.

GO enrichment analysis of Profile 17 highlighted significant involvement in biological processes related to protein metabolism, energy regulation, and cell growth. Major BP terms included ubiquitin-dependent protein catabolic process (e.g., ATE1, RNF11, RNF14), protein refolding (e.g., DNAJA1, DNAJB2, HSPA8, HSPA2), and chaperone-assisted refolding (e.g., DNAJB1, HSPA8, HSPA5, HSPA2), underscoring the module’s role in maintaining proteostasis. Genes were also implicated in the insulin-like growth factor receptor signaling pathway, cellular response to amino acid starvation, and negative regulation of TOR signaling, reflecting functions in nutrient sensing, energy metabolism, and muscle growth regulation. Key genes such as HSPA8, HSPA2, HSPA5, RAF1, MAPK3, and PRKAA2 appeared across multiple processes, indicating coordinated regulation of protein turnover, energy balance, and stress response, impacting muscle development and metabolic balance.

GO analysis of Profile 4 revealed strong enrichment in translation, ribosome biogenesis and function, and cell cycle regulation. BP terms such as translation (RPL and RPS families), cytoplasmic translation (RPL family), translational initiation (EIF3 series), and G1/S transition of the mitotic cell cycle (e.g., CCNA2, CCNB1, CCND1, EIF4EBP1, E2F1, RCC1) emphasized the profile’s core role in protein synthesis and proliferative control. The profile was also associated with heart development, angiogenesis, the Wnt signaling pathway, and negative regulation of apoptosis, suggesting broader roles in tissue development and cellular homeostasis. Key genes—including ribosomal proteins (RPL and RPS families), translation initiation factors (EIF3 series), and cell cycle regulators (CCNA2, CCNB1, CCND1, CDK1)—co-occurred across multiple pathways, collectively regulating protein synthesis and cell division, which are indispensable for muscle growth and progression.

4. Discussion

This study investigated the transcriptomic dynamics underlying breast muscle growth in Qingyuan partridge chickens, a high-quality indigenous breed, across four postnatal stages (days 1, 35, 70, and 105). The central hypothesis was that key genes and pathways regulating growth could be systematically identified by integrating multidimensional bioinformatics analyses. This approach, which combined differential expression analysis, WGCNA, STEM time-series clustering, and functional enrichment, successfully revealed critical molecular networks associated with body and breast muscle weight development. Our findings not only validate this integrative strategy but also provide novel insights into the genetic basis of growth in indigenous poultry.

Our phenotypic data clearly delineated the growth curve of Qingyuan partridge chickens. Both BW and BMW increased significantly with age. BW maintained a high rate of weight gain between 1 and 105 days. BMW also exhibited a sustained and significant increase with age, indicating that breast muscle development progressed in tandem with overall growth.

Analysis of genes specifically expressed in 1-day-old chicks revealed a highly active molecular landscape. The substantial set of 3521 genes specifically highly expressed at D1 was significantly enriched in pathways related to ribosome biogenesis, cytoskeleton regulation, focal adhesion, and the cell cycle. This result is consistent with biological expectations. The immediate post-hatch period is critical for the rapid development and maturation of muscle tissue, during which myoblasts undergo vigorous proliferation, migration, and fusion to form myotubes that ultimately mature into muscle fibers [30,31]. The abundant enrichment of ribosomal pathways indicates an intense state of protein synthesis in the breast muscle of 1-day-old chicks, providing the material basis for rapid muscle formation. Concurrently, the reprogramming of cytoskeleton and focal adhesion pathways is essential for morphological changes and cell–cell interactions among muscle fibers [32,33]. These findings collectively indicate that at one day of age, Qingyuan chickens are undergoing fundamental construction and intense morphogenesis, laying a solid structural foundation for subsequent rapid weight gain.

In this study 26 key genes related to BW regulation were identified through a multistep screening strategy. CALM2, HSP90AA1, CHRND, SLC8A3, ISCU, GNMT, and NPR2 form a complex regulatory network. Among them, CALM2 (calmodulin 2) is a core sensor in the calcium signaling pathway, which plays a central role in muscle formation, cell proliferation and differentiation, and energy metabolism [34]. Calmodulin binds Ca^2+^ with high affinity [35] and is considered the most significant Ca^2+^ signal transducer in cells [36]. If Ca^2+^ is not transported to muscle fibers, its concentration in the sarcoplasmic reticulum may increase [37,38], affecting numerous biological processes, including muscle contraction, oxidative stress, inflammation, and glycogen metabolism [35,39,40]. Disruption of calcium signaling pathways has been associated with various myopathies in humans [41] and chickens [37,42]. Furthermore, reduced CALM2 levels have been shown to result in inadequate intracellular calcium transport, which can lead to a defective muscle contraction [35]. HSP90AA1 (heat shock protein 90 alpha family class A member 1) is a crucial molecular chaperone that regulates cell cycle, apoptosis, and stress responses by stabilizing the conformation of various signaling proteins, and is vital for muscle development [43]. CHRND (cholinergic receptor nicotinic delta subunit) functions at the neuromuscular junction, potentially influencing muscle innervation and function. CHRND overexpression might lead to abnormal cell proliferation or cell communication that promotes excessive fibrous connective tissue deposition, the direct cause of wooden texture in muscles [44]. SLC8A3 (solute carrier family 8 member A3) is involved in maintaining intracellular calcium homeostasis. ISCU (iron-sulfur cluster assembly enzyme) is a key factor in mitochondrial energy metabolism, while GNMT (glycine N-methyltransferase) and NPR2 (natriuretic peptide receptor 2) are associated with amino acid and cardiovascular metabolism. The functional diversity of these proteins indicates that weight regulation is a complex process involving multiple biological layers: muscle structural development, energy metabolic homeostasis, intracellular signal transduction, and neuroregulation. Protein–protein interaction network analysis further categorized these key genes into two functional modules. One module is centered around CALM2 and HSP90AA1, suggesting the pivotal pivotal hub role of calcium signaling and protein homeostasis regulation in the growth of Qingyuan partridge chickens. Functional enrichment analysis of the two key STEM expression profiles, Profiles 17 and 4, to which these genes belong, supported the above findings at the pathway level. Genes in Profile 17 (upregulated trend) were significantly enriched in autophagy and the ErbB, FoxO, mTOR, and insulin signaling pathways. These are well-established core pathways regulating cell growth, proliferation, and metabolism. The mTOR pathway is a central controller of protein synthesis in response to nutrient and energy status [45,46]. The insulin and ErbB pathways activate downstream pro-growth cascades upon receiving external stimuli [47]. Under energy deficiency or stress, the FoxO pathway and autophagy negatively regulate mTOR and promote protein degradation and energy recycling to maintain homeostasis. The synergy of these pathways ensures the balance between the synthesis and degradation of muscle proteins. Conversely, Profile 4 (downregulated trend) was enriched in ribosome, RNA transport, and cell cycle pathways. This aligns with the post-D1 developmental trend, where rapid cell proliferation and protein synthesis activities gradually stabilize as muscle tissue matures.

Our results show both consistency and uniqueness compared with findings in other chicken breeds. For instance, we confirmed the conserved role of the IGF signaling pathway and the MYH (myosin heavy chain) gene family (significantly expressed among DEGs) in muscle growth. However, key genes identified in our study, including CHRND, CALM2, and HSP90AA1, have not been emphasized as core regulators in previous studies on BW in commercial broilers. This discrepancy may stem from the distinct genetic background of Qingyuan partridge chickens as an indigenous Chinese breed, which differs fundamentally from that of foreign commercial breeds like Arbor Acres Chicken. Therefore, our research fills a gap in the study of molecular mechanisms governing growth in Qingyuan partridge chickens. Furthermore, the discovered genes hold value as potential targets for molecular breeding of local chickens. Finally, qPCR validation of six genes showed significant correlation with BW, demonstrating the reliability and accuracy of the sequencing data in this study.

This study systematically revealed the weight-related transcriptomic dynamics during development in Qingyuan partridge chickens. We propose the following model: weight gain in Qingyuan partridge chickens is regulated by a multigene network. Pathways related to fundamental cell proliferation and structural formation are highly activated after hatching (D1). Subsequently, the regulatory focus shifts to fine-tuning protein metabolism and energy homeostasis, predominantly governed by signaling pathways such as mTOR, insulin, and FoxO. Among the activated genes, CALM2, HSP90AA1, and others occupy central positions in this regulatory network. These findings deepen the fundamental understanding of poultry growth biology and provide valuable genetic resources and novel candidate genes for marker-assisted and genomic selection breeding in Qingyuan partridge chickens. Potentially, these resources could help develop new lines that retain high-quality meat characteristics while exhibiting superior growth performance.

5. Conclusions

This study systematically investigated the molecular mechanisms governing body weight (BW) growth in Qingyuan partridge chickens by integrating phenotypic and transcriptomic analyses across developmental stages. Our findings highlight a distinct early developmental program, with 3521 genes specifically expressed in breast muscle at D1 enriched in ribosome biogenesis, cytoskeletal regulation, and cell proliferation—processes crucial for rapid tissue expansion and protein synthesis post-hatch.

Through integrated WGCNA and STEM analyses, we identified 26 core BW-associated genes, including CALM2, HSP90AA1, and CHRND. Functional analysis indicates that these genes collectively modulate the balance between muscle protein synthesis and degradation via key signaling pathways such as mTOR, insulin, and FoxO. These results provide novel molecular insights into the genetic regulation of growth in this indigenous breed and establish a functional link between transcriptional dynamics and weight-related traits.

Beyond confirming known growth-related pathways, this study proposes a coordinated regulatory network centered on protein turnover and cellular remodeling during postnatal development. The identified core genes represent promising candidates for molecular breeding aimed at optimizing growth performance in Qingyuan partridge chickens. Future studies should validate the roles of these candidate genes using functional assays and explore their potential interactions with environmental or nutritional factors to further refine selection strategies.