Genetic Architecture of Feeding Behavior and Feed Efficiency in a Duroc Pig Population

Increasing feed efficiency is a major goal of breeders as it can reduce production cost and energy consumption. However, the genetic architecture of feeding behavior and feed efficiency traits remains elusive. To investigate the genetic architecture of feed efficiency in pigs, three feeding behavior traits (daily feed intake, number of daily visits to feeder, and duration of each visit) and two feed efficiency traits (feed conversion ratio and residual feed intake) were considered. We performed genome-wide association studies (GWASs) of the five traits using a population of 1,008 Duroc pigs genotyped with an Illumina Porcine SNP50K BeadChip. A total of 9 genome-wide (P < 1.54E-06) and 35 suggestive (P < 3.08E-05) single nucleotide polymorphisms (SNPs) were detected. Two pleiotropic quantitative trait loci (QTLs) on SSC 1 and SSC 7 were found to affect more than one trait. Markers WU_10.2_7_18377044 and DRGA0001676 are two key SNPs for these two pleiotropic QTLs. Marker WU_10.2_7_18377044 on SSC 7 contributed 2.16 and 2.37% of the observed phenotypic variance for DFI and RFI, respectively. The other SNP DRGA0001676 on SSC 1 explained 3.22 and 5.46% of the observed phenotypic variance for FCR and RFI, respectively. Finally, functions of candidate genes and gene set enrichment analysis indicate that most of the significant pathways are associated with hormonal and digestive gland secretion during feeding. This study advances our understanding of the genetic mechanisms of feeding behavior and feed efficiency traits and provide an opportunity for increasing feeding efficiency using marker-assisted selection or genomic selection in pigs.


INTRODUCTION
Pork is an important meat source for humans, accounting for nearly 40% of all meat consumed by the world population (Wang et al., 2015). The share of feed cost, which is the highest of the total production cost, remains high, ranging from 50 to 85% (Young et al., 2011). The key to reducing feed cost is to increase feed efficiency. Increasing feed efficiency not only reduces feed consumption while decreasing farming cost and energy use, but also lowers manure production and the total amount of potential greenhouse gas emission (O'Shea et al., 2012;Bartos et al., 2016).
In animal breeding programs, feed efficiency traits are difficult to improve by direct selection because feed efficiency cannot be measured directly (Hoque and Suzuki, 2008). Feed conversion ratio (FCR) and residual feed intake (RFI) are two traits that have been used to evaluate feed efficiency (Lu et al., 2017). FCR (the ratio of feed intake to output) is widely used to estimate feed efficiency in pig breeding because of its simplicity in calculation and its correlation with growth rate and body weight (BW) (Fan et al., 2010;Do et al., 2013a). However, previous genetic studies also show that FCR is not always effective (Webb and King, 1983). For example, pigs with low feed intakes and undesirably low gains may also have high FCRs (Webb and King, 1983;Iwaisaki, 1989;Smith et al., 1991). RFI is defined as the difference of feed intake between the actual feed eaten and the expected feed intake required for production and maintenance (Koch et al., 1963). RFI appears to be a better indicator of feed efficiency as compared to FCR. However, computation of RFI varies widely, depending on the predicted feed requirement for production and maintenance (Hoque et al., 2009;Do et al., 2013aDo et al., , 2014, which adds to the difficulty of comparing different studies. Feeding behavior is one of the most important factors affecting feed efficiency. Several studies have shown a strong genetic and phenotypic correlation between feeding behavior and feed efficiency traits in swine (Hoque et al., 2009). For instance, daily feed intake (DFI) has a positive genetic and phenotypic correlation with FCR (0.65 and 0.67, respectively) and RFI (0.95 and 0.90, respectively) (Do et al., 2013a). Therefore, to understand the molecular mechanism and genetic basis of feed efficiency, it is likely helpful to consider both feeding behavior and feed efficiency traits.
Among a total of 26,076 quantitative trait loci (QTL) associated with 647 different traits reported in pigs (PigQTLdb 1 ), 639 QTLs have been identified for feeding behavior and feed efficiency traits. Most of these QTL were identified using linkage mapping. Because of the large intervals of QTLs, directly using them for genetic improvement remains difficult (Tabor et al., 2002). In recent years, with the advent of dense marker panels, association mapping has become a powerful strategy for the detection of genetic variants associated with complex traits. It has been widely used in humans (McCarthy et al., 2008) and domestic animals (Andersson, 2009;Ma et al., 2014;Guo et al., 2016).
Duroc pig population is widely used as the terminal male parent of the DLY (Duroc × Landrace × Yorkshire) commercial pigs thanks to its excellent performance on growth traits. In previous studies, several single nucleotide polymorphisms (SNPs) that were significantly associated with FCR were detected on SSC 12 in a Canadian Duroc population by genome-wide association study (GWAS) . However, due to the limited sample size (<400), we detected only a few significant SNPs for feeding behavior and feed efficiency traits. In particular, the cost to measure feed efficiency traits has historically been the primary limitation to population-wide selection to improve feed efficiency in pigs (Bai et al., 2017). Here, we perform a GWAS in a larger American Duroc population to identify genetic variants associated with feeding behavior and feed efficiency traits.

Ethics Statement
The experimental procedures used in this study met the guidelines of the Animal Care and Use Committee of the South China Agricultural University (SCAU) (Guangzhou, China). The Animal Care and Use Committee of the SCAU approved all the animal experiments described in this study.

Animals and Phenotype
During the period of 2013-2016, phenotypic data were collected for Duroc pigs (n = 1,008) from the Guangdong Wen's Foodstuffs Group, Co., Ltd. (Guangdong, China) using the Osborne FIRE Pig Performance Testing System (Osborne, KS, United States) as previously described . Briefly, pigs were subjected to uniform feeding conditions for measurement of traits during the fattening period (approximately 11 weeks) from 30 to 100 kg live weight. Each animal was labeled a unique electric identification tag on its ear that could be captured by the automatic feeder. The time, duration, feed consumption, and BW of each individual were recorded at every visit to the feeder. Back fat (BF) of each animal was evaluated with a PIGLOG 105B ultrasound machine (SFK Technology, Søborg, Denmark) at the end of the test. The following feeding behavior and feed efficiency traits were defined and recorded for each pig: average DFI (kg/d), total daily time spent in feeder (TPD, min), number of daily visits to feeder (NVD), and FCR (Do et al., 2013a(Do et al., ,b, 2014Ding et al., 2017). RFI was computed using methods similar to those used by Cai et al. (2008). In the model, predicted DFI was estimated using linear regression of DFI on metabolic BW at mid-test (MWT), average daily gain from 30 to 100 kg (ADG), and BF. MWT was equal to [(BW at on-test + BW at off-test)/2] 0.75 . In summary, 1,008 pigs had phenotypic data for feeding behavior traits (DFI, TPD, and NVD), 981 for FCR and 971 for RFI.

Genotyping and Quality Control
Genotyping was performed as described by Ding et al. (2017). Briefly, genomic DNA was extracted from ear tissue samples using the phenol-chloroform method. Genotyping was performed using the Geneseek Porcine 50K SNP Chip (Neogen, Lincoln, NE, United States), which contains 50,703 SNPs across autosomes and sex chromosomes. Quality control of the SNP data was performed using PLINK software (Purcell et al., 2007). Animals with call rates > 0.95, and SNPs with call rates > 0.99, minor allele frequency > 0.01, and P-value > 10 −6 for the Hardy-Weinberg equilibrium test were included. Only autosomal SNPs were considered for subsequent analyses.

GWAS and Genetic Analyses
GEMMA was used to perform GWAS with a univariate linear mixed model Stephens, 2012, 2014). Prior to GWAS, GEMMA was used to estimate the n × n standardized relatedness matrix (K) between the individuals. The following statistical model in matrix form was used: y = Wα + xβ + u + ε , where y is the vector of phenotypic values for all pigs; W is the incidence matrices of covariates (fixed effects) including sex, pig pen, and year-season effects; α is the vector of corresponding coefficients including the intercept; x is the vector of marker genotypes, and β is the corresponding effect size of the marker; u is the vector of random effects, with u ∼ MVN n (0, λ τ −1 K); ε is the vector of random residuals with ε ∼ MVN n (0, τ −1 I n ); τ −1 is the variance of the residual errors; λ is the ratio between the two variance components; K is a known n × n relatedness matrix; and I is an n × n identity matrix. MVN n denotes the n-dimensional multivariate normal distribution.
Genome-wide significance was determined using the Bonferroni method by dividing the desired type I error level by the number of SNPs tested. To include additional candidate genes and enable gene set enrichment analysis, we also set a more lenient threshold by multiplying the Bonferroni threshold by a constant of 20 (Yang et al., 2005).
The GCTA tool was used to compute the genomic heritability by dividing the estimated genetic variance by the total variance measured, and phenotypic variances contributed by significant SNPs for each trait (Yang et al., 2010(Yang et al., , 2011. Genetic correlation was estimated using GCTA in the bivariate mode. A number of SNPs that were significantly associated with the target trait by GWAS were detected based on their strong linkage with highly causal mutants. To demarcate the independence of all significant signals in a putative region, conditional and LD analyses were performed in univariate linear mixed models by fitting the genotypes of peak SNPs as covariates (Yang et al., 2012). Moreover, to further detect candidate regions associated with feeding behavior and feed efficiency traits, PLINK (Purcell et al., 2007) and Haploview (Barrett et al., 2005) were utilized for haplotype block analysis. LD blocks were defined using the solid spin algorithm by the criteria of Gabriel et al. (2002).

Annotation of SNPs
All SNP location on the Sus scrofa 10.2 genome version were downloaded from Ensembl. The Ensembl annotation of the S. scrofa 10.2 genome version was employed to find genes that were nearest to the significant SNPs 1 . To annotate significant SNP located in previously mapped QTLs in pigs, all QTL data in pigs were downloaded from http://www.animalgenome. org/cgi-bin/QTLdb/SS/download?file=gbpSS_10.2 (accessed on December 18, 2017) (Hu et al., 2013). For functional annotation, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology analyses were used for the identification of related pathways. KEGG and GO analyses were performed on KOBAS 3.0 (Xie et al., 2011). Fisher's exact test was used to assess the significance of the enriched terms, and P < 0.05 was selected to explore the genes involved in biological processes (Xing et al., 2016;Wang et al., 2017).

Quantitative Genetics of Feeding Behavior and Feed Efficiency Traits in Pigs
We considered five feeding behavior and feed efficiency traits, including DFI, TPD, NVD, FCR, and RFI. Summary statistics for these traits, and their heritabilities are presented in Table 1.
All phenotypic data conformed to the Gaussian distribution based on the Shapiro test before GWAS analysis (Theune, 1973). There existed substantial phenotypic variation, with coefficient of variation (CV) ranging from 10 to 27% for the five traits (Table 1). We partitioned the phenotypic variance into genetic and environmental components using genetic relationship matrix computed from genotypes. The genomic heritabilities of the traits are moderate, ranging from 0.28 (DFI) to 0.38 (NVD) ( Table 2). Bivariate analysis indicated that the traits are positively correlated with each other, both phenotypically and genetically, except between DFI and NVD. In particular, DFI are strongly positively correlated with FCR and RFI with a genetic correlation of 0.75 and 0.98, respectively ( Table 2). Corresponding SD behind the mean heritabilities on diagonal (bold), phenotypic correlations above diagonal, genetic correlations below diagonal.

Genome-Wide Association Studies
To understand the genetic architecture of the feeding behavior and feed efficiency traits, we performed GWAS for each trait. After a series of filtering steps for quality control, 32,446 SNPs and 1,008 pigs were available for subsequent GWAS. The number of SNPs on each chromosome and the average distances between pairs of SNPs after QC are provided in Supplementary Table S1. The average physical distance between two neighboring SNPs on the same chromosome was approximately 75.3 kb and ranged from 64.5 (SSC11) to 112.3 kb (SSC1). Single marker tests using mixed model were performed to identify genetic markers associated with the traits. At a stringent genome-wide Bonferroni threshold P < 1.54E-06 (0.05/32,446), 9 SNPs on SSC1 were associated with RFI. At a more lenient threshold (P < 3.08E-05) for suggestive associations, 8 SNPs were associated with DFI, 15 with FCR, 5 with NVD, 2 with TPD, and 14 with RFI (Figure 1 and Table 2). QQ plots of P-values and the computed genomic inflation factors (λ) indicated no evidence of population stratification (Pearson and Manolio, 2008;Utsunomiya et al., 2013) (Supplementary Figure S1).
Multiple SNPs in close proximity were found to be associated with the same traits, possibly due to their linkage and/or linkage disequilibrium. Indeed, LD block analysis showed that the multiple significant SNPs on SSC1 associated with RFI were located within a haplotype block that spanned 2,169 kb (Figure 4). To test whether LD caused the associations, we performed conditional analyses with DFI and RFI in which the lead SNP WU_10.2_7_18377044 on SSC7 was fitted in the model as a covariate and the conditional P-values for other SNPs in the vicinity were obtained. While many SNPs in high LD with the lead SNP were significant in the GWAS for DFI (Figure 2A), their significance almost diminished completely after the lead SNP was included as a fixed effect in the model ( Figure 2B). The same pattern was also observed for the same lead SNP for RFI (Figures 2C,D). This is not surprising because the two DFI and RFI were almost perfectly (r = 0.98) correlated genetically (Table 2).
Similarly, the lead SNP DRGA0001676 significantly associated with RFI and FCR in SSC1 explained the association between multiple SNPs with the traits in the vicinity of DRGA0001676 (Figure 3).

Comparison With Previously Mapped Pig QTLs
To evaluate whether SNPs associated with the feeding behavior and feed efficiency trait in this study replicate any previously known QTLs, we search the pigQTLdb based on SNP and QTL locations. A total of 13 SNPs associated with FCR and/or RFI were identified within the genomic regions where QTLs for DFI have been previously mapped in pigs (Table 3). Eight SNPs associated with DFI and RFI were located on previously reported QTL regions for time spent drinking. One SNP on SSC13 associated with TPD was located on previously reported QTL regions for FCR in pigs. Moreover, 10 SNPs were located in the genomic regions where QTLs were previously detected by GWAS or linkage mapping for average daily gain and growth-related traits in pigs. Given the high genetic correlations ( Table 2), overlaps were considered with QTLs for correlated traits as evidence for replication.

Candidate Genes and Functional Analysis
A total of 16 functional genes that were within or near the identified tag SNPs were detected based on annotations of the Sus scrofa 10.2 genome assembly ( Table 2). Many candidate genes appear to have biochemical and physiological roles that were are relevant to feeding behavior and feed efficiency traits, including neurensin 1 (NRSN1) and doublecortin domaincontaining 2 (DCDC2) for DFI; ADAM metallopeptidase domain 12 (ADAM12) for TPD; phospholipase C beta 1 (PLCB1), GNAS complex locus (GNAS), and ephrin B2 (EFNB2) for FCR; prolactin (PRL) for both FCR and RFI; leucine-rich repeat and fibronectin type III domain-containing 5 (LRFN5), and multiple EGF-like domains 11 (MEGF11) for RFI and FCR.
To uncover genes and pathways involved in the biology of feed efficiency and feeding behavior traits, we performed gene set enrichment analysis of the 16 candidate genes for all five traits. Interestingly, at an FDR = 0.05, several KEGG pathways and GO terms are significantly enriched for the candidate genes, including pathways related to hormone metabolic process, secretion of digestive enzymes, among others (Supplementary Table S2).

DISCUSSION
In this study, we performed a GWAS for five feeding behavior and feed efficiency traits. We identified a number of genetic markers and genes associated with the traits. Gene set enrichment analysis revealed pathways that appear to be consistent with the underlying biology of the traits. This represents a first step toward understanding the genetic basis of feeding behaviors and feed efficiency and developing informative genetic markers for efficient breeding programs.
Despite similar design and analytical methods, none of the QTLs identified in this study replicated QTLs in a previous study . There could be a number of reasons. First, this study has a substantially larger sample size (1,008 vs. 338). Second, while both studies used Duroc pigs, the genetic backgrounds have subtle differences. For example, Duroc boars in the previous study  was of Canadian origin while pigs in the present study were of American origin. It is well-established that genetic backgrounds can have substantial influence on single marker associations. Moreover, this study detected QTLs that overlapped previously identified QTLs (Table 4). Indeed, over 87% of SNPs (28/32) significantly associated with feeding behavior and feed efficiency traits were located within genomic regions where QTLs for feeding behavior and feed efficiency were previously reported in pigs (pigQTLdb 2 ). Derivative traits are correlated with each other. For example, RFI and DFI has a highly correlated at 0.81 phenotypically and 0.98 genetically. Indeed, GWAS signals overlapped to some extent for these two traits (Figures 1A,B). Signals are detected for both traits on SSC7. And while the strengths of associations differ and may not reach thresholds for both traits, SNPs associated with RFI on SSC1 also tended to have an effect for DFI (Figures 1A,B). However, because of the less than unity correlations, in particular the correlation at the phenotypic level which were used for GWAS mapping, full overlap is not expected.
Because of the low to moderate heritabilities of these traits, informative markers are particularly useful in selection programs. Our analyses indicated that the markers WU_10.2_7_18377044 and DRGA0001676 were two important polymorphisms Feeding is one of the most conserved activities of animals, and regulating feed intake is a fundamental process for animal survival (Bader et al., 2007). Physiological modulation is often accomplished through the regulation of the nervous system (Destexhe and Marder, 2004). In mammals, the hypothalamus controls appetite and satiety by integrating neuropeptides or hormonal signals (Berthoud, 2002;Leibowitz and Wortley, 2004;Fang et al., 2011). Our GWAS implicated several genes that have roles in neural development. The gene NRSN1 associated with DFI plays an important role in neural organelle transport and in the transduction of nerve signals or in nerve growth (Yu and Kaang, 2016;Lencer et al., 2017). Three significant SNPs located within the DCDC2 gene are associated with DFI. DCDC2 is involved in the conduction of neural signals and the development of neurons (Meng et al., 2005). Takasu et al. (2002) reported that EFNB2 protein, which is associated with FCR, is localized at excitatory synapses, regulating the development and remodeling of neural signals. MEGF11, another gene associated with FCR, is highly expressed in the central nervous system, whose homology MEGF10 in Drosophila plays an important role in maintaining the normal functioning of the brain (Draper et al., 2014). Finally, LRFN5, also a FCR associated gene, is a key neurodevelopmental gene that is associated with developmental delay (Mikhail et al., 2011).
Daily feed intake is a major concern of breeders because it is highly correlated with growth rate and feed efficiency (Do et al., 2013a). The most significant locus associated with DFI and RFI on SSC 7, WU_10.2_7_18377044, is close to the PRL gene, which encodes the anterior pituitary hormone prolactin. The secreted hormone is a growth regulator for many tissues (McAtee and Trenkle, 1971;Kennaway et al., 1982). Yayou et al. (2011) reported that administration of prolactin-releasing peptide in the central nervous system will result in reduced feed intake and abnormal feeding behavior in steers. In goats, prolactin has also been found to affect feed intake (Sergent et al., 1988). FIGURE 4 | The LD block in the significantly associated region on SSC1. LD blocks are marked with triangles. Values in boxes are LD (r 2 ) between SNP pairs and the boxes are colored according to the standard Haploview color scheme: LOD > 2 and D = 1, red; LOD > 2 and D < 1, shades of pink/red; LOD < 2 and D = 1, blue; LOD < 2 and D < 1, white (LOD is the log of the likelihood odds ratio, a measure of confidence in the value of D ). Annotated genes in the chromosomal region retrieved from the Ensemble genome browser (www.ensembl.org/Sus_scrofa/Info/Index).
Growth rate and feed intake were major influencing factors of FCR. Among the candidate genes identified by GWAS, PLCB1 is a protein-coding gene near the SNP ALGA0093563. Meng et al. (2017) reported that this gene is significantly associated with ADG in a Yorkshire pig and may be involved in pig growth and development. The GNAS gene is near SNP WU_10.2_17_66292358. The imprinted GNAS is involved in obesity, energy metabolism, feeding behavior, and viability. Eaton et al. (2012) reported that mutations in this gene can lead to preweaning growth retardation and incomplete catch-up growth in mice.
Feeding time is an important feeding behavior trait that directly affects the growth rate and efficiency. ADAM metallopeptidase domain 12 (ADAM12) is located proximal to the SNP WU_10.2_14_14709705, which is associated with TPD. ADAM12 had been identified as a protease to insulin-like growth factor (IGF)-binding proteins. This gene may have a regulatory function in controlling the amount of free bioactive IGF (Cowans and Spencer, 2007). In chickens, insulin has been consistently shown to affect feed intake and feed behavior (Shiraishi et al., 2008(Shiraishi et al., , 2011. Functional annotation revealed a number of pathways and biological processes that are significantly overrepresented among the 16 positional candidate genes for feeding behavior and feed efficiency traits (Supplementary Table S2). Most of the significantly enriched pathways are associated with hormonal and digestive gland secretion during feeding such as thyroid hormone and gastric acid secretion. Given the principal roles of hormone and digestive gland secretion in feeding and metabolism, their involvement in the feeding behavior and feed efficiency traits are conceivable.
Several GWAS for feeding behavior and feed efficiency traits in Duroc populations with different genetic backgrounds have been reported (Do et al., 2013a,b;Ding et al., 2017;Quan et al., 2017). We previously found, by GO enrichment analysis, that most of the potential candidate genes were involved in the development of the hypothalamus . This agrees with our findings that feeding behavior is mainly controlled by the central nervous system. Results in the present study further extend to suggest that the hormonal and digestive glands are also involved. However, functional studies remain to be performed to delineate the mechanisms of how the constellations of genes implicated by the GWAS affect feeding.

AUTHOR CONTRIBUTIONS
ZW, JY, and WH: conceived and designed the experiments. RD, MY, XW, JQ, ZZ, SZ, and SL: performed the experiments. RD, JY, and WH: analyzed the data and wrote manuscript. MY, ZX, EZ, GC, and DL: collected the samples and recorded the phenotypes. ZW: contributed materials. All authors reviewed and approved the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene. 2018.00220/full#supplementary-material FIGURE S1 | Quantile-quantile (Q-Q) plots of genome-wide association studies for feeding behavior and feed efficiency in male Duroc pigs. Q-Q plots show the observed versus expected negative log10 P-values. On the vertical axis, Q-Q plot for (A) residual feed intake (RFI), (B) total daily feed intake (DFI), (C) feed conversion ratio (FCR), (D) total daily time spent in feeder (TBD), and (E) number of visits to feeder (NVD), respectively.