Preliminary comparative transcriptomic profiling of leg muscle in Pekin and Hanzhong Ma ducks

Introduction

Duck farming represents an important component of the poultry sector because of its contribution to meat and egg production, adaptability to diverse agro-climatic conditions and suitability for low-input farming systems. In many Asian countries, particularly China and India, ducks are valued for their ability to efficiently utilize natural feed resources and maintain productivity under semi-intensive and scavenging production systems. Improvement of economically important traits such as muscle growth, feed efficiency and adaptability has therefore become a major objective in modern duck breeding programs. Skeletal muscle is one of the most important tissues influencing growth performance and meat yield in poultry species. Apart from its essential role in locomotion and body support, skeletal muscle also contributes significantly to metabolic regulation and energy homeostasis within the body (Mitchell et al., 2002) [14] . Muscle development is regulated through myogenesis, a highly coordinated biological process involving the proliferation, differentiation and fusion of myogenic precursor cells. These processes are controlled by complex gene regulatory networks and signaling pathways that influence muscle fiber formation, hypertrophy and metabolic activity (Buckingham, 2017) [4] . Variations in the expression of genes associated with these pathways can lead to substantial differences in muscle growth and production efficiency among breeds. The concept of the transcriptome emerged following the development of gene expression profiling approaches such as Serial Analysis of Gene Expression (SAGE) and Expressed Sequence Tags (ESTs), which enabled largescale identification and characterization of expressed genes (Adams et al., 1992; Velculescu et al., 1995) [1,15] . The term "transcriptome" was later introduced by Velculescu et al.

(1997) [16] during studies on gene expression profiling in yeast. More recently, RNA sequencing (RNA-Seq) has become a powerful and widely adopted technique for studying gene expression patterns, identifying differentially expressed genes and understanding biological pathways associated with economically important traits. Comparative transcriptomic studies have been extensively used to investigate muscle development and growth-related traits in poultry species. Davis et al. (2015) [7] [17] reported that modern meat-type Pekin ducks and native Pekin ducks exhibit distinct differences in muscle growth and lipid metabolism despite sharing similar genetic backgrounds, suggesting that intensive selection has significantly shaped their transcriptomic profiles. Furthermore, Hu et al. (2021) [8] analysed skeletal muscle transcriptomes of Hanzhong Ma ducks during embryonic and postnatal developmental stages and identified several candidate genes associated with muscle growth and development. These findings highlight the importance of transcriptomic analysis for understanding breed-specific molecular mechanisms regulating skeletal muscle formation and metabolism.

Despite the availability of transcriptomic information in ducks, comparative studies focusing on molecular differences between commercial meat-type breeds and indigenous dual-purpose breeds remain limited.

Understanding these differences is important for developing sustainable breeding strategies aimed at improving productivity while maintaining adaptability and efficient resource utilisation. Therefore, the present study was undertaken to compare the leg muscle transcriptomes of Pekin and Hanzhong Ma ducks using publicly available RNA-Seq datasets retrieved from the National Center for Biotechnology Information database. Differential gene expression analysis and functional enrichment analysis were performed to identify genes and biological pathways associated with muscle growth, energy metabolism and adaptive regulation. The findings of this study may contribute to the identification of potential molecular markers for sustainable and efficient duck production systems.

Materials and Methods

Quality assessment of raw reads

The quality of raw sequencing reads was evaluated using FastQC (Andrews, 2010) [3] . Several quality parameters, including per-base sequence quality, GC content distribution, sequence duplication levels, adapter contamination and overrepresented sequences, were examined to assess the suitability of the sequencing data for downstream transcriptomic analysis.

Read preprocessing and trimming

Raw reads were subjected to preprocessing using fastp (v0.20.1) (Chen et al., 2018) [5] to remove low-quality reads, adapter sequences and ambiguous bases. Trimming procedures were carried out to improve the overall quality of the sequencing data before genome alignment. Only highquality filtered reads obtained after preprocessing were retained for subsequent analyses.

Reference genome preparation

The duck reference genome assembly GCF_015476345.1 and its corresponding annotation files were downloaded from the National Center for Biotechnology Information database. The reference genome index required for transcriptome alignment was constructed using HISAT2.

Alignment of sequencing reads

High-quality filtered reads from each sample were aligned against the duck reference genome using HISAT2. HISAT2 is a splice-aware aligner designed for efficient mapping of RNA-Seq reads to reference genomes and was used to generate alignment files in Sequence Alignment Map (SAM) format. The resulting SAM files were converted into Binary Alignment Map (BAM) format and sorted using SAMtools for efficient downstream analysis and transcript quantification.

Assessment of alignment quality

Alignment quality and mapping statistics were evaluated using SAMtools flagstat analysis. Parameters assessed included total reads, mapped reads, properly paired reads, singleton reads and reads aligned to different chromosomes. These analyses were performed to verify the reliability and overall quality of the transcriptome alignment.

Quantification of gene expression

Gene-level read counts were generated from the aligned BAM files using transcript annotation information corresponding to the duck reference genome. The generated count matrix was used for downstream differential gene expression analysis.

Differential gene expression analysis

Differential gene expression analysis between Pekin and Hanzhong Ma ducks was performed using the DESeq2 package (Love et al., 2014) [12] implemented in R statistical software. Read counts were normalized to account for differences in sequencing depth and library size among samples. For statistical comparison, Hanzhong Ma ducks were considered the reference group, while Pekin ducks were used as the comparison group. Differentially expressed genes (DEGs) were identified based on an adjusted p-value threshold of less than 0.01 and an absolute log₂ fold change greater than 2.

Functional enrichment analysis

Functional annotation and enrichment analyses of differentially expressed genes were carried out in DAVID (Huang et al., 2009) [9] to identify significantly enriched biological processes and metabolic pathways associated with skeletal muscle growth and metabolism. Gene Ontology (GO) analysis was performed to classify genes into biological process, molecular function and cellular component categories. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was conducted to identify metabolic and signalling pathways enriched among the differentially expressed genes.

Protein-protein interaction network analysis

Protein-protein interaction (PPI) network analysis was performed to investigate functional interactions among differentially expressed genes and to identify central regulatory genes within the interaction network. Networks were visualised in Cytoscape (v3.10.3) and hub genes were identified using the cytoHubba plugin based on the Maximal Clique Centrality (MCC) algorithm (Chin et al., 2014) [6] . Hub genes were identified based on their interaction connectivity and biological significance in muscle growth, metabolism and cellular regulation.

Data visualization

Visualization analyses including volcano plots, heatmaps and principal component analysis (PCA) plots were generated using R packages to evaluate transcriptomic variation and expression patterns between the two duck breeds. Volcano plots were used to visualize significantly differentially expressed genes, whereas PCA analysis was performed to assess clustering patterns among samples based on gene expression profiles.

Ethical statement

The present investigation was based entirely on publicly available RNA-Seq datasets obtained from the National Center for Biotechnology Information database. No live animals or biological samples were directly used in this study; therefore, separate institutional ethical approval was not required.

Results

RNA-Seq data quality assessment and alignment statistics

High-throughput RNA sequencing generated transcriptomic datasets from the leg muscle tissues of Pekin and Hanzhong Ma (HZM) ducks. Following quality assessment using FastQC and preprocessing with fastp, clean paired-end reads were aligned to the duck reference genome (GCF_015476345.1) using HISAT2. All four samples showed mapping percentages greater than 92%. The HZM sample SRR12701866 showed the highest mapping rate (97.06%), while Pekin sample SRR12702256 showed the lowest mapping rate (92.61%). Properly paired reads ranged from 83.03% to 85.43%. Alignment statistics for all samples are summarized in Table 1.

Identification of differentially expressed genes

Differential expression analysis was performed using DESeq2 with Hanzhong Ma ducks considered as the reference group and Pekin ducks as the comparison group. Genes with adjusted p-value <0.01 and |log₂FoldChange| >2 were considered significantly differentially expressed. A total of 3,398 differentially expressed genes (DEGs) were identified between Pekin and Hanzhong Ma ducks. Among these, 1,823 genes were significantly upregulated in Pekin ducks, whereas 1,523 genes were downregulated relative to Hanzhong Ma ducks. The volcano plot displayed the distribution of significantly upregulated and downregulated genes between the two breeds (Figure 1).

Protein-protein interaction network analysis of upregulated DEGs

A protein-protein interaction (PPI) network was constructed for 919 significantly upregulated DEGs using STRING with a high-confidence interaction score (>0.9). The network contained 919 nodes and 756 edges, with an average node degree of 1.65 and an average local clustering coefficient of 0.186. The observed number of edges (756) was higher than the expected number of edges (209), with a PPI enrichment p-value <1.0e -16 . The PPI interaction network of upregulated genes is presented in Figure 2.

GO and KEGG enrichment analysis of upregulated DEGs

Gene Ontology (GO) enrichment analysis of upregulated DEGs identified significant enrichment of biological process, cellular component and molecular function categories. The significantly enriched biological process terms included mitochondrial respiratory chain complex I assembly, mitochondrial electron transport from NADH to ubiquinone, proton motive force-driven ATP synthesis, tricarboxylic acid cycle, glycolytic process, reactive oxygen species metabolic process and fatty acid beta-oxidation. Cellular component enrichment included mitochondrion, mitochondrial inner membrane, mitochondrial matrix, mitochondrial large ribosomal subunit and protontransporting ATP synthase complex. The major molecular function categories included NADH dehydrogenase activity, proton transmembrane transporter activity and ATP synthase activity. KEGG pathway analysis identified oxidative phosphorylation, metabolic pathways, carbon metabolism, citrate cycle (TCA cycle), glycolysis/gluconeogenesis, pyruvate metabolism, PPAR signalling pathway, fatty acid degradation and muscle contraction pathways as significantly enriched pathways. GO and KEGG enrichment results for upregulated DEGs are presented in Figure 4a and4b.

Protein-protein interaction network analysis of downregulated DEGs

A PPI network was constructed for 1,067 significantly downregulated DEGs using STRING with a highconfidence interaction score (>0.9). The downregulated network consisted of 1,067 nodes and 328 edges, with an average node degree of 0.615 and an average local clustering coefficient of 0.166. The observed number of edges (328) exceeded the expected number of edges (78), with a PPI enrichment p-value <1.0e -16 . The PPI interaction network of downregulated genes is shown in Figure 5. The major hub genes identified among the downregulated DEGs included CDC45, MCM2, MCM5, MCM3, MCM10, ORC1, BUB1B, BUB1, KIF11, NUF2, MCM6, MAD2L1, NDC80, CENPI, TTK, CENPM, CENPK, CENPN, CENPU and CDK1. The top 20 hub genes identified from the downregulated network are presented in Figure 6.

GO and KEGG enrichment analysis of downregulated DEGs

GO enrichment analysis of downregulated DEGs identified significant enrichment in biological process, cellular component and molecular function categories. The major enriched biological process categories included extracellular matrix organization, cell division, mitotic cell cycle, chromosome segregation, microtubule-based movement, collagen fibril organization, spindle assembly checkpoint signalling, DNA replication initiation and nervous system development. Cellular component enrichment included extracellular matrix, neuron projection, glutamatergic synapse, spindle microtubules, condensin complex, kinetochore, chromosome centromeric regions and collagencontaining extracellular matrix. The significantly enriched molecular function categories included microtubule binding, microtubule motor activity, calcium ion binding, single-stranded DNA helicase activity, extracellular matrix structural constituent activity and metalloendopeptidase activity. KEGG pathway enrichment analysis identified motor proteins, cytoskeleton in muscle cells, cell cycle regulation, ECM-receptor interaction and neuroactive ligand-receptor interaction as significantly enriched pathways. GO and KEGG enrichment results for downregulated DEGs are shown in Figure 7a and7b.

Discussion

Several studies have investigated meat quality in many domesticated animals, including pig, cattle and chicken. The present study compared the leg muscle transcriptomic profiles of Pekin and Hanzhong Ma (HZM) ducks to identify molecular mechanisms associated with muscle growth and production-related traits. Differential expression analysis was performed using a stringent threshold of |log₂ fold change| > 2 and false discovery rate (FDR) < 0.05, following criteria adopted in previous transcriptomic studies (Luo et al., 2021) [13] . Using these parameters, a large number of significantly differentially expressed genes (DEGs) were identified between the two duck breeds, indicating substantial transcriptional divergence associated with their contrasting phenotypes. One of the major findings of this study was the strong enrichment of mitochondrial and oxidative metabolism-related genes in Pekin ducks. Upregulated genes were predominantly associated with oxidative phosphorylation, mitochondrial electron transport, ATP synthesis, glycolysis, tricarboxylic acid (TCA) cycle and fatty acid β-oxidation pathways. In particular, the top hub genes identified from the protein-protein interaction (PPI) network, including NDUFB9, NDUFB5, NDUFB4, NDUFB8, NDUFA6 and NDUFS7, were mainly components of mitochondrial respiratory chain complex I. These genes encode subunits of mitochondrial respiratory chain Complex I and are essential for electron transport and ATP synthesis during oxidative phosphorylation (Alston et al., 2018) [2] . The enrichment of these pathways indicates that Pekin ducks prioritize mitochondrial energy metabolism and oxidative ATP generation, which likely supports rapid muscle hypertrophy, higher feed conversion efficiency and superior meat yield. Previous transcriptomic studies in ducks have also identified genes associated with muscle growth, ATP generation, lipid metabolism and muscle organization, supporting the role of metabolic regulation in breed-specific growth performance (Wang et al., 2017) [17] .

In contrast, genes relatively more expressed in Hanzhong Ma ducks were enriched in pathways related to cell cycle regulation, DNA replication, extracellular matrix (ECM) organization, chromosome segregation, spindle assembly and neural development. Hub genes such as CDC45, MCM2, MCM5, ORC1, BUB1 and CDK1 are known regulators of mitotic progression and cellular proliferation. The enrichment of ECM organization, collagen fibril organization and cytoskeletal pathways suggests active tissue structuring and developmental regulation in Hanzhong Ma ducks. Neural-associated pathways, including axon guidance and nervous system development, were also significantly enriched. These transcriptomic characteristics may align with the indigenous adaptive features of Hanzhong Ma ducks, including structural resilience, slower but sustained growth and egg-laying performance. The comparative transcriptomic patterns observed between the two breeds reflect differences arising from long-term breeding objectives. Pekin ducks have undergone intensive artificial selection for rapid growth, muscle mass accumulation and commercial meat production, whereas Hanzhong Ma ducks have traditionally been maintained as indigenous dual-purpose ducks with emphasis on adaptability and reproductive performance. The present findings indicate that Pekin ducks rely predominantly on metabolic efficiency and mitochondrial energy production, while Hanzhong Ma ducks exhibit greater activity in pathways related to cellular organization and proliferative development. Protein-protein interaction analysis further highlighted the biological significance of the identified DEGs. The upregulated PPI network in Pekin ducks showed greater interaction density compared to the downregulated network, indicating coordinated regulation of metabolic pathways associated with oxidative phosphorylation and energy production. The identification of the top 20 hub genes provides important molecular targets that may be useful for future breeding and genomic selection programs. These genes could potentially be leveraged in molecular breeding strategies aimed at improving muscle growth efficiency, feed utilization, and adaptability in duck populations. Although the study identified biologically meaningful transcriptomic differences, several limitations should be acknowledged. First, the analysis was based on a relatively small sample size with only two biological replicates per breed. Consequently, the present findings should be considered exploratory rather than definitive. Several similar comparative transcriptome studies with sample size (n = 2) per breed was done (Karimi et al., 2024; King et al., 2026) [10,11] . A limited sample size reduces statistical power and may not fully represent within-breed variation. Second, the study relied on publicly available RNA-Seq datasets and factors such as environmental conditions, nutrition, management practices and genetic diversity could not be independently standardized. In addition, only leg muscle tissue from a single age group was examined, whereas gene expression patterns may vary across developmental stages and tissues. Despite these limitations, the study provides valuable transcriptomic insights into the molecular basis of breedspecific growth and muscle characteristics in ducks. The identification of mitochondrial metabolism-related genes in Pekin ducks and proliferation-and ECM-related genes in Hanzhong Ma ducks contributes to understanding the biological mechanisms underlying growth performance and adaptation. From a broader perspective, these findings may support sustainable poultry production strategies by identifying genetic markers associated with efficient resource utilization and breed adaptability.

Conclusions

The present study also highlights the potential relevance of integrating molecular breeding approaches with sustainable and eco-friendly farming practices. The identified hub genes and enriched pathways may serve as candidate markers for resource-efficient duck breeding programs, particularly for improving indigenous and regionally adapted duck populations in India. By combining transcriptomic insights with traditional farming systems and AYUSH-based agricultural values emphasizing sustainability and environmental harmony, future breeding programs may contribute toward improved poultry productivity, climate resilience and enhanced rural livelihoods.

References

1. Adams MD, Dubnick M, Kerlavage AR, Moreno R, Kelley JM, Utterback TR, et al. Sequence identification of 2,375 human brain genes. Nature. 1992;355(6361):632-634.
2. Alston CL, Heidler J, Dibley MG, Kremer LS, Taylor LS, Fratter C, et al. Bi-allelic mutations in NDUFA6 establish its role in early-onset isolated mitochondrial complex I deficiency. American Journal of Human Genetics. 2018;103:592-601.
3. Babraham Bioinformatics; 2010.
4. Buckingham M. Gene regulatory networks and cell lineages that underlie the formation of skeletal muscle. Proceedings of the National Academy of Sciences. 2017;114(23):5830-5837.
5. Chen S, Zhou Y, Chen Y, Gu J. Fastp: an ultra-fast allin-one FASTQ preprocessor. Bioinformatics. 2018;34:i884-i890.
6. Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: Identifying hub objects and sub-networks from complex interactome. BMC Systems Biology. 2014;8(Suppl 4):S11.
7. Davis RVN, Lamont SJ, Rothschild MF, Persia ME, Ashwell C, Schmidt CJ. Transcriptome analysis of posthatch breast muscle in legacy and modern broiler chickens reveals enrichment of several regulators of myogenic growth. PLOS ONE. 2015;10(3):e0122525.
8. Hu Z, Cao J, Zhang J, Ge L, Zhang H, Liu X. Skeletal muscle transcriptome analysis of Hanzhong Ma duck at different growth stages using RNA-Seq. Biomolecules. 2021;11(3):3152.
9. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nature Protocols. 2009;4(1):44-57.
10. Karimi P, Bakhtiarizadeh MR, Salehi A, Izadnia HR. Transcriptome analysis reveals the potential roles of long non-coding RNAs in feed efficiency of chicken. Scientific Reports. 2022;12(1):2558.
11. King JS, Seena TX, Thomas M, Rojan PM, Sukumaran J, Bindu KA. Comparative transcriptomic profiling of long non-coding RNAs reveals regulatory networks underlying growth divergence in White Leghorn and Black Australorp chickens. Gene Reports. 2026;43:102498.
12. Love M, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-Seq data with DESeq2. Genome Biology. 2014;15:550.
13. Luo L, Wang Q, Ma F. RNA-Seq transcriptome analysis of ileum in Taiping chicken supplemented with the dietary probiotic. Tropical Animal Health and Production. 2021;53(1):131.
14. Mitchell PO, Mills ST, Pavlath GK. Calcineurin differentially regulates maintenance and growth of phenotypically distinct muscles. American Journal of Physiology-Cell Physiology. 2002;282:C984-C992.
15. Velculescu VE, Zhang L, Vogelstein B, Kinzler KW. Serial analysis of gene expression. Science. 1995;270(5235):484-487.
16. Velculescu VE, Zhang L, Zhou W, Vogelstein B, Kinzler KW. Characterization of the yeast transcriptome. Cell. 1997;88(2):243-251.
17. Wang L, Li X, Ma J, Zhang Y, Zhang H. Integrating genome and transcriptome profiling for elucidating the mechanism of muscle growth and lipid deposition in Pekin ducks. Scientific Reports. 2017;7(1):3837.