The Ningdu yellow population used in this study comprised 399 purebred male chickens. All of the chickens were born and raised for 1 day by Jiangxi Nanshi Science and Technology Co., Ltd (Nanchang, ChinaNanchang, China). They were then randomly selected from the same batch and transferred to a farm in the city of Nanchang, ChinaNanchang, China. The chickens were fed with the same diet under a standardized feeding and management regimen, and given ad libitum access to water. At the age of 16 weeks, all of chickens were slaughtered at a commercial abattoir in Nanchang. After the Ningdu yellow chickens were slaughtered, the pituitary gland was carefully harvested for RNA isolation. The tissues were put into the sterile and frozen cryopreservation tubes and dipped into liquid nitrogen, and then conserved in −80°C ultra-freezer until RNA extraction. Moreover, a total of 17 growth and carcass traits were measured using “The Poultry Production Performance Terms and Measurement Statistics Method (NY/T823-2004),” including body length, live weight, carcass weight, carcass rate, semi-eviscerated weight, semi-eviscerated rate, eviscerated weight, eviscerated rate, breast muscle weight, breast muscle rate, chest circumference, average daily gain, chest depth, birth weight, back width, chest angle, and chest width. Each trait was measured by the same researcher.

A total of 12 pituitary glands of Chinese Ningdu yellow chickens exhibiting extreme growth and carcass phenotypes, with either a high live weight and carcass weight (H group) or a low live weight and carcass weight (L group) (Supplementary Table S1), were randomly divided into H and L groups with six chickens per group. The total RNA was isolated with TRIzol (Invitrogen, USA) by following the manufacture's instruction. The residual DNA was cleared away from total RNA by incubation with RNase-free DNase I (New England Biolabs, UK) at 37°C for 30 min. The quality of total RNA was assessed by a 2100 Bioanalyzer (Agilent, UK) and agarose gel electrophoresis. Three μg RNA per sample was used as input material for RNA preparation. Finally, 12 samples with RNA integrity number (RIN) values above eight were used for libraries construction. Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) by following manufacturer's recommendations, and index codes were added to attribute sequences to each sample. Briefly, mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer. First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (NaseH–). Subsequently, second strand cDNA synthesis was performed using DNA Polymerase I and RNase H. The remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext Adaptor with hairpin loop structure was ligated to prepare for hybridization. In order to select cDNA fragments of preferential ~200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA). Then 3 μl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min, followed by 5 min at 95°C before PCR. Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed using the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) according to the manufacturer's instructions. After cluster generation, the library preparations were sequenced on an Illumina Novaseq platform and 150 bp paired-end reads were generated.

The raw reads with N base contents more than 10% in a read and the percentage of low quality (Q ≤ 5) bases more than 50% were removed (14). The chicken reference genome Gallus_gallus-6.0/galGal6 was downloaded from genome website (http://asia.ensembl.org/Gallus_gallus/Info/Index) directly. The clean reads were mapped with TopHat (15) version 2.0.13. TopHat incorporates the Bowtie2 version 2.2.5 algorithm to perform the alignment (16). The aligned read files were processed by Cufflinks (17) version 2.2.1. Cufflinks uses the normalized RNA-Seq fragment counts to measure the relative abundances of transcripts. The unit of measurement is Fragments Per Kilobase of exon per Million fragments mapped (FPKM). Confidence intervals for FPKM estimates were calculated using a Bayesian inference method (18). Cuffdiff, a part of the Cufflinks package, takes the aligned reads from two or more conditions and reports genes and transcripts that are differentially expressed using a rigorous statistical analysis (19). Partial Least Squares-Discriminant analysis (PLS-DA) was performed to evaluate the whole gene transcripts between the H and L groups (20). Finally, Genes with a false discovery rate of <5% (i.e., FDR <0.05) and the difference ratio of FPKM between the H and L groups no less than two (|log2 fold change|≥ 1.0) as the significant threshold to identify differential expression.

Gene Ontology (GO) enrichment analysis and KEGG pathway analysis of the DEGs were implemented using DAVID Bioinformatics Resources (2021 Update) (21). Hypergeometric test was applied to map all DEGs to terms in GO (http://www.geneontology.org/GO.database.shtml) and KEGG pathway (http://www.genome.jp/kegg/pathway.html) database and search significantly enriched GO terms and KEGG pathway in DEGs comparing to the genome background. The FDR <0.05 was taken as a threshold of significance for GO terms and KEGG pathway.

cDNA was synthesized from 1 μg of total RNA using a PrimerScript RT reagent Kit (Takara, Japan) by following the standard manual. qRT-PCR was performed in a 10 μl reaction system containing 1 μl of 2.5-fold diluted cDNA, 5 μl of Power SYBR Green PCR Master Mix (Applied Biosystems, USA), 2 pM of each forward and reverse primer (Supplementary Table S2), and 3.6 μl water. PCR was conducted on a Bio-Rad CFX96 instrument (Bio-Rad) under the following cycling conditions: 10 min at 95°C, followed by 40 cycles at 95°C for 15 s and 60°C for 50 sec. GAPDH was used as the endogenous control. The quantification of gene expression levels was performed using the comparative Ct (2^−ΔΔCt) method. All assays were performed in triplicate and the average values were calculated.

Primer (Supplementary Table S2) was designed using Primer 5.0 software based on chicken PRKG2 gene sequences (SSC4: 45635577 – 45636527) and synthesized by GENERAY Biotech Co., Ltd (Shanghai, China). PCR was performed in a 30 μl reaction mixture containing 1 μl of DNA, 1 μl of forward and reverse primers (10 μM), 15 μl of 2 × Taq PCR MasterMix, and 12 μl ddH₂O. All amplifications included an initial denaturing step of 3 min at 95°C, followed by 35 cycles of 20 s at 94°C, 20 s at an optimized annealing temperature and 40 s at 72°C, with a final extension of 5 min at 72°C. The PCR products were sequenced by the GENERAY Biotech Co., Ltd (Shanghai, China). The ABI 3730xl DNA Sequencer was used for sequencing by the Sanger method. Long fragments were sequenced by bi-directional sequencing and then assembled using DNAStar software. SNPs were identified by comparative analysis of the complete PRKG2 sequence using Seqman software. Association analysis between the SNP and phenotypic traits was performed using the general linear model procedure in Plink v 1.07 (22). The slaughter batch was fitted as fixed effect in the model.

The associations between gene expression levels and the phenotypic values of growth and carcass traits were evaluated with Pearson correlation coefficient by using the cor function in the R software program. According to the previously reported QTT analyses in Drosophila melanogaster (23), pig (24) and human obesity (25), a conservative P < 0.0005 was set as the significance threshold.

In this study, 399 Ningdu yellow chickens were phenotyped across 17 growth and carcass traits. The summary statistic results for the 17 traits are shown in Table 1. Results indicated that most of the tested growth and carcass traits had low coefficients of variation, with the highest value of 25.01% for breast muscle weight and the lowest value of 5.86 for body length.

Correlational analyses were then conducted among 17 traits for growth and carcass traits (Figure 1, Supplementary Tables S3, S4). Results indicated that most growth and carcass traits had high correlations. More specifically, the carcass weight, semi-eviscerated weight, and breast muscle weight were all significantly, positively correlated with other growth and carcass traits. For example, the carcass weight was significantly, positively correlated with all other growth and carcass traits, especially for semi-eviscerated weight and eviscerated weight (r values were 0.95; P < 0.001).

To identify the genes associated with growth and carcass traits, genome-wide expression levels of genes in the pituitary gland were determined by RNA-seq for 12 Chinese Ningdu yellow chickens. After the adaptor sequences, ambiguous nucleotides, and low-quality sequences were removed, a total of 543,865,996 (81.60 Gb of data bulk) clean reads were generated from the 12 libraries and were divided into two groups with either a high live weight and carcass weight (H group) or a low live weight and carcass weight (L group). The summary of reads and the mapped genomic features are depicted in Supplementary Table S5. The results showed that about 85.59–90.94% of the clean reads mapped to the Gallus gallus reference genome, of which 1.35–1.75% have multiple alignments.

After the repeated and short-length sequences were eliminated, a total of 21,380 gene transcripts were aligned to the reference genome. An obvious shift in the whole gene transcripts was observed between the H and L groups (Figure 2). Specifically, the differential gene expression analysis identified a total of 428 (|log2 fold change| ≥ 1.0, FDR <0.05) DEGs between the H and L groups. In this study, the DEGs with higher expression levels in the H group as compared to the L group were denoted as “up-regulated,” while those with lower expression levels in the H group as compared to the L group were denoted as “down-regulated.” The expression levels of 110 among 428 genes were determined to be up-regulated in the H group, while the other 318 genes were down-regulated (Figure 3, Supplementary Table S6). For example, SLCO1C1 (log2 fold change = 3.68, FDR = 3.62E-04), EGR4 (log2 fold change = 3.30, FDR = 9.86E-03), and RLN3 (log2 fold change = 2.07, FDR = 4.46E-02) were found to be significantly upregulated; and MC1R (log2 fold change = −4.53, FDR = 2.43E-03), MYB (log2 fold change = −3.43, FDR = 1.11E-03), and IRF4 (log2 fold change = −3.40, FDR = 5.88E-06) were found to be significantly downregulated.

The NCBI web-based functional annotation tool DAVID Bioinformatics Resources (2021 Update) was used to investigate the functional associations of the changes in gene expression identified in the pituitary glands between the two groups (21, 26). GO biological process, molecular function, and cellular component were selected as the annotation categories for clustering. The software identifies enriched ontologies for a particular gene list and clusters statistically significant terms for their constituent genes. The gene list used in this analysis contained the 428 DEGs. The results showed that these DEGs were enriched (FDR <0.05) in many more GO terms than expected (Supplementary Table S7), including 134 terms significantly enriched in GO biological process and 15 terms significantly enriched in GO cellular component. For example, immune system process, immune response, and response to stimulus were enriched in GO biological process (Figure 4A), and plasma membrane part, plasma membrane, and cell periphery were enriched in GO cellular component (Figure 4B). Specifically, nine terms were significantly enriched in GO molecular function. Such as, receptor activity, molecular transducer activity, and transmembrane receptor activity (Figure 4C).

Furthermore, the KEGG Pathway Database, which documents many types of influences and interactions, was used to help us generate an overview of the network pathway. The results showed that neuroactive ligand-receptor interaction including 14 genes, FDR = 6.00E-03, and cytokine-cytokine receptor interaction including 10 genes, FDR = 8.10E-03 were the only two terms significantly enriched in the KEGG pathway. The DEGs enriched in two KEGG pathways are shown in Supplementary Table S8. Specifically, the differentially expressed GH (growth hormone) gene and GHRHR (growth hormone releasing hormone receptor) gene, which are related to growth, were included in neuroactive ligand-receptor interaction pathway.

To further validate the DEGs result, the expression patterns of six randomly selected DEGs were quantified by qRT–PCR. The six genes were GH, IRF4 (interferon regulatory factor 4), NPBWR2 (neuropeptides B and W receptor 2), NRROS (negative regulator of reactive oxygen species), PRKG2 (protein kinase cGMP-dependent 2), and ENSGALG00000051068. A total of 12 RNA samples (six for each group) used for RNA sequencing analysis, were also used for qRT–PCR analysis. The results showed that all six genes were differentially expressed between the H and L groups, supporting a high reproducibility of RNA-seq results (Figure 5).

To identify the candidate genes for chicken growth and carcass traits from the DEGs, the QTT analysis was performed to test the correlation between gene expression levels and live weight or carcass weight values. Seven DEGs were identified to be significantly associated with live weight and three were identified to be significantly associated with carcass weight (Table 2). For example, two were positively correlated (GRIA1 (r = 0.92, P = 1.81E-05) and PRKG2 (r = 0.85, P = 4.82E-04)) and five were negatively correlated [i.e., NRROS (r = −0.91, P = 5.00E-05) and IRF4 (r = −0.85, P = 4.31E-04)] with live weight. Furthermore, the eight DEGs were further investigated to determine if growth-related phenotypes or diseases had been reported in their corresponding knockout mice using the Mouse Genome Informatics database (http://www.informatics.jax.org/). For three of the candidate genes, GRIA1, IRF4, and PRKG2, the results supported a role for these genes in growth and carcass traits (Table 2).

We next investigated the variants of the three candidate genes through RNA-seq data using GATK (v4.1.1.0) software. A total of 39 polymorphisms (35 SNPs and four indels) were identified in the 12 Chinese Ningdu yellow chickens. Of these mutations, the genotypes of the SNP rs16400745 was the only SNP consistent with the genotypic patterns of the phenotypic values of the aforementioned 12 chickens, which was considered to the candidate SNP. The SNP rs16400745 was 3 prime UTR variant in the PRKG2 gene. We further genotyped the SNP rs16400745 in a total of 399 Chinese Ningdu yellow chickens using the Sanger method (Primer: FP: TAAAGACTCCGAAACTCACT, RP: ACGCACCATAGACTCATT). The results showed that the SNP rs16400745 had three genotypes. The frequencies of CC, CT, and TT in rs16400745 were 0.40, 0.50, and 0.10, respectively.

An association analysis between the SNP rs16400745 and 17 growth and carcass traits was performed, and the results showed that the SNP rs16400745 was significantly associated with carcass weight (P = 9.68E-06), eviscerated weight (P = 3.04E-05), semi-eviscerated weight (P = 2.14E-04), breast muscle weight (P = 2.81E-04), live weight (P = 4.34E-04), body length (P = 4.95E-03), chest circumference (P = 8.33E-03), and average daily gain (P = 8.80E-03). The effects of the SNP rs16400745 on growth and carcass traits are illustrated in Table 3. For example, the SNP rs16400745 (T allele) was associated with an increased carcass weight and live weight by 37.50 and 33.00 g, respectively.