Exploring Genomic Variants Related to Residual Feed Intake in Local and Commercial Chickens by Whole Genomic Resequencing

Improving feed efficiency is a major goal in poultry production to reduce production costs and increase profitability. The genomic variants and possible molecular mechanisms responsible for residual feed intake (RFI) in chickens, however, remain poorly understood. In this study, using both local and commercial breeds, genome re-sequencing of low RFI and high RFI chickens was performed to elucidate the genomic variants underlying RFI. Results showed that 8,505,214 and 8,479,041 single nucleotide polymorphisms (SNPs) were detected in low and high RFI Beijing-You chickens, respectively; 8,352,008 and 8,372,769 SNPs were detected in low- and high-RFI Cobb chickens, respectively. Through a series of filtering processes, 3746 candidate SNPs assigned to 1137 genes in Beijing-You chickens and 575 candidate SNPs (448 genes) in Cobb chickens were found. The validation of the selected 191 SNPs showed that 46 SNPs were significantly associated with the RFI in an independent population of 779 Cobb chickens, suggesting that the method of screening associated SNPs with whole genome sequencing (WGS) strategy was reasonable. Functions annotation of RFI-related genes indicated that genes in Beijing-You were enriched in lipid and carbohydrate metabolism, as well as the phosphatase and tensin homolog (PTEN) signaling pathway. In Cobb, however, RFI-related genes were enriched in the feed behavior process and cAMP responsive element binding protein (CREB) signaling pathway. For both breeds, organismal development physiological processes were enriched. Correspondingly, NOS1, PHKG1, NEU3 and PIP5K1B were differentially expressed in Beijing-You, while CDC42, CSK, PIK3R3, CAMK4 and PLCB4 were differentially expressed in Cobb, suggesting that these might be key genes that contribute to RFI. The results of the present study identified numerous novel SNPs for RFI, which provide candidate biomarkers for use in the genetic selection for RFI. The study has improved knowledge of the genomic variants and potential biological pathways underlying RFI in chickens.


Introduction
As with most animal production systems, the cost of feed in chicken meat production is a high proportion of the total farming expense, being nearly 70% of the total cost of poultry production [1]. The selection of more efficient animals reduces production costs, decreases the area of land required for feed production and reduces environmental impact of nitrogen and phosphorus excretion resulting from the digestion/metabolic processes.
Residual feed intake (RFI), defined as the difference between the actual feed intake (FI) and predicted requirements based on growth and maintenance, was first proposed as an alternate measure of feed efficiency (FE) by Koch et al. [2]. RFI is a sensitive and accurate indicator that is increasingly accepted as an alternative measure of FE of livestock. RFI has a moderate heritability (0.21 to 0.49) in broilers and responds to selection [3]. Elucidation of the genomic variants for RFI could potentially be important for marker-or gene-assisted selection [4].
Current whole genomic sequencing methodologies are important tools to unravel the mechanisms of complex traits, facilitate new understanding of the genetic regulation of phenotype and allow the identification of potential biomarkers for early or more accurate genetic prediction [5]. In recent years, some progress has been gained in attempts to identify the molecular mechanisms underlying RFI. Transcriptome analysis of mRNA and microRNA in skeletal muscle in pigs showed that genes involved in mitochondrial energy metabolism were downregulated, while those involved in skeletal muscle differentiation and proliferation were upregulated in low RFI (LRFI), as compared to high RFI (HRFI) animals [6]. The duodenal transcriptome architecture of extreme RFI phenotypes of six brown-egg dwarf hens showed that significant differentially-expressed genes were involved in several specific biological functions related to the processes of digestibility, metabolism, biosynthesis and energy homeostasis [7]. Combination analysis of genome-wide association and transcriptome sequencing in chickens showed that the effective single-nucleotide polymorphisms (SNPs) related to energy utilization were located in a 1-Mb region (16.3-17.3 Mb) of chicken chromosome 12, and liver transcriptomes may contribute to RFI, which is directly influenced by appetite, cell activities and fat metabolism [8]. Genome-level analyses of genetic markers and candidate genes for RFI in pigs by genome-wide associations and systems genetic analyses revealed that the SNPs of MAP3K5, PEX7 and DSCAM might be interesting markers for RFI. Functional annotation of genes indicated that regulation of protein and lipid metabolism, gap junction formation, inositol phosphate metabolism and the insulin signaling pathway were significant biological processes and pathways associated with RFI [4]. Genome-wide association analysis of FI and RFI in Nellore cattle identified several markers located on chromosomes 4, 8, 14 and 21 in regions near genes regulating appetite and ion transport [9].
There are more than one hundred domestic chicken breeds available across China, and collectively, they have roughly the same market share in China as do introduced commercial lines. Beijing-You chickens, one of the representative domestic chickens, are known for their peculiar characteristics including higher intramuscular fat content and chewy texture [10], which are appealing palatability factors for Chinese consumers. The Cobb line has been intensely selected for high feed efficiency and growth rate and is an important commercial breed world-wide. To gain insights into the molecular alterations underlying differences between HRFI and LRFI chickens, whole genome sequencing (WGS) was conducted using HRFI and LRFI birds of local and commercial breeds.

Ethics Statement
All of the animal experiments were conducted in accordance with the Guidelines for Experimental Animals established by the Ministry of Science and Technology (Beijing, China). Animal experiments were approved by the Science Research Department (in charge of animal welfare issue) of the Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (IAS-CAAS) (Beijing, China).
Ethical approval on animal survival was given by the animal ethics committee of IAS-CAAS (No. IASCAAS-AE-03).

Animals and Sample Collection
A population of 400 Beijing-You chickens (200 males and 200 females) with a similar body weight, aged 70 days, was selected from the base population of the Institute of Animal Sciences (CAAS), and these birds came from 75 male and 153 female families (Table 1) where b 0 is the intercept and b 1 and b 2 are partial regression coefficients of FI on BW 0.75 and WG, respectively. RFI values were generated by the regression procedure using SAS 8.0 for Windows statistical software (SAS Institute, Cary, NC, USA). Genomic DNA was isolated from whole blood samples collected from 48 HRFI and 48 LRFI chickens (the details of the group information are shown in Table 1), within each breed, using the phenol-chloroform method. DNA quality was determined by gel electrophoresis, and the 48 individual DNA samples from each group were used to construct three pools, each of 16 samples. For Beijing-You chickens, each pool consisted of 8 males and 8 females. In brief, there were three pools and 16 chickens of each pool for each sub-group. Paired-end sequencing libraries were constructed using the Nextera DNA Library Preparation Kit (Illumina Inc., San Diego, CA, USA) according to the manufacturer's standard protocol. At the end of the feeding studies (Day 98), birds were killed, and weight percentages of breast and thigh muscles from the two groups (n = 48, HRFI and LRFI, in Beijing-You chickens) were measured by standard dissection procedures. Breast muscle and liver samples from subsamples (n = 8) of HRFI and LRFI male Beijing-You chickens were cut into small pieces, rapidly snap-frozen and stored at −80 • C for quantitative real-time PCR (qRT-PCR) analyses. Likewise, samples of breast muscle and hypothalamus from HRFI and LRFI Cobb chickens (each n = 8) were collected and saved the same way.

Genome Sequence Assembly and Data Analysis
All libraries were sequenced on the Hiseq 2500 platform (Illumina) instrument with 125-bp paired-end reads. Reads of extremely low quality (>10 consecutive nucleotides with Phred scores < 10), with adaptor contamination or without a quality control-passed paired read were discarded using NGSQC toolkit (v2.3.3), resulting in a qualified clean dataset. More than 20× coverages were determined for each pool ( Table 1). All clean sequenced reads were mapped to a reference genome (Galgal 4.0) using the BWA tool (v0.7.10) [11] with default parameters. Output SAM files were converted to BAM files with command Samtools (v0. 1.1.18), and all converted BAM files were sorted with command Samtools (v0.1.1.18). PCR duplications were removed with the "rmdup" argument in Samtools (v0.1.18) [12] after mapping. SNPs were identified and genotyped for each sample pool with the mpileup function in Samtools (v0.1.18), as well as the VarScan.jar tool (v2.3.7) written in Java [13]. Only highly confident variants supported by both methods were used for further downstream analyses.
The sequence data reported in this paper have been deposited in the Genome Sequence Archive of BIG Data Center, Beijing Institute of Genomics (BIG), Chinese Academy of Sciences and is publicly accessible at http://gsa.big.ac.cn/preview/preview.action?code=WkQSr8M0 (accession no. PRJCA000311).

SNP Detection and Function Annotation
To identify SNPs with frequencies that consistently differed across replicates (pools) of these divergent populations, the divergence of allelic frequency between two phenotypically-extreme groups in the two breeds of chickens was statistically tested using analysis of variance and F-tests. SNPs with Benjamini-adjusted p-values of <0.05 were analyzed further. Pairwise fixation index (F ST ) values were calculated using the F ST -test tool in the PoPoolation2 package. SNPs with F ST values in the top 5% and allelic frequency changes greater than 35% (greater than the median of allelic frequency divergence in the top 5% SNPs) were identified as trait-related SNPs. Genes harboring these SNPs (50 kb of upstream and downstream of the SNPs) were selected for pathway and function enrichment analyses using the Ingenuity Pathway Analysis (IPA) system (http://www.ingenuity.com) according to the default parameters. Functional and canonical pathway analyses were used to identify significant biological functions and pathways; associated networks were also constructed.

SNP Validation by Sanger Sequencing
Nine randomly-chosen SNPs were subjected to validation by PCR analysis and Sanger sequencing using 16 phenotypically-divergent chickens. Targeted sequences were first PCR-amplified using 50 ng of genomic DNA with a Taq PCR MasterMix (Tiangen Biotech Co., Ltd., Beijing, China). Verification of PCR was performed by 1% agarose gel electrophoresis. Amplicons were then purified according to the EasyPure PCR Purification Kit's standard protocol (Transgen Biotech Co., Ltd., Beijing, China) and sent to Tianyi Huiyuan (Beijing, China) for Sanger sequencing using primers described in Table S1. Results were analyzed using DNAMAN 8.0 (Lynnon Biosoft, San Ramon, CA, USA) (http://www. lynnon.com/), and ratios of bases occurring at SNP locations were recorded.

Association Analysis of the SNPs from Whole Genome Sequencing in the Validation Population of Cobb
In the validation population, there were 779 Cobbs (345 males and 434 females), which were the offspring of the generation for WGS bird screening. The RFIs of these chickens were measured by the methods described above. One hundred ninety SNPs were randomly selected from the candidate SNPs obtained from the WGS of Cobb. Genotyping was done on Axiom ® SNP arrays using the Affymetrix GeneTitan ® system according to the procedure described by Axiom 2.0 Target Prep 384 Samples Protocol (Affymetrix; Santa Clara, CA, USA) (https://www.thermofisher.com/). Genetic analyses were performed using ASREML [14] packages under the R v3.1.1 environment [15]. Breeding values of Genes 2018, 9, 57 5 of 18 RFI were estimated using the BLUP based on an animal model with a relationship matrix. The models used in matrix notation were as follows: where y is the vector of observations; b, the vector of fixed effects of sex; a, the vector of random direct genetic effects; e, the vector of random residual effects; X and Z are incidence matrices relating the observations to the respective fixed and direct genetic effects.
The association analyses of estimated breeding values (EBVs) and SNPs were conducted using the Wald test under the R v3.1.1 environment [15]. The associations between SNP marker genotypes and RFI were estimated by using Proc GLM in SAS 9.4 (SAS Institute Inc., Cary, NC, USA). Additive genetic effects were estimated by pair-wise comparison of the two homozygous genotypes, and the dominance effects were calculated as the deviation of the heterozygote effect from the average of the two homozygous genotypes.

Quantitative Reverse Transcription-PCR Validation of Expression of Candidate Genes
To investigate whether the expression profiles of the variant genes differed between HRFI and LRFI groups, key genes (Table S2) that were enriched in the canonical pathways were chosen for qRT-PCR analyses. Total RNA was extracted from hypothalamus, breast muscle and liver using the RNAsimple Total RNA Kit (Tiangen) according to the manufacturer's recommendations. For complementary DNA (cDNA) synthesis, 1 µg of total RNA was reversed transcribed with high capacity cDNA reverse transcription kits according to the manufacturer's protocol (TaKaRa, Dalian, China). Zero-point-five microliters of cDNA serving as a template in a 20-µL PCR mixture were subsequently used for qRT-PCR analyses with an ABI 7500 Detection System (Applied Biosystems, Foster City, CA, USA) and primers designed using Primer Premier 5.0 software (PREMIER Biosoft, Palo Alto, CA, USA), as listed in Table S2. Eight independent replications were used for each assay. The mRNA abundance of candidate genes was determined using the KAPA SYBR ® FAST qPCR Master Mix (2×) Universal Cocktail (KAPA Biosystems, Boston, MA, USA). A corrected Ct (∆Ct) was calculated by subtracting the β-actin Ct value from that of the target genes for each sample. Relative differences from the control sample (the HRFI groups) were then calculated according to the 2 −∆∆CT method. Differences between the two groups were analyzed using Student's t-test for independent samples in SAS 8.0 for Windows.

Characteristics of Birds in High and Low Phenotypic Groups
A summary of the phenotypic measurements of Beijing-You and Cobb chickens is presented in Table 2. There was no significant difference between HRFI and LRFI sub-groups within breeds in the initial body weights (BW 70 of Beijing-You and BW 28 of Cobb) and the final BW (BW 98 of Beijing-You and BW 42 of Cobb), nor the average daily gain (ADG). Notice, as expected, how different the means are for the Beijing-You and Cobb chickens. There were, however, significant differences (p < 0.01) in daily feed intake (DFI) between the phenotypically-divergent groups in both breeds. Consequently, the differences in mean RFI values and the feed conversion ratio (FCR) between LRFI and HRFI chickens in each breed were highly significant (p < 0.01). The significance of differences between LRFI and HRFI birds was determined using Student's t-test. RFI: residual feed intake; DFI: daily feed intake; ADG: average daily gain over the assessed feeding period; BW: body weight; FCR: feed conversion ratio = DFI/ADG.

Genomic Variants Related to Residual Feed Intake in Beijing-You Chickens
In LRFI and HRFI Beijing-You chickens, 8,505,214 and 8,479,041 SNPs were detected by applying an advanced detection algorithm and filtering criteria (Table S3). There were 26,125 SNPs with F ST values in the top 5% shown in Figure 1. More significant SNPs were enriched on chicken (Gallus gallus) chromosomes (GGA) 1, GGA2, GGA4 and GGA Z. After filtering with allelic frequency divergence >35% between the HRFI and LRFI groups, there were 6288 SNPs identified as candidates for RFI traits (Figure 2A, Table S4). Consistently, over 4000 SNPs were found on four chromosomes GGA 1, GGA2, GGA4 and GGA Z. Among them, 2542 SNPs were located in intergenic regions and 3746 in genes, including upstream and downstream genes. The SNPs with synonymous and missense mutation were only 85 and 27 ( Figure 3A). The 3746 SNPs were then annotated to 1137 genes (Table S4).

Genomic Variants Related to Residual Feed Intake in Beijing-You Chickens
In LRFI and HRFI Beijing-You chickens, 8,505,214 and 8,479,041 SNPs were detected by applying an advanced detection algorithm and filtering criteria (Table S3). There were 26,125 SNPs with FST values in the top 5% shown in Figure 1. More significant SNPs were enriched on chicken (Gallus gallus) chromosomes (GGA) 1, GGA2, GGA4 and GGA Z. After filtering with allelic frequency divergence >35% between the HRFI and LRFI groups, there were 6288 SNPs identified as candidates for RFI traits (Figure 2A, Table S4). Consistently, over 4000 SNPs were found on four chromosomes GGA 1, GGA2, GGA4 and GGA Z. Among them, 2542 SNPs were located in intergenic regions and 3746 in genes, including upstream and downstream genes. The SNPs with synonymous and missense mutation were only 85 and 27 ( Figure 3A). The 3746 SNPs were then annotated to 1137 genes (Table S4).
Functional annotation and pathway analysis were performed using IPA. Specific for Beijing-You, phosphatase and tensin homolog (PTEN) signaling, as well as lipid and carbohydrate metabolism were enriched ( Figure 4, Table 4, Figure S1). Correspondingly, hepatic abundance of NEU3, PHKG1 and PIP5K1B transcripts that participate in lipid and carbohydrate metabolism were upregulated (p < 0.05) in the HRFI compared to the LRFI birds ( Figure 5A).       Functional annotation and pathway analysis were performed using IPA. Specific for Beijing-You, phosphatase and tensin homolog (PTEN) signaling, as well as lipid and carbohydrate metabolism were enriched ( Figure 4, Table 4, Figure S1). Correspondingly, hepatic abundance of NEU3, PHKG1 and PIP5K1B transcripts that participate in lipid and carbohydrate metabolism were upregulated (p < 0.05) in the HRFI compared to the LRFI birds ( Figure 5A).

Genomic Variants Related to Residual Feed Intake in Cobb Chickens
In Cobb chickens, 8,352,008 and 8,372,769 SNPs were detected in the LRFI and HRFI sub-groups (Table S3). There were 15,350 SNPs with FST values in the top 5% shown in Figure 1. More significant SNPs were enriched on GGA1, GGA2, GGA3 and GGA4. With the same data filtering criteria as for Beijing-You chickens, 1001 SNPs were obtained with >600 on GGA 1, GGA2, GGA3, GGA4 and GGA Z ( Figure 2B, Table S5). Of these, 575 SNPs were located within genes, including 10 synonymous mutations and two missense mutations ( Figure 3B). The 575 SNPs were then annotated to 448 genes

Genomic Variants Related to Residual Feed Intake in Cobb Chickens
In Cobb chickens, 8,352,008 and 8,372,769 SNPs were detected in the LRFI and HRFI sub-groups (Table S3). There were 15,350 SNPs with F ST values in the top 5% shown in Figure 1. More significant SNPs were enriched on GGA1, GGA2, GGA3 and GGA4. With the same data filtering criteria as for Beijing-You chickens, 1001 SNPs were obtained with >600 on GGA 1, GGA2, GGA3, GGA4 and GGA Z ( Figure 2B, Table S5). Of these, 575 SNPs were located within genes, including 10 synonymous mutations and two missense mutations ( Figure 3B). The 575 SNPs were then annotated to 448 genes (Table S5). The number of SNPs identified in Beijing-You was more than twice that of Cobb.
Specific for Cobb, function annotation showed that genes were enriched in behavior-related physiological process (Table 3). Feed behavior played a critical role in feed intake, and variation in feed intake was associated with variation in RFI. cAMP responsive element binding protein (CREB) signaling is a key factor in the control of appetite and feed intake by the hypothalamus. CREB signaling in neurons as the canonical pathway was significantly enriched in Cobb chickens ( Figure 6, Table 4). Correspondingly, the expression levels of CAMK4, PLCB4 and FGFR4 in the hypothalamus were higher in HRFI chickens than in LRFI chickens (p < 0.05, Figure 5B).  T cell receptor signaling 6.38 × 10 −3 9.1%, 7/77 Neuropathic pain signaling in dorsal horn neurons 1.02 × 10 −2 8.3%, 7/84 Rac signaling 1.37 × 10 −2 7.9%, 7/89 Actin cytoskeleton signaling 1.42 × 10 −2 8.3%, 7/84 A summary of result from Ingenuity ® IPA software. PTEN: phosphatase and tensin homolog; CNTF: ciliary neurotrophic factor; CREB: cAMP responsive element binding protein; Rac: Ras-related C3 botulinum toxin substrate.

Common Features of Both Breeds
Only 127 enriched genes were common to the two breeds ( Figure 7). There were, however, several physiological processes commonly enriched in the two breeds, such as the development of nervous system, connective tissue and the skeletal and muscular systems (Table 3).

Common Features of Both Breeds
Only 127 enriched genes were common to the two breeds ( Figure 7). There were, however, several physiological processes commonly enriched in the two breeds, such as the development of nervous system, connective tissue and the skeletal and muscular systems (Table 3). In Beijing-You, the enrichment of nNOS signaling in skeletal muscle cells as the significant pathway supported the important role played by organismal development for the RFI trait ( Figure 8, Table 4). Consistently, the abundance of NOS1 transcripts was upregulated (p < 0.01) in the breast muscle of LRFI chickens compared to HRFI chickens ( Figure 5A). In Beijing-You, the enrichment of nNOS signaling in skeletal muscle cells as the significant pathway supported the important role played by organismal development for the RFI trait ( Figure 8, Table 4). Consistently, the abundance of NOS1 transcripts was upregulated (p < 0.01) in the breast muscle of LRFI chickens compared to HRFI chickens ( Figure 5A). Figure 7. Numbers of genes near significant SNPs associated with RFI in Beijing-You and Cobb. There were 1137 genes and 448 genes identified in Beijing-You and Cobb, respectively, of which 127 were commonly found in both breeds.
In Beijing-You, the enrichment of nNOS signaling in skeletal muscle cells as the significant pathway supported the important role played by organismal development for the RFI trait ( Figure 8, Table 4). Consistently, the abundance of NOS1 transcripts was upregulated (p < 0.01) in the breast muscle of LRFI chickens compared to HRFI chickens ( Figure 5A).  Again, in Cobb, the enrichment of Rac and actin cytoskeleton signaling was consistent with the important role played by organismal development for the RFI trait (Table 4, Figures 9 and 10). Transcript abundance of the key genes (CDC42, CSK and PIK3R3) in Rac and actin cytoskeleton signaling were consistently higher (p < 0.01) in breast muscle of LRFI than HRFI chickens ( Figure 5B).
Genes 2018, 9, 28 11 of 18 Again, in Cobb, the enrichment of Rac and actin cytoskeleton signaling was consistent with the important role played by organismal development for the RFI trait (Table 4, Figures 9 and 10). Transcript abundance of the key genes (CDC42, CSK and PIK3R3) in Rac and actin cytoskeleton signaling were consistently higher (p < 0.01) in breast muscle of LRFI than HRFI chickens ( Figure 5B).   Figure 6. The red symbols show the genes that were differentially expressed in breast muscle. Figure 9. Rac signaling pathway in Cobb chickens. Molecular interaction and symbols are the same as the description in Figure 6. The red symbols show the genes that were differentially expressed in breast muscle. Figure 10. Actin cytoskeleton signaling pathway in Cobb chickens. Molecular interaction and symbols are the same as the description in Figure 6. The red symbols show the genes that were differentially expressed in breast muscle.  Figure 6. The red symbols show the genes that were differentially expressed in breast muscle.

SNP Validation
The genome sequencing used DNA from 16 chickens in each of the 12 pools. Nine SNPs were randomly chosen from the candidate SNPs and were subjected to verification by PCR analysis and Sanger sequencing using DNA from individual birds. The results showed (Table S6) that the average similarity was about 75% between the WGS and Sanger sequencing methods, indicating the reliability of the WGS results from the pooling strategy used here.

Association of the SNPs with the Breeding Value of RFI in the Cobb Population of 779
The association analysis to validate the 191 SNPs was done using the validation population of 779 offspring from the Cobb chickens making up the WGS resource (Table S7). There were 46 SNPs that were demonstrated to be significantly associated with RFI (p < 0.05). The details of these significant SNPs, including their positions in the genome, the nearest known genes and the p-values, are given in Table 5. The association between 46 SNPs genotypes and RFI was estimated. Of these, twelve SNPs, the genotypes showed significant differences. The effects of the genotypes and additive and dominance effects are shown in Table 6.

Discussion
Genetic improvement of economically-important traits, such as growth rate and feed efficiency, has occurred for about 100 years in the Cobb chickens. There has been no systematic genetic selection in the conservation population of the local Beijing-You chickens. Consistent with the current findings, more genomic variants exist in the Beijing-You than Cobb chickens, because the high intensity of genetic selection in the latter would be expected to result in genomic regions with increased homozygosity. There were numerous SNPs associated with RFI in the Beijing-You and a larger proportional difference in FCR (HRFI/LRFI = 1.2), suggesting a high potential for genetically-increasing feed efficiency in Beijing-You chickens. Obviously, FCR is already much lower in Cobb chickens, and there was less difference (HRFI/LRFI = 1.08) between the phenotypic extremes.
The two breeds examined here showed some unique characteristics, and only 127 genes associated with RFI were identified in both breeds. These represented a small proportion of the total (127/1137) for Beijing-You, but much higher (127/448) for Cobb chickens. PTEN signaling, as well as lipid and carbohydrate metabolism were only enriched in Beijing-You chickens. The PTEN/PI3K/AKT signaling pathway has been shown to regulate cell metabolism, growth and survival [16,17]. PTEN/PI3K/AKT signaling can mediate the insulin-stimulated activation of many major pathways of hepatic lipid metabolism, including lipogenesis, fatty acid β-oxidation, fatty acid esterification to triglycerides and very low density lipoprotein-triglyceride (VLDL-TG) secretion [18,19]. PTEN, a negative regulator, dephosphorylates PIP3 to inhibit the activity of PI3K/AKT signaling [20]. When the PTEN was overexpressed, insulin signaling effects were inhibited, resulting in decreased RAC-alpha serine/threonine-protein kinase (AKT) activity and solute carrier family 2, facilitated glucose transporter member 4 (GLUT4) translocation to the cell membrane, which consequently contributed to insulin resistance and progression of non-alcoholic fatty liver disease [21][22][23]. In addition to controlling the activity of signaling intermediates at the cell membrane, PI3Ks can drive gene expression programs in the nucleus. The forkhead box protein O1(FOXO) transcription factors are coupled to PI3K activity through AKT-dependent phosphorylation [24], which are considered effectors of the pathway regulating hepatic glucose production. Matsumoto et al. [25] showed that delivery of FOXO1 to mouse liver resulted in steatosis arising from increased triglyceride accumulation and decreased fatty acid oxidation. In the present study, the genes (PHKG1, PIP5K1B and NEU3) that participate in carbohydrate and lipid metabolism [26][27][28] were upregulated in the livers of HRFI Beijing-You chickens. Previous studies [5,29] have demonstrated that HRFI animals required more energy to maintain biological activities, such as protein turnover, stress response and muscle activity. Beijing-You chickens are characterized by high fat deposition in muscle and abdominal adipose [30]. Therefore, lipid and carbohydrate metabolism may play crucial roles in RFI through the PTEN/PI3K/AKT signaling pathway in Beijing-You chickens.
For Cobb chickens, the results showed that feed behavior and the related pathway (CREB signaling in neurons) were enriched. Feed behavior plays a key role in feed intake, and variation in feed intake is associated with variation in RFI. The CREB transcription factor has been widely investigated as a key factor in the control of appetite and FI by the hypothalamus [31]. A previous report showed that an increase and phosphorylation of CREB in the hypothalamus, as a result of injection of NPY into the rat hypothalamus or fasting for 48 h, suggested that CREB phosphorylation might be involved in the regulation of feeding behavior [32]. Calcium/calmodulin-dependent protein kinase (CAMK) plays a critical role in activation of CREB via Ser-133 phosphorylation [33]. PLCB4 and CAMK4 mRNA levels were higher in the hypothalamus of HRFI than LRFI chickens. The upregulation of CAMK4 mRNA levels may increase FI through induction of CREB phosphorylation and activation of downstream signaling pathways, and upregulation of PLCB4 expression may affect activation of CAMKs through changes in Ca 2+ concentrations. Therefore, the CREB signaling pathway may contribute to RFI in Cobb chickens by affecting feed behavior. Compared to the Beijing-You, the Cobb is a specialized broiler breed as a result of long-term selection for increased feed efficiency and growth rate, which may have resulted in part from changes in feed behavior and its related genetic basis.
Organismal development as the common physiological process was enriched in both breeds. The concept of RFI was initially developed to determine differences among animals in feed conversion into body tissue. A previous study in cattle showed a variation in RFI of 37%, which was explained by tissue metabolism, protein turnover and stress response [34]. Global RNA expression in breast muscle obtained from a male broiler line phenotyped for high or low feed efficiency showed that significant differences in gene expression profiles were associated with muscle fiber development, muscle function and cytoskeletal organization [35]. In the present study, skeletal and muscular system development and function were enriched in both Beijing-You and Cobb breeds. Phenotypic measurement showed that the percentages of breast and thigh muscle were higher in LRFI than HRFI Beijing-You chickens ( Table S8), suggesting that the development and function of the skeletal and muscular systems might control efficient conversion of dietary nutrients to lean tissue and growth. The canonical pathways that participated in the organismal development were also enriched in both Beijing-You and Cobb breeds. nNOS plays an important role in the regulation of many muscle functions, such as blood flow, contraction and metabolism [36]. For example, nitric oxide (NO), which is physiologically synthesized in a tightly-controlled manner by NO synthases (NOSs), directly inhibits mitochondrial respiration by interfering with electron transport chain complexes, which consequently affects skeletal muscle metabolism and stimulates glucose transport by activating upstream signaling that leads to increased production of GLUT4 [37]. Rac and actin cytoskeleton signaling has been identified in various morphogenetic events associated with tissue organization during organ development. For instance, the conserved role of paxillin in promoting signals from RAC to the actin cytoskeleton is crucial for normal tissue morphogenesis during wing and leg development of Drosophila [38]. Ras-related C3 botulinum toxin substrate 1 (Rac1) and Cell division control protein 42 homolog (CDC42) have also been implicated in neural development and myogenesis [39]. CSK-1 was also essential in the organization of pharyngeal muscle filaments [40]. In the present study, expression in muscle of CDC42, CSK and PIK3R was significantly upregulated in LRFI chickens, as compared to HRFI chickens. Therefore, all of these results showed that although there were small breed differences in the participating pathways, organismal development as the common physiological process contributed to the RFI in the two breeds.
In order to further validate the reliability of the WGS results, 191 SNPs were selected for genotyping and association analysis. The results showed that there were 46 SNPs (48 genes) significantly associated with the RFI, which demonstrated that the WGS result was reliable, and these SNPs could prove to be interesting markers for RFI.

Conclusions
The present study provides novel insights into genomic variants related to RFI in two very distinct breeds of chicken. A total of 3746 reliable SNPs included in 1137 genes in Beijing-You chickens and 448 genes with 575 SNPs in Cobb chickens were identified as significant variant markers associated with RFI. Forty-six SNPs were significantly associated with the RFI in an independent population of 779 Cobb chickens, suggesting that the method of screening the associated SNP with WGS strategy was reasonable. Lipid and carbohydrate metabolism and the PTEN signaling pathway were enriched in the Beijing-You, while feed behavior and the CREB signaling pathway were enriched in Cobb chickens. Organismal development as the common physiological process was enriched in the two breeds. NOS1, PHKG1, NEU3 and PIP5K1B were differentially expressed in Beijing-You, while CDC42, CSK, PIK3R3, CAMK4 and PLCB4 were differentially expressed in Cobb, suggesting that these might be key genes that contribute to RFI. This present study, with chickens of very distinct backgrounds and characteristics, offers important knowledge of biomarkers, candidate genes and potential biological pathways underlying RFI.