All experimental procedures were performed according to the Guide for Animal Care and Use of Laboratory Animals of the Institutional Animal Care and Use Committee of Ningxia University. The animal use protocol listed below has been reviewed and approved by the Animal Ethical and Welfare Committee (NXUC3127858).

The Jingyuan chickens used in this study were from the Zhaona Chicken Breeding Center in Pengyang County, Ningxia. All the chickens were grown under the same management conditions and given the same feeding. Thereafter, 45 chickens (15 males and 30 females) were slaughtered when they were 180 days old. Tissue samples were rapidly collected from the breast muscles of the slaughtered chickens; these tissues were immediately snap-frozen in liquid nitrogen tanks and stored in the laboratory in −80°C and −20°C refrigerators for subsequent experiments.

Shear force was determined according to the Warner-Bratzler method [12]. For this purpose, the meat was placed in a constant temperature water bath at 80°C and when the temperature of the center of the meat sample reached 75°C, it was removed and cooled down to the room temperature. Samples of 1.0 cm×1.0 cm×3.0 cm were cut perpendicular to the direction of the muscle fibers with a sampler and the tenderness of these meat samples was measured using an RH-N50 tenderness tester (Tenovo, Beijing, China). The same process was repeated three times for each sample, and the average value was noted. The contents of water holding capacity (WHC) and cooking loss were determined using the method described by Dalgaard et al [12]. The WHC of the meat samples was measured by Bulader-M10 meat water loss tester (Bulader, Beijing, China). A meat sample of 1.0 cm in diameter and 1.0 cm thickness was excised from the breast muscle in the direction of the muscle fibers with a sharp circular sampler. The initial weight of the sample was marked as W1. Thereafter, it was placed between two layers of gauze with 20 layers of filter paper on top and bottom and pressurized to 35 kg for 5 min. The sample was again weighed, and this weight was marked as W2. The process was repeated three times for each sample was repeatedly measured three times and the average value was noted. The WHC of the meat samples was expressed as a measure of the water loss as follows: water loss (%) = (W1–W2)/W1×100. To determine the cooking loss, the breast muscle tissue was cut into 2.0 cm×2.0 cm×3.0 cm meat samples; the fat and connective tissues were removed from the surface of the muscle and the initial weight of the sample was marked as Z1. Thereafter, the tissue samples were placed in boiling water in an aluminum steamer for 30 min, kept in a cool ventilated place for 15 min, and again weighed; this weight was marked as Z2. The same process was repeated three times for each sample, and the average value was noted. The cooking loss was expressed as, cooking loss (%) = (Z1–Z2)/Z1×100.

The contents of IMP and inosine were determined using the method described by Jung et al [7]. Breast muscle samples, 10 g, were ground in a meat grinder and each sample was divided into three portions weighing 2 g, accurate to 0.0001 g. The samples were placed in a 15 mL centrifuge tube with 5 mL of 6% perchloric acid. Subsequently, these samples were homogenized at high speed, followed by centrifugation of the homogenate at 8,000 rpm for 10 min; the centrifuged samples were filtered into a 100 mL volumetric flask. The filter residues were again homogenized with 5 mL of 6% perchloric acid, and the homogenate was centrifuged and filtered. The two supernatants were thus combined and adjusted to pH 6.5 with 5.0 mol/L and 0.5 mol/L NaOH, and then mixed well. This sample was then filtered through a filter screen, and the filtrate was collected for high-performance liquid chromatography analysis.

The chromatographic conditions were as follows: Waters Atlantis dC18 chromatographic column (5 μm, 4.6 mm×150.0 mm) was used and the column temperature was 25°C. The mobile phase consisted of a 15:85 mixture of liquid A and liquid B, where liquid A was methanol and liquid B was 0.05 mol/L of a 1:1 mixture of dipotassium hydrogen phosphate and potassium dihydrogen phosphate (liquid B was prepared on the spot, filtered by 0.45 μm membrane, and degassed by ultrasonication). The flow rate was 1 mL/min, injection volume was 10 μL, detection wavelength was 250 nm.

The Trizol method was used to extract the total RNA from nine breast muscle tissue samples excised from the Jingyuan chickens. The total RNA was tested for integrity by 1% gel electrophoresis, and the concentration and purity of the RNA samples were detected by Nanodrop ND-2000 (Therom, Waltham, MA, USA) and Agilent 2100 bioanalyzer (Agilent, Santa Clara, CA, USA). All the RNA samples with RNA integrity numbers >8.0 and absorbance 260:280 ratios >1.9 were selected for library construction and deep sequencing. A total amount of 10 μg RNA per sample was used as input material for the RNA sample preparations. Sequencing libraries were generated using NEB Next UltraTM RNA Library Prep Kit for Illumina (NEB, Beijing, China) and following the procedures and standards of the kit’s manual. For this process, mRNA was purified from the total RNA using Oligo (dt) magnetic beads, followed by fragmentation of the mRNA under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5×). First strand cDNA was synthesized using a random hexamer primer and reverse transcriptase. Subsequently, second strand cDNA synthesis was performed using DNA polymerase I (DNA-pol I) and ribonuclease H (RNase H). After poly-A of the 3′ ends of the DNA fragments, NEBNext adapter with a hairpin loop structure was ligated to prepare the DNA for hybridization. AMPure XP beads (Beckman Coulter, Beverly, MA, USA) were used to select cDNAs around 350 to 400 bp in length. Subsequently, the second strand of U-containing cDNA was degraded using the USER enzyme (NEB, China), and finally the cDNA samples were amplified by polymerase chain reaction (PCR). The cDNA library quality was assessed with an Agilent 2100 bioanalyzer. Finally, 9 libraries were sequenced using an Illumina HiSeqTM2500 (Illumina, San Diego, CA, USA).

To get clean reads, the raw sequencing data were quality controlled using FastQC software, that is, clean reads were obtained after removing the reads containing adapters, low quality reads, and reads containing poly-N ≥10%. The Q20, Q30, and content of guanine and cytosine contents of the clean reads were calculated based on the abovementioned criteria. The clean reads that were obtained after the quality control process were mapped to the chicken reference genome via Hisat2 and Bowtie2 (ftp://ftp.ensembl.org/pub/current_fasta/gallus_gallus/dna/Gallus_gallus.Gallus_gallus-5.0.dna.toplevel.fa.gz). The gene expression level was estimated by the number of normalized fragments per kilobase of transcript per million fragments (FPKM) method. Differentially expressed genes (DEGs) among the samples were analyzed using the edgeR package. Genes that had an adjusted p-value <0.05 and a fold change (|log₂FC|)>1 were considered as DEGs.

The obtained DEGs required an insight into their biological functions to uncover important regulatory factors; this is usually analyzed using gene ontology (GO) (http://www.geneontology.org/) in bioinformatics analysis [13]. GO functional enrichment analysis can locate gene functions in three aspects: cellular components, biological processes, and molecular functions. Transcripts with similar annotations can be integrated and categorized using the GO annotation enrichment analysis. Biological pathway enrichment analysis relies heavily on the Kyoto encyclopedia of genes and genomes database (KEGG) ( http://www.genome.jp/kegg/) [13]. Pathway enrichment analysis can pinpoint the major metabolic transductions and signaling pathways of the DEGs. GO and KEGG enrichment analyses were performed using the GOseq based on Wallenius non-central hyper-geometric distribution and clusterprofiler R package, with a setting of p<0.05 as the significant enrichment threshold.

Protein interaction networks were analyzed by String protein interaction database ( http://string-db.org/) for proteins encoded by DEGs and edited visually using the Cytoscape v3.7.0 software [14].

Gene set enrichment analysis, construction of molecular transformation networks, and metabolic mapping were carried out through the Reactome pathway database ( https://reactome.org/) [15]. Network cluster analysis was performed by the String protein interaction database.

Based on the published gene sequences of ethanolamine-phosphate phospho-lyase (ETNPPL), acyl-CoA synthetase long-chain family member 4 (ACSL4), phosphodiesterase 10A (PDE10A), glutamate-ammonia ligase (GLUL), ATP synthase F1 subunit alpha W chromosome (ATP5A1W), ATP synthase F1 subunit alpha Z chromosome (ATP5A1), histidine triad nucleotide binding protein W (HINTW), and beta-actin (β-actin) on GeneBank, primers were designed using the primer 5.0 software (Table 1) and sent to the Sangon Company for synthesis. Total RNA was extracted from the pectoral muscle samples, using β-actin as the housekeeping gene. cDNA was synthesized using PrimeScript RT reagent Kit with gDNA Eraser (Takara, Dalian, China). Quantitative Real-time PCR (qRT-PCR) amplification was performed using SYBR premix Ex TaqTM II kit (Takara, China). Three replicates were performed for each sample. The relative gene expression levels were calculated using the 2^−ΔΔCt method [16].

One-way analysis of variance (ANOVA) was performed using IBM SPASS Statistics 23.0 software, and the results were expressed as mean±standard error; the correlation of IMP content with shear force and the relative gene expression was assessed using bivariate correlation analysis.

Shear force, WHC, and cooking loss were calculated for the tissue samples excised from the breast muscles of Jingyuan chickens reared under the different feeding management conditions. A one-way ANOVA of these results showed that the breast muscle shear force in caged chickens was significantly higher than that in free-range chickens (p<0.05), while the cooking loss in meat of caged chickens was extremely significantly higher than that in the meat of free-range chickens (p<0.01). With respect to the different sex of the Jingyuan chickens, breast muscle shear force was extremely significantly higher in hens than in cocks (p<0.01), while WHC was significantly higher in hens than in cocks (p<0.05) (Table 2).

IMP and inosine content were measured in the tissue samples excised from the breast muscles of Jingyuan chickens reared under the different feeding management conditions. A one-way ANOVA of these results showed that the breast muscle IMP content in free-range chickens was extremely significantly higher than that of caged chickens (p<0.01), while the inosine content in caged chickens was extremely significantly higher than that of free-range chickens (p<0.01). Moreover, with respect to the different sex of the Jingyuan chickens, breast muscle IMP content in hens was higher than that in cocks, while the breast muscle inosine content in cocks was higher than that in hens, but neither of these were significantly different (Table 3).

Bivariate correlation analysis of the IMP content in the breast muscle tissue samples and the shear force of the breast muscles showed that the correlation in case of caged hens is R² = −0.278 (p>0.05), in case of free-range hens it is R² = −0.522 (p<0.05), and in case of caged cocks it is R² = −0.38 (p>0.05). The scatter plot (regression line) of IMP content in breast muscle tissues with shear force of breast muscles visually reflects the relationship (Figure 1). The IMP content in the breast muscle tissues of free-range hens was significantly negatively correlated with the breast muscle shear force (p<0.05), and the breast muscle shear force was negatively correlated with the IMP content in breast muscles tissues of caged hens and cocks.

Sequencing of cDNA libraries from nine breast muscle tissue samples showed over 40.12 million raw reads per sample sequenced. After the quality control of the raw reads and the subsequent removal of the low quality reads, over 38.94 million clean reads were obtained per sample. The nucleobase content distribution showed that the GC content was stable between 55.86% and 56.72%, and the nucleobase composition was stable and balanced. The number of clean reads obtained per sample was above 97% for Q20 and above 93% for Q30. More than 70% of the clean reads mapped to the chicken reference genome and over 65% map to unique loci. Furthermore, to ensure the accuracy of sequencing results and exclude errors caused by abnormal samples, the pearson product moment correlation coefficient was evaluated among the samples. The biological replicate R²≥0.85, the expression pattern similarity between the samples was high, and the RNA-seq results were reliable and could be used for subsequent analyses.

In this study, the expression levels of genes were quantified by FPKM values, and 149 DEGs were identified between the two groups Jingyuan chickens of different feeding patterns; 60 of these genes were upregulated, and 89 genes were downregulated (Figure 2A). There were 107 DEGs identified between the two groups of different sexes of Jingyuan chickens, between the two groups Jingyuan chickens of different sexes, 53 of these genes were upregulated and 54 genes were downregulated (Figure 2B). Clustering analysis of the DEGs helped to visualize the expression levels of these DEGs among the samples, and it highlighted the reproducibility and credibility of the data. The results of the clustering analysis of the DEGs between Jingyuan chickens of different feeding patterns and different sexes groups were good (Figure 3A, B).

To understand the biological functions involved with the DEGs, GO, and KEGG signaling pathway enrichment analyses were performed. The DEGs among the Jingyuan chickens of different feeding patterns were mainly enriched in 993 GO terms in three major categories, namely biological processes, molecular functions, and cellular components, out of which 105 terms were significantly enriched (p<0.05) (Figure 4A); these included purine nucleoside metabolic process and glutamine family amino acid metabolic process (Table 4). KEGG signaling pathway enrichment analysis results showed that the DEGs were enriched in alanine (Ala), aspartate (Asp), and glutamate (Glu) metabolism as well as purine metabolism. The DEGs among the Jingyuan chickens of different sexes were mainly enriched in 696 GO terms in three major categories, namely biological processes, molecular functions, and cellular components, out of which 83 terms were significantly enriched (p<0.05) (Figure 4B). However, in this case, the purine nucleoside monophosphate metabolic process was not significantly enriched (Table 4), and the KEGG signaling pathway enrichment analysis results showed that the DEGs were not enriched in the IMP-related pathways.

In transcriptional regulation, protein network regulatory analysis can identify the key candidate genes that potentially regulate important traits. Therefore, we combined the results of GO and KEGG signaling pathway enrichment analyses to further construct protein network interaction analysis maps for the screened key genes PDE10A and GLUL (Figure 5). The results of the protein network interactions of these two DEGs showed that PDE10A interacted with seven genes in purine metabolism (confidence 0.700), including the ADSL gene in the IMP de novo synthesis pathway (Figure 6A-J), while GLUL interacted with eight genes in glutamine metabolism (confidence 0.700), including PPAT in the IMP de novo synthesis pathway. Therefore, the differential deposition of IMP in the breast muscles of Jingyuan chickens with different feeding patterns may be related to the expressions of GLUL and PDE10A genes.

Reactome pathway analysis showed that the PDE10A gene is involved in eight reactome pathways, including the purine salvage. GLUL gene is involved in five reactome pathways, including the synthesis, interconversion, and transamination of amino acid as well as the metabolism of amino acids and their derivatives. The purine salvage involves IMP synthesis (Figure 7), while the synthesis and interconversion of amino acids as well as the metabolism of amino acids and their derivatives include the metabolism of Asp, asparagine (Asn), Glu, and glutamine (Gln). Network cluster analysis (confidence 0.900) showed that the PDE10A gene is involved in nucleotide metabolism, IMP metabolism, and nicotinate and nicotinamide metabolism (Table 5). On the contrary, the GLUL geng encoded protein is involved in glutamine synthetase, carbon metabolism, amino sugar, and nucleotide sugar metabolism as well as the metabolism and synthesis of many amino acids (Table 6).

Analysis of the qRT-PCR validation results for the seven genes showed that the melting curves of all genes were single-peaked, the primers were well-designed, and the validation results were consistent with the trend of RNA-seq results (Figure 8). Hence, the sequencing results were reliable and could be used for subsequent gene function validation.

Bivariate correlation analysis between the PDE10A gene expression and breast muscle IMP content in free-range chickens showed the correlation coefficient, R² = 0.616 (p<0.05). A similar analysis for caged chickens showed that the correlation coefficient, R² = 0.567 (p<0.05). Similarly, the correlation coefficient, R² = 0.509 (p>0.05) for the correlation between the GLUL gene expression and breast muscle IMP content in free-range chickens, and in caged chickens it is R² = 0.44 (p>0.05). The scatter plot (regression line) was drawn to visualize the relationship between the IMP content and gene expression. (Figure 9). PDE10A gene expression was significantly and positively correlated with the IMP content in breast muscle tissues (p<0.05), and GLUL gene expression was positively correlated with the IMP content in the breast muscle tissues.