Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 https://doi.org/10.1186/s12711-019-0490-6 Ge n e t i c s Se l e c t i o n Ev o l u t i o n RESEARCH ARTICLE Open Access Integrating genome‑wide co‑association and gene expression to identify putative regulators and predictors of feed efficiency in pigs Yuliaxis Ramayo‑Caldas1, Emilio Mármol‑Sánchez2, Maria Ballester1, Juan Pablo Sánchez1, Rayner González‑Prendes2, Marcel Amills2,3 and Raquel Quintanilla1* Abstract Background: Feed efficiency (FE) has a major impact on the economic sustainability of pig production. We used a systems-based approach that integrates single nucleotide polymorphism (SNP) co-association and gene-expression data to identify candidate genes, biological pathways, and potential predictors of FE in a Duroc pig population. Results: We applied an association weight matrix (AWM) approach to analyse the results from genome-wide associa‑ tion studies (GWAS) for nine FE associated and production traits using 31K SNPs by defining residual feed intake (RFI) as the target phenotype. The resulting co-association network was formed by 829 SNPs. Additive effects of this SNP panel explained 61% of the phenotypic variance of RFI, and the resulting phenotype prediction accuracy estimated by cross-validation was 0.65 (vs. 0.20 using pedigree-based best linear unbiased prediction and 0.12 using the 31K SNPs). Sixty-eight transcription factor (TF) genes were identified in the co-association network; based on the lossless approach, the putative main regulators were COPS5, GTF2H5, RUNX1, HDAC4, ESR1, USP16, SMARCA2 and GTF2F2. Fur‑ thermore, gene expression data of the gluteus medius muscle was explored through differential expression and mul‑ tivariate analyses. A list of candidate genes showing functional and/or structural associations with FE was elaborated based on results from both AWM and gene expression analyses, and included the aforementioned TF genes and other ones that have key roles in metabolism, e.g. ESRRG​, RXRG, PPARGC1A, TCF7L2, LHX4, MAML2, NFATC3, NFKBIZ, TCEA1, CDCA7L, LZTFL1 or CBFB. The most enriched biological pathways in this list were associated with behaviour, immunity, nervous system, and neurotransmitters, including melatonin, glutamate receptor, and gustation pathways. Finally, an expression GWAS allowed identifying 269 SNPs associated with the candidate genes’ expression (eSNPs). Addition of these eSNPs to the AWM panel of 829 SNPs did not improve the accuracy of genomic predictions. Conclusions: Candidate genes that have a direct or indirect effect on FE-related traits belong to various biological processes that are mainly related to immunity, behaviour, energy metabolism, and the nervous system. The pituitary gland, hypothalamus and thyroid axis, and estrogen signalling play fundamental roles in the regulation of FE in pigs. The 829 selected SNPs explained 61% of the phenotypic variance of RFI, which constitutes a promising perspective for applying genetic selection on FE relying on molecular-based prediction. *Correspondence: raquel.quintanilla@irta.cat 1 Animal Breeding and Genetics Program, Institute for Research and Technology in Food and Agriculture (IRTA), Torre Marimon, 08140 Caldes de Montbui, Spain Full list of author information is available at the end of the article © The Author(s) 2019. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat​iveco​mmons​.org/licen​ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/ publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 2 of 17 Background The objective of our study was to use a systems-based Improving feed efficiency (FE) has become a relevant approach that integrates information from several FE- but challenging focus for pig breeding selection schemes related phenotypes, SNP co-associated networks, and due to its strong influence on the economic sustainabil- gene-expression data to disentangle the molecular mech- ity and environmental impact of pig production. FE is a anisms that underlie FE in pigs. Our aim was also to complex phenotype that depends on genetics [1–3], on identify candidate genes and their potential regulators in the health and physiological status of the animals [4], on order to develop a panel of markers that can be used as environmental factors [5, 6], and on the gut microbial predictors of an individual’s genetic potential for FE. composition [7–9]. An additional complexity for genetic improvement of FE is the definition of adequate selec- Methods tion criteria relative to feed intake and other produc- In this study, we integrated several genome and gene tion traits, such as growth, which in turn show a mutual expression analyses that aimed at identifying candidate dependency. Besides the conventional feed conversion genes, putative regulators, and predictors of FE in pigs. ratio (FCR), the most widely used measure of FE dur- The different approaches used are described below and ing the last decade is residual feed intake (RFI), i.e. the outlined in Additional file 1: Figure S1. deviation of the animal’s feed intake from the amount of feed predicted to be required for maintenance, growth, and back fat deposition [10]. Estimates of the heritabil- Animals and phenotypes ity of RFI range from 0.14 to 0.53 [1–3], whereas the A population of 350 Duroc barrows from five paternal reported accuracies of RFI genomic prediction range half-sib families was used in this experiment. Pigs were from 0.40 to 0.53 [11, 12]. Selection for RFI has proven raised under intensive standard conditions at the Institut to be a successful strategy for improving FE in pigs [1, de Recerca i Tecnologia Agroalimentàries (IRTA) Experi- 2, 13] but requires recording of individual feed intake, mental Pig Farm in Monells (Girona). They were distrib- body weight gain, and back fat, which is expensive and uted across four fattening batches in a partially balanced time-consuming. Thus, identifying candidate genes but connected design: one batch included the offspring of and potential regulators of FE that are predictive of the all five sires, while the remaining batches contained the animal’s genetic potential for this phenotype is of para- offspring of four of the five sires. All animals were sub- mount interest. jected to the same management procedures, with ad libi- Several studies at both the genomic and transcriptomic tum access to feed, with two standard diets with energy levels have been performed to identify candidate sin- densities of 10.27 MJ/kg until the animals reached 90 kg gle nucleotide polymorphisms (SNPs) associated to FE of live weight (▲150 days of age), and of 9.94 MJ/kg dur- and to unravel the genetic architecture of this complex ing the last ▼40 days before slaughter. A detailed descrip- trait in pigs [14–16], including genome-wide association tion of the experimental population and management studies (GWAS) [14, 17] and whole-genome expression conditions is in Gallardo et al. [29] and Quintanilla et al. analyses [15, 16, 18]. Among the tissues that are relevant [30]. Animal care and experimental procedures were per- for FE, several studies have focused on the muscle tran- formed by following national and institutional guidelines scriptome [19–21]. Muscle plays a central role in main- for the Good Experimental Practices and were approved taining overall energy balance by controlling the storage by the IRTA Ethical Committee. Pigs were weighed indi- of lipids and carbohydrates [22, 23]. However, the results vidually at ~ 65 days of age and every 3 weeks during the from different studies are not always consistent and most fattening period, plus on the day of slaughtering (~ 190 studies ignore gene-by-gene interactions and fail to inte- days of age). Backfat thickness (BF) was also measured grate structural and functional genomics data. Holistic every 3 weeks using PIGLOG 105 ultrasound equipment. ers located in each pen ­(IVO®-feeding station, Insentec, approaches that combine multiple sources of informa- Individual feed intake was recorded by electronic feed- tion increase the power to identify candidate genes [24, 25] and provide a more complete picture of the biological Marknesse, The Netherlands). The average daily feed processes under investigation. To date, a limited number intake (ADFI) of each individual during the trial was of studies targeting FE in pigs and other livestock species computed. Two measures of individual FE during fatten- have implemented integrative approaches, such as system ing were computed: FCR, measured as the simple ratio genetics- or gene network-based methods [26–28]. of ADFI and average daily gain (ADG) (kg/kg), and RFI, computed as the residual of the following model: Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 3 of 17 where ADFIij is the average daily feed intake of individual (target) phenotype in the AWM procedure, whereas the i (in batch j ) during the whole fattening period (from ~ 70 other traits were selected based on their association with to ~ 190 days of age); bj is the effect of batch level j(four FE, but also on their relevance in pig production. In a first batches); Ai is the age of individual i at the midpoint of step, a GWAS was performed for each of the nine afore- the analysed period (135 days on average but ranging mentioned traits by using the genome-wide complex trait from 121 to 148 days), and α is the corresponding regres- analysis (GCTA) software [36]. The additive effects of sion coefficient; MWi , ADGij and BFi are, respectively, each SNP (k ) on each trait were estimated according to the metabolic weight (computed as body weight 0.75) at the following model: the midpoint of the trial, the ADG during the period, and the BF at the end of the period for individual i ; γ(j), yij = bj + βagei + ui + sik ak + eij , and δ(j), are the corresponding partial regressions coef- ficients nested within batch; and RFIij is the residual feed where yij corresponds to the phenotype of the i th indi- intake of individual i . vidual in the j th batch; bj corresponds to the j th batch Pigs were slaughtered at an approximate age of 190 days effect (4 levels); agei is the covariate of age at slaughter of (average live weight of 122 kg). After recording live body individual i , and β is the corresponding regression coef- weight and BF in vivo, pigs were slaughtered according ficient; ui is the  infinitesimal genetic effect of individual i , with u ∼ N 0, Gσu2 , where G is the genomic relation-  to a commercial protocol. Carcass weight (CW) was reg- istered, the killing out percentage (KO, %) computed, ship matrix (GRM) calculated using the filtered autoso- and lean percentage (LEAN) was inferred based on fat mal SNPs based on the methodology of Yang et al. [18], and muscle thickness data measured with an AutoFOM and σu2 is the additive genetic variance; sik is the genotype ultrasound device. Finally, the percentage of intramus- (coded as 0,1,2) of individual i for the k th SNP, and ak is cular fat content (IMF) was determined on a sample of the allele substitution effect of SNP k on the trait under study; and eij is the residual term. tance (NIT, Infratec ® 1625, Tecator Hoganas, Sweden), gluteus medius (GM) muscle by near infrared transmit- Following a previously published procedure [33], in as described in [30, 31]. order to build the AWM, we retained the SNPs that were associated (nominal p-value < 0.05) with RFI (target trait) Genotype information and/or with three or more of the remaining eight pheno- Genome-wide SNP genotyping of the 350 Duroc pigs was types. Considering both cis-action and the extent of link- performed using the porcine SNP60 BeadChip (Illumina, age disequilibrium (LD) in pigs, only the SNPs that were San Diego, CA), which contains 62,163 SNPs. SNPs with located within or less than 10 kb from the nearest anno- a minor allele frequency (MAF) lower than 5%, a rate of tated gene (Sscrofa11.1 assembly) were retained. Next, missing genotypes higher than 10%, and those that did we used the z-scores of the estimated allele substitution not conform to Hardy–Weinberg expectations (thresh- effects of the SNPs to build the AWM matrix of dimen- old set at a p value of 0.001) were filtered out. We also sion number of retained SNPs (rows) per number of traits excluded SNPs that did not map to the porcine reference (columns). Hierarchical clustering of traits from allele genome (Sscrofa11.1 assembly) and that were located substitution effects of SNPs was estimated and visualized on the X chromosome. After these filtering steps, we using the ‘hclust’ R function. The AWM matrix allowed obtained a subset of 30,096 SNPs that were used in the us to explore both correlations between traits based on GWAS and in the expression GWAS (eGWAS). Quality SNP additive effects (column-wise) and gene-by-gene control of genotypes and the filtering steps were per- interactions (row-wise). Gene-by-gene interactions were formed with the GenomeStudio (Illumina) and PLINK predicted using the partial correlations and information [32] programs, respectively. theory (PCIT) algorithm [37]. The resulting network was used to identify potential regulators by focusing on the transcription factors (TF) within the network. Once the Association weight matrix (AWM) and network analysis TF were identified, we applied an information lossless The association weight matrix (AWM), which has been approach [34] to the inferred co-association network, applied in previous studies [33–35], allows gene co-asso- which explored the combinations (trios and quartets) of ciation networks with regulatory significance to be gener- TF that spanned most of the network topology with min- ated by combining GWAS results with network inference imum redundancy. algorithms. In this study, we used the AWM approach to identify candidate genes and regulators that underlie FE. Nine phenotypes were considered in the analysis: RFI, Gene expression data FCR, ADFI, ADG, CW, IMF, KO, BF and LEAN. Since Gene expression profiles in muscle were obtained for the most accepted FE measure is RFI, it was set as the key three groups of 35 pigs with high, medium, and low lipid Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 4 of 17 metabolism, respectively, which were selected from the difference between the logarithms of mean expression lev- 350 pigs in the population by a principal components els in LFE and HFE pigs, i.e. a positive FC corresponds to analysis, as described in Ref. [31]. The objective, here, higher expression in the LFE group compared to their HFE was to select animals with divergent profiles regarding counterparts. Genes were considered to be DE when |FC| fat deposition and lipid metabolism, but these animals was higher than 1.5 and the q-value lower than 0.05, after covered the whole spectrum of the population’s vari- adjusting for multiple-testing with the false discovery rate ability regarding RFI and the other analysed traits. GM method [42]. muscle samples from these 104 pigs were immediately The same dataset of the expression level of 7007 genes on collected after slaughter and snap-frozen in liquid nitro- 20 extreme pigs for RFI was used in a multivariate frame- gen until storage at − 80 °C. Total RNA was extracted by work to perform a sparse partial least squares discrimi- using the acid/phenol method [38] implemented in the nant analysis (sPLS-DA) [43] in order to identify a subset Ribopure isolation kit (Ambion, Austin, TX). The mRNA of genes that discriminated samples according to RFI clas- expression profile of each sample was characterized by sification (HFE vs. LFE). In a first step, we determined the hybridization to the GeneChip Porcine Genome Array classification error rates for the sample group assignation (Affymetrix Inc., Santa Clara, CA), which includes 23,998 with respect to the number of variables (genes) selected for probes, in two laboratories (66 and 38 samples in each each component. The classification error rates and optimal laboratory). Details about RNA isolation and microar- number of variables to select for each component are rep- ray hybridisation procedures are in [31]. Pre-processing, resented in Additional file 1: Figure S2. The classification background correction, normalization, and log-transfor- performance of the final model was assessed in a fivefold mation of the expression data were performed by com- cross-validation repeated 500 times, by a function of the puting a robust multi-array average (RMA) per probe maximum distance between overall misclassification error [39]. The gene intensity significance level for detecting rate and balanced error rate that considers the proportion expressed probes was calculated by using the MAS 5.0 of incorrectly classified samples weighted by the number of algorithm [39]. Control probes and probes for which the samples per group. expression level was lower than the detection threshold Finally, a regularized canonical correlation analysis in 75% of the pigs were discarded for further analyses. (rCCA) was performed using the whole expression data- The remaining probes were mapped to the Sus scrofa set for all 104 individuals. The rCCA is an unsupervised genome assembly (Sscrofa11.1) using the Biomart data- multivariate approach to identify subsets of canonical vari- base available at the Ensembl repository (https​://www. ables that maximize the correlation between two datasets ensem​bl.org/bioma​rt/martv​iew/). Expression values for X and Y, of sizes (n × p) and (n × q), respectively [43]. In probes that mapped to the same locus were averaged in our analysis, X was the matrix of phenotypes for all 104 order to obtain a global estimate of transcript expression animals for the traits that were most directly associated at the gene level. Probes that failed to map to a known with FE (i.e. RFI, FCR, ADG and ADFI) and Y was a matrix gene were also removed. The effect of the laboratory with gene expression values for all 104 animals. The shrink- where microarrays were assayed on gene expression lev- age method was used to tune the regularization param- els was estimated to be an additive effect, thus, a system- eters λ1 and λ2, with shrinkage values of λ1 = 0.121138 atic laboratory effect was included in subsequent analyses and λ2 = 0.128887. Instead of considering all genes that of the gene expression data. were included in the first canonical component (CC1), we decided to apply a more conservative approach and keep Differential expression and multivariate analyses as candidates only the genes for which the correlation of expression data between gene expression and FE related traits (RFI, FCR, Twenty pigs with extreme RFI were selected for differential ADG and ADFI) was higher than 0.29 (median + 2*SD). expression (DE) analysis: 10 highly feed-efficient (HFE) pigs Expression‑based genome‑wide association studies with low RFI and 10 lowly feed-efficient (LFE) pigs with (eGWAS) high RFI. Offspring of four sires were present in both the HFE and LFE groups and all 20 animals had different dams. The aim of the eGWAS was to identify SNPs that are asso- The DE analysis between HFE and LFE animals was done ciated with the expression of the candidate genes identi- by following the limma-trend pipeline recommendations fied in the AWM approach, as well as genes reported by at [40, 41], fitting a model with batch and laboratory effects, least two of the three following methods mentioned before: in addition to FE-group. The limma’s empirical Bayes pro- rCCA, sPLS-DA and DE. The association of each SNP with cedure was modified to incorporate a mean–variance trend the expression of each gene was estimated by fitting the fol- that models the relationship between variance and gene lowing model using the GCTA software [36]: signal intensity. Fold-change (FC) was computed as the yijl = bj + ll + ui + sik ak + eijl , Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 5 of 17 where yijl corresponds to the gene expression of individ- was run for 100,000 iterations with a burn-in of 10,000 ual i raised in batch j and processed in laboratory l ; bj , rounds, and then saving one out each 10 samples. ui , a and eijl are as defined in the previous GWAS model; The accuracy for predicting RFI phenotype based on and ll is the fixed effect of the laboratory where the the different sources of genetic (pedigree) and genomic microarray was assayed (2 levels). After multiple testing (all SNPs, and AWM or AWM-eGWAS subsets of SNPs) adjustment, the cut-off for a significant association at the information was assessed by cross-validation. The cross- whole-genome level was established at a q-value ≤ 0.05. validation scheme comprised 20 random replicates. In each replicate, the whole dataset was split randomly into Gene functional classification and canonical pathway training and validation datasets that contained approxi- analyses mately 88 and 12% of records, respectively. The training Functional classification and pathway analyses of the dataset was used to predict the genetic additive effects list of candidate genes were carried out using the Inge- of SNPs by solving the mixed model equations (blupf90) nuity Pathways Analysis software (IPA; Ingenuity Sys- with variance components estimates that were obtained tems, http://www.ingen​uity.com). Significance levels for with the complete dataset. Subsequently, phenotypes in enrichment of each canonical pathway in the list of can- the validation dataset were predicted from model solu- didate genes were calculated using Fisher’s exact test and tions obtained in the training set and the prediction accu- the resulting p-values were corrected for multiple-test racy was defined as the correlation coefficient between using the Benjamini and Hochberg algorithm [42]; the predicted and observed records in the validation dataset. cut-off for considering an enrichment as significant was The accuracy of each model/SNP subset was computed established at a corrected p-value < 0.05. by averaging correlations across replicates. Proportion of phenotypic variance explained by the SNPs Results and their prediction accuracy Phenotypes analysed The proportion of phenotypic variance of RFI that was In our study, we focused on the genetic regulation of explained by the SNPs identified in the previous analyses FE, but we considered nine traits that are directly or was estimated by a Bayesian Gibbs sampling approach, indirectly associated with FE, such as growth rate, feed using the gibbs2f90 program included in the Blupf90 intake, carcass traits, and fat deposition. The analysed package [36, 44]. We performed three estimations of pigs belonged to a commercial Duroc line that is used to genomic variance by considering different subsets of produce highly cured products and characterised by its SNPs: (1) SNPs that were selected by the AWM proce- high IMF depot. Summary statistics of the phenotypes dures (829 SNPs); (2) the former AWM-SNPs plus the considered for the analysed population (350 pigs) are in eSNPs identified by the eGWAS (1078 SNPs); and (3) all Table 1. Due to being a residual term, RFI has a mean of SNPs that passed quality controls (~ 31K SNPs). The sin- zero, whereas the mean for FCR indicated that, on aver- gle-step method [45] implemented in Blupf90 was used. age, the analysed Duroc pigs consumed 3.16 kg of feed Subsequently, we also performed a pedigree-based esti- for 1 kg of growth. The values for the remaining pheno- mation of RFI heritability, in order to define a baseline for types are consistent with the general characteristics of comparing the proportion of variance explained by the this Duroc line: carcasses weighing ~ 95 kg with a KO SNPs. percentage of ~ 75%; high subcutaneous fat deposition All estimates of variance components were obtained (mean BF of 24 mm) and intramuscular compartments with the following animal model: (mean IMF in GM reaches values > 5%), and a low lean RFIij = bj + βagei + ui + eij , percentage in the carcass of 40.8% compared with other where all terms are defined as previously described for breeds [47]. It should be noted that age and weight at the GWAS. For the Bayesian implementation, the addi- slaughter were substantially greater than the values nor- tive genetic effect ui was assumed to follow a normal dis- mally reached in commercial production conditions. tribution with a mean of zero and different (co)variance Phenotypic means in the two groups selected for extreme matrices depending on the analysis: (1) Kσu2 for genomic- RFI (HFE and LFE) are also in Table 1. The two groups based estimations, where K is the genomic relation- diverged significantly for all FE phenotypes, particularly ship matrix computed for the different SNP subsets, as in the classification criterion RFI (mean values of − 0.30 described in [46]; and (2) Aσu2 for pedigree-based estima- vs. 0.26 kg/d in the HFE vs. LFE groups) but also in FCR tion, where A is the pedigree-based numerator relation- (mean values of 2.72 vs. 3.30 for the HFE vs. LFE pigs, ship matrix. A prior uniform distribution was assumed p-value < 10−5). These two groups of animals also differed for batch and age effects. The Gibbs sampler algorithm in feed consumption and production traits, with the HFE Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 6 of 17 Table 1 Mean (standard deviation) of the analysed residual term from a regression function that encom- phenotypes in the population of 350 individuals passed ADG and BF and was, therefore, independent and in the two extreme groups (10 individuals each) from these two phenotypes. The strongest phenotypic for RFI, denoted as HFE and LFE (high and low feed relationships were observed between feed intake (ADFI), efficiency groups, respectively), and the significance growth (ADG), and CW ­(rP ranging from 0.71 to 0.81). (p-value) of differences between HFE and LFE groups Subcutaneous fat deposition (BF) showed a moderate to Phenotype Mean (SD) Mean Mean p-valuea high correlation with all traits, except with RFI and KO, N = 350 HFE LFE whereas IMF had a low to moderate, but always signifi- (N = 10) (N = 10) cant, correlation with all analysed traits, including RFI. RFI (kg/day) 0.00 (0.17) − 0.30 0.26 < 10−11 LEAN displayed negative associations with fat deposition FCR (kg feed/kg 3.16 (0.31) 2.72 3.30 < 10−5 traits (BF and IMF), but also with ADG, ADFI and RFI. gain) Finally, the phenotypic correlations of KO with the other ADFI (kg/day) 2.81 (0.37) 2.31 3.06 < 10−4 traits were negligible, except with CW. ADG (kg/day) 0.89 (0.11) 0.85 0.93 n.s. CW (kg) 94.8 (10.6) 89.0 98.6 < 10−1 KO (%) 74.8 (2.4) 73.4 76.6 < 10−1 Association matrix, gene co‑association network, BF (mm) 24.0 (4.9) 22.4 24.5 n.s. and potential regulators for FE Lean (%) 40.8 (4.5) 42.6 36.8 < 10−1 The GWAS results served as the basis for the AWM IMF (%) 5.23 (2.06) 4.83 7.68 < 10−2 approach. After the SNP selection process, 829 SNPs RFI, residual feed intake; FCR, feed conversion ratio; ADFI, average daily feed were retained to build the AWM co-association matrix intake; ADG, average daily gain; BF, back fat thickness; CW, carcass weight; KO, killing out percentage; LEAN, lean percentage; IMF, intramuscular fat content; with RFI plus the eight other analysed phenotypes. Con- n.s., non-significant differences between groups sistently with the AWM procedure, most selected SNPs a Student t test of LS means for group effect (620 out of 829) were associated with the key pheno- type RFI and with on average two other traits, while the remaining 209 SNPs were associated with at least three of animals displaying lower feed intake (0.70 kg less per the other traits but not with RFI. Annotation of the 829 day), smaller weight at slaughter (carcasses 9 kg lighter), SNPs that were selected in the AWM procedure identi- higher lean content (42.6 vs. 37.0%) and lower IMF (4.8 fied 879 genes, since several SNPs were annotated to vs. 7.7%) than the LFE pigs. more than one gene. The list of selected SNPs and the Correlation coefficients between the analysed pheno- corresponding annotated genes is in Additional file 2: types ­(rP) are in Table 2. High but less than 1 phenotypic Table S1. Among the genes that were associated with sev- correlations ­(rP = 0.68) were observed between the two eral of the analysed traits, we would like to mention those FE indicators RFI and FCR. Both these traits were sig- related to the nervous system, such as sidekick cell adhe- nificantly correlated with feed intake ­(rP = 0.46 and 0.43 sion molecule 1 (SDK1), neuronal pentraxin 1 (NPTX1), for RFI and FCR, respectively). Conversely, only FCR neuronal guanine nucleotide exchange factor (NGEF), and was associated with ADG and BF, because RFI was the catenin delta 2 (CTNND2). We also identified genes that Table 2 Estimates of phenotypic (above the diagonal) and AWM-genomic (below the diagonal) correlations between the analysed phenotypes and their significance RFI FCR ADFI ADG CW KO BF Lean IMF RFI 0.677*** 0.461*** − 0.004 ns − 0.001 ns 0.055 ns 0.025 ns − 0.079* 0.162** *** *** *** ns ns ** ns FCR 0.763 0.425 − 0.276 − 0.075 0.025 0.191 − 0.050 0.138** *** *** *** *** ns *** *** ADFI 0.435 0.177 0.738 0.709 0.087 0.687 − 0.231 0.302*** ns ** ** *** ns *** ** ADG 0.011 − 0.261 0.771 0.811 0.044 0.569 − 0.191 0.205*** *** *** *** ** *** CW − 0.061 + − 0.282 0.720 0.927 0.338 0.600 − 0.141 + 0.188** *** *** ns *** ** ns ns KO − 0.204 0.422 0.039 0.219 0.422 0.044 0.039 0.029+ ns ns *** *** *** *** BF 0.015 − 0.038 0.618 0.740 0.675 0.060+ − 0.279 0.352*** *** *** *** *** ** ** ** Lean − 0.304 − 0.241 − 0.415 − 0.295 − 0.172 0.194 − 0.462 − 0.182*** *** *** *** *** *** + *** *** IMF 0.242 0.162 0.481 0.402 0.325 − 0.014 0.607 − 0.451 RFI, residual feed intake; FCR, feed conversion ratio; ADFI, average daily feed intake; ADG, average daily gain; BF, back fat thickness; CW, carcass weight; KO, killing out percentage; LEAN, lean percentage; IMF, intramuscular fat content; ns non-significant *** p < 0.0001; **p < 0.01; *p < 0.05; + p < 0.1 Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 7 of 17 are linked to the immune system such as tyrosine-protein kinase JAK1 (JAK1) and DnaJ heat shock protein family (Hsp40) member C6 (DNAJC6), which are among the genes that were associated with a large number of traits. Using the AWM columns, correlations between traits were also computed based on the estimates of the stand- ardized allele substitution effects of the set of 829 SNPs across traits. Sample size and experimental design did not allow us to properly estimate genetic correlations between traits; using pedigree information only, the uncertainty regions for these parameters covered the whole parametric space. Therefore, phenotypic cor- relations (Table 2) were used as a basis to compare the relationships between traits based on SNP effects. The correlation coefficients obtained from AWM were in general consistent with but larger than the estimated phenotypic correlations between traits (Table 2). For instance, the AWM-derived correlation coefficient between RFI and FCR increased to 0.76 (vs. ­rP = 0.68), and that between BF and IMF was 60% higher than the corresponding phenotypic correlation (­rAWM = 0.61 vs. ­rP = 0.39). In contrast, the correlation between ADG and ADFI captured by AWM was only slightly higher than that the estimated phenotypic correlation (­rAWM = 0.77 vs. ­rP = 0.74). The map of associations between phe- notypes based on the AWM additive effects of selected SNPs was not identical to the observed phenotypic rela- tionships between traits, and this result is reflected by the hierarchical tree cluster in Fig. 1, in which traits that Fig. 1 Hierarchical tree cluster of the nine analysed phenotypes obtained from the standardized additive effects of 829 SNPs are mainly associated with FE, i.e. RFI and FCR, cluster identified by the AWM procedure. RFI, residual feed intake; FCR, feed together. This FE block was associated with a second conversion ratio; ADFI, average daily feed intake; ADG, average daily block that included the remaining traits distributed in gain; BF, back fat thickness; CW, carcass weight; KO%, killing out two groups: Group 1 included growth, feed intake, and percentage; LEAN, lean percentage; IMF, intramuscular fat content fat deposition, which was further subdivided in two blocks: ADFI + ADG + CW and BF + IMF; and Group 2 included LEAN and KO, which clustered together in interactions, each). Other TF genes from the RUNX spite of their negligible phenotypic correlation. and FOX (forkhead box protein) families, jointly with The network derived from the co-association analy- the bromodomain PHD finger transcription factor sis obtained by PCIT gathered 829 nodes, which corre- (BPTF), Wilms tumor 1 (WT1), general transcription sponded to the formerly mentioned 829 SNPs and that factor IIF subunit 2 (GTF2F2) or COP9 signalosome were connected by 57,718 significant edges that repre- subunit 5 (COPS5) genes, can also be included in the sented the significant interactions occurring between list of top TF based on their number of interactions, them. Sixty-eight SNPs mapped to genes classified as since they all showed a number of connections larger TF (indicated as TF in Additional file 2: Table S1) and than average (139). were thus candidates as potential regulators of the net- The information lossless approach allowed the iden- work. As expected, TF genes were the most interacting tification of the combinations of TF that spanned genes in the network, i.e. the hub in the co-associa- most of the network topology with minimum redun- tion network topology. Among these genes, ubiquitin dancy. These combinations (trios or quartets) of TF are carboxyl-terminal hydrolase 16 (USP16), runt-related in Table 3. The top trios of regulators spanned a net- transcription factor 1 (RUNX1), and SWI/SNF-related work that gathered between 519 and 521 nodes (out matrix-associated actin-dependent regulator of chro- of 829). When a combination of four TF was consid- matin A2 (SMARCA2) accumulated the largest num- ered, the regulated network expanded from 603 to 613 ber of interactions within the network (more than 220 nodes. Both types of combinations, trios and quartets, Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 8 of 17 Table 3 Combination of regulators (trios and quartets) variance component estimation, only 20% of the RFI vari- in the list of 829 SNPs identified by the AWM procedure ance was captured by the additive genetic effects of SNPs. that spanned most of the network topology obtained Finally, the pedigree-based estimation yielded a herit- by PCIT with minimum redundancy ability of 0.51 for RFI (Table 4). Although their estimated Combinations Regulators Number errors are remarkably large, these latter figures allow us of target genes to define a baseline to assess the relevance of the propor- Trios COPS5-GTF2H5-RUNX1 521 tion of RFI variance explained by the SNPs identified by COPS5-HDAC4-RUNX1 520 AWM. ESR1-SMARCA2-USP16 519 The expected accuracy for predicting RFI using the Quartets COPS5- ESR1-GTF2F2- USP16 613 SNP panel selected by AWM was assessed by cross-vali- COPS5- ESR1-RUNX1- WT1 603 dation. The correlation between actual and predicted RFI using the effects of the 829 SNPs was 0.65 (Table 4). This COPS5- GTF2F2- HDAC4- USP16 603 value is far from that obtained when RFI was predicted by considering all 31K SNP genotypes or when using only pedigree information in a best linear unbiased predic- included at least one of the top TF mentioned above tion (BLUP) procedure (prediction accuracies of 0.12 and (USP16, RUNX1, SMARCA2 or COPS5; and in quartets 0.20, respectively). Finally it is worth mentioning that the also GTF2F2 or WT1), combined with other TF that mean prediction accuracy obtained with the 620 (out of had a smaller number of but less redundant interac- 829) SNPs that were directly associated to RFI was equal tions, such as general transcription factor IIH subunit 5 to 0.61, a value that was slightly but not significantly (GTF2H5), histone deacetylase 4 (HDAC4) or estrogen lower than that obtained when all 829 SNPs identified by receptor 1 (ESR1), which was present in three out of six AWM were considered (Table 4). top combinations (Table 3). The co-association network that spanned the maximum number of nodes (613) was linked to the combination of COPS5, ESR1, GTF2F2 Gene expression and feed efficiency (DE and multivariate and USP16 TF. Other relevant TF included in the net- analyses) work (although not taking part in the top combinations Two groups of animals selected for extreme RFI (HFE or regulators), that had a large number of interactions and LFE, 10 individuals for each) were used in the DE (> 100) and that displayed associations with FE were analysis of muscle. After quality control, this analysis the LIM homeobox protein 4 (LHX4), mastermind-like was performed with the expression levels of 7007 genes protein 2 (MAML2), and transcription factor 7 like 2 that were expressed in the GM muscle. A total of 991 (TCF7L2) genes. were differentially expressed between the HFE and LFE groups (|FC| > 1.5; q-value < 0.05), of which 892 genes RFI prediction based on the SNPs identified with AWM had higher expression in the LFE samples (Fig. 2 and We evaluated the usefulness of the set of SNPs identified Additional file 3: Table S2). The estrogen related recep- in the AWM procedure to predict RFI phenotype. Fol- tor gamma (ESRRG​) gene showed the largest differ- lowing a single-step animal model, the additive effects of ences in expression between the two groups of extreme the 829 SNPs yielded by the AWM procedure explained FE (FC = 4.94; q-value = 8.2*10−5). Regulator genes that about 61% of the phenotypic variance of RFI (Table 4); are involved in energy metabolism and IMF, such as please note the low error estimated for this parameter. PPARGC1A, were also upregulated in the LFE group, Noteworthy, when all available (31K) SNPs were used for as well as the heat shock protein DNAJC2 gene, which Table 4 Proportion of variance in residual feed intake explained by additive genetic effects (heritability) and correlations between observed and predicted records (prediction accuracy) based on either pedigree or genomic data using three sets of SNPs Pedigree All SNPs AWM SNPs AWM + eGWAS SNPs a Heritability 0.51 ± 0.21 0.20 ± 0.11 0.61 ± 0.06 0.66 ± 0.06 Prediction ­accuracyb 0.20 (0.13) 0.12 (0.11) 0.65 (0.15) 0.60 (0.09) All SNPs, the complete 31K SNP dataset; AWM SNPs, SNPs identified by the association weight matrix (AWM) approach (829 SNPs); AWM + eGWAS SNPs, SNPs identified by AWM and eGWAS (1078 SNPs) a Marginal posterior mean ± marginal posterior standard deviation b Average (standard deviation) across 20 replicates of cross-validation Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 9 of 17 pigs belonged to pathways related to protein ubiqui- tination, valine and isoleucine degradation, IL-6 and IL-8 signalling, glucocorticoid receptor, and estrogen receptor signalling. Regarding the multivariate analyses that were per- formed with sPLS-DA, the first principal component (PC1) combined the expression pattern of 200 genes, which explained 24% of the total variance in gene expres- sion and allowed a clear discrimination between the two extreme feed efficiency groups (Fig. 3). Supporting this accurate classification, a low balanced error rate ranging from 0.13 (PC1) to 0.08 (PC2) was observed (see Addi- tional file 1: Figure S3). The list of 200 genes included in the PC1 and their corresponding contribution to PC1 is in Additional file 4: Table S3. Importantly, 17 out of the 18 most discriminating genes (LAMC2, RF00278, MFSD1, ULK2, HP1BP3, ZNF276, BBS4, RBBP6, PER3, PRKAA2, TCEA3, NFATC3, MTUS1, SNTB1, EIF3F, Fig. 2 Volcano plot of differential gene expression between the two groups of animals with high and low feed efficiency (HFE and LFE, SLC16A5, and UHRF1BP1) were also identified in the DE respectively) and/or in the rCCA analyses described below. The rCCA procedure allowed us to explore the gene expression and phenotype joint (co)variation by using the drives the cellular response to heat stress. In contrast, whole expression dataset (104 individuals). This analy- nuclear receptor genes involved in myogenesis such sis yielded 350 genes (see Additional file 5: Table S4) as retinoid X receptor gamma (RXRG) were upregu- that were included in the first canonical component lated in the more efficient group (HFE). In general, we (CC1), which, in turn, displayed correlations higher observed that the genes that were upregulated in LFE than 0.29 with the FE traits. Among these genes, it is a b 2195 1910 1468 1456 2716 8693 1528 2152 1624 8694 1452 2151 2734 2552 2805 2546 2737 2550 2807 2709 Fig. 3 a Representation of the samples belonging to animals with high (green) and low (red) feed efficiency according to the two first components of gene expression levels obtained with the discriminant analysis performed by sPLS-DA. b Clustering of samples obtained with the 200 genes included in the first component Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 10 of 17 worth highlighting the ESRRG​ gene, which was previ- ously shown to have the highest DE between groups of extreme FE. Other relevant genes identified within CC1 were ALDH1A2, NSMAF, ARMC6, PI15, CD163, C4BPA, PLCXD2, MYOD1, PIK3R1, NNAT, ALDH18A1, ARHGAP29, TMEM158, ART3, SYBU, TMEM98, and CCDC71. Genes included in the CC1 were also used to establish a correlation network between gene expression and phenotype variation for the four traits that were most associated with FE: RFI, FCR, ADG and ADFI. We found that most of the genes from CC1 were correlated with RFI, and that the muscle expression pattern for these genes allowed the clustering of RFI and FCR (Fig. 4), as was the case when co-associated SNPs were used. In summary, 57 common genes were reported by all three approaches, i.e. DE, sPLS-DA, and rCCA, and this number increased to 221 when genes detected by at least two approaches were considered (Fig. 5). Among the genes for which expression was associ- ated with FE indicator traits (RFI and FCR) and with Fig. 5 Overlapping genes from the three approaches used for gene expression analysis the other production traits, several TF genes were identified as potential regulators of FE, including ESRRG​, ZNF473, NFATC3, RXRG, PPARGC1A, NFK- BIZ, TCEA1, CDCA7L, ZFP64, LZTFL1, RBL2, and CBFB. Finally, it is worth mentioning the concordance between the results obtained with different approaches in gene expression analyses. This way, a strong and positive correlation (r = 0.64) was found between the loadings of the 57 common genes in the first principal component (PC1) and the first canonical component (CC1) that were obtained with the sPLS-DA and rCCA Fig. 6 Correlation between the loading factors for the gene expression of the common genes identified with the rCCA and sPLS-DA procedure procedures, respectively (Fig. 6). Similarly, genes that were up-regulated in the LFE group showed negative loadings in PC1 and were positively correlated with RFI and FCR, whereas genes with a positive weight on PC1 Fig. 4 Clustering of RFI, FCR, ADG and ADFI obtained with the 350 were up-regulated in the HFE group and negatively cor- genes identified by rCCA​ related with RFI and FCR. Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 11 of 17 Synapc Long Term Depression Reelin Signaling in Neurons tRNA Splicing Role of BRCA1 in DNA Damage Response Amyotrophic Lateral Sclerosis Signaling Branched-chain α-keto acid Dehydrogenase Complex Aldosterone Signaling in Epithelial Cells Cardiac β-adrenergic Signaling Synapc Long Term Potenaon Sperm Molity CREB Signaling in Neurons tRNA Charging Gustaon Pathway Relaxin Signaling Acetyl-CoABiosynthesis I (Pyruvate Dehydrogenase Complex) Dopamine-DARPP32 Feedback in cAMP Signaling GPCR-Mediated Nutrient Sensing in Enteroendocrine Cells Wnt/Ca+ pathway Protein Ubiquinaon Pathway Renoate Biosynthesis I RAR Acvaon 0 1 2 3 4 5 6 -log(p-value) Fig. 7 Biological pathways that are over-represented in the list of candidate genes for feed efficiency. The x-axis represents the −log(p-value) Candidate genes associated to FE and functional by taking into account their mRNA levels in muscle (only classification 290 of the genes identified by AWM were expressed in Our results obtained in the structural and functional GM). genomic analyses were combined to compile a list of can- The eGWAS showed significant associations at the didate genes, which included the 879 genes that were pre- genome-wide level (q-value < 0.05) between 269 SNPs sent in the co-association network built by AWM and the (eSNPs) and the expression of 16 of the 497 genes 221 genes that were identified in at least two of the three (Table 5 and Additional file 7: Table S6). These 269 eSNPs approaches used for gene expression analyses. Functional were distributed across 30 intervals (Table 5) that were annotation (Fig. 7 and Additional file 6: Table S5) showed located on nine chromosomes (SSC for Sus scrofa chro- that these genes belong to a wide variety of biological mosome, SSC1, SSC2, SSC6, SSC7, SSC8, SSC9, SSC12, processes related to behaviour, immunity, nervous sys- SSC14, and SSC17). The largest number of associations tem signalling, and neurotransmitters. was found for the ENSSSCG00000024596 and ENS- SSCG00000032907 genes, which included 13.6% (39) Identification of eSNPs (eGWAS) for the candidate genes and 11.1% (32) of the eSNPs, respectively. In addition, we An eGWAS was performed to identify SNPs that were found that 40.5% of these SNPs (117 eSNPs) regulated associated with expression (eSNPs) of candidate genes gene expression in cis- (i.e. the genome position of eSNP for FE. We hypothesized that SNPs that regulate the and the target gene map differed by ± 1 Mb) and the expression of genes that are associated with FE may con- remaining 172 eSNPs were associated with the expres- tribute to improve the molecular-based prediction of RFI. sion level of genes located either on a distant genome Thus, the final goal was to evaluate the potential increase region (distance > 1 Mb) or on another chromosome. in predictive ability of RFI that was derived from adding Finally, to verify whether the inclusion of these these polymorphisms to the panel of selected SNPs iden- putative regulatory markers improves RFI predic- tified by AWM. The eGWAS was performed on 497 genes tion, the 269 eSNPs that affect genes associated to selected from the aforementioned list of candidate genes FE were added to the former SNP panel of 829 SNPs. Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 12 of 17 Table 5 Genome-wide eQTL for 497 candidate genes expressed in the gluteus medius (GM) muscle of Duroc pigs Ensembl ID SSC N SNP Region (Mb) q-value B δ ± SE A1 MAF ENSSSCG00000004064 1 27 rs81289665 5.69–10.57 4.34E−03 1.30E−02 − 0.45 ± 0.09 G 0.42 ENSSSCG00000004071 1 6 rs81278169 12.26–13.05 8.75E−04 8.75E−04 0.51 ± 0.09 A 0.40 ENSSSCG00000014432 2 1 rs80929229 149.45–149.45 3.71E−03 2.97E−02 − 0.47 ± 0.10 A 0.10 ENSSSCG00000014432 2 9 rs81290831 150.27–151.35 2.72E−03 1.91E−02 − 0.46 ± 0.09 A 0.11 ENSSSCG00000005976 4 6 rs80965603 14.33–15.13 1.55E−02 3.10E−02 0.24 ± 0.05 C 0.28 ENSSSCG00000005976 4 1 rs81261174 29.68–29.68 4.40E−02 3.08E−01 0.23 ± 0.05 C 0.29 ENSSSCG00000032907 4 3 rs80955389 89.79–89.95 4.01E−02 1.00E+00 − 0.24 ± 0.06 C 0.48 ENSSSCG00000032907 4 32 rs80977488 90.00–99.43 4.25E−04 4.25E−04 0.35 ± 0.06 G 0.36 ENSSSCG00000032907 4 4 rs80889119 103.32–103.47 4.01E−02 1.00E+00 0.31 ± 0.08 G 0.45 ENSSSCG00000033497 6 4 rs81388022 68.99–69.13 1.85E−02 7.41E−02 0.34 ± 0.07 C 0.23 ENSSSCG00000027659 6 25 rs81389074 80.81–82.56 8.02E−04 1.36E−02 − 0.30 ± 0.06 A 0.15 ENSSSCG00000038994 6 19 rs81389074 80.81–82.15 6.32E−03 1.14E−01 0.85 ± 0.18 A 0.15 ENSSSCG00000027659 6 2 rs81309503 111.39–111.40 3.98E−02 1.00E+00 − 0.26 ± 0.06 G 0.15 ENSSSCG00000025087 6 2 rs81393449 156.12–156.13 4.27E−02 8.55E−02 − 0.42 ± 0.09 A 0.09 ENSSSCG00000003844 6 3 rs81345502 157.59–158.65 1.13E−06 1.13E−06 0.30 ± 0.05 G 0.43 ENSSSCG00000003844 6 1 rs81333128 162.74–162.74 2.37E−02 7.10E−02 0.21 ± 0.04 A 0.34 ENSSSCG00000001455 7 12 rs339209635 23.40–26.06 9.64E−03 1.98E−02 0.60 ± 0.12 A 0.49 ENSSSCG00000024596 8 1 rs81237638 37.96–37.96 1.60E−02 6.55E−01 − 0.33 ± 0.08 G 0.33 ENSSSCG00000024596 8 22 rs81402024 83.19–89.67 1.58E−02 1.23E−01 − 0.37 ± 0.08 C 0.37 ENSSSCG00000024596 8 39 rs81402439 90.31–99.95 1.58E−02 2.99E−01 0.34 ± 0.08 G 0.33 ENSSSCG00000024596 8 11 rs81310686 100.02–101.17 1.58E−02 4.81E−01 0.34 ± 0.08 A 0.31 ENSSSCG00000024596 9 1 rs81420190 14.32–14.32 4.59E−02 1.00E+00 − 0.30 ± 0.08 C 0.24 ENSSSCG00000015129 9 4 rs81408352 24.50–28.89 6.62E−03 1.28E−01 − 0.27 ± 0.06 G 0.10 ENSSSCG00000015129 9 28 rs81318704 35.15–39.98 6.62E−03 1.39E−01 − 0.25 ± 0.06 A 0.12 ENSSSCG00000015129 9 10 rs81410722 45.90–47.65 2.64E−05 5.28E−05 − 0.25 ± 0.04 G 0.17 ENSSSCG00000015129 9 3 rs81226448 54.24–54.49 2.97E−02 7.72E−01 − 0.15 ± 0.03 C 0.28 ENSSSCG00000015129 9 1 rs81412811 65.10–65.10 4.49E−02 1.00E+00 − 0.19 ± 0.05 C 0.12 ENSSSCG00000017895 12 9 rs81436670 50.91–53.86 4.50E−07 1.36E−06 0.40 ± 0.06 G 0.33 ENSSSCG00000023829 14 1 rs80875259 31.27–31.27 1.24E−06 3.74E−02 − 0.36 ± 0.07 G 0.42 ENSSSCG00000014791 17 1 rs81466737 47.44–47.44 9.17E−07 2.76E−02 0.33 ± 0.07 A 0.22 SSC, Sus scrofa chromosome; N, number of SNPs significantly associated with the gene under study; SNP, corresponds to the most significant associated SNP; Region (Mb), region containing the SNPs significantly associated with the gene under study; P-value, nominal P-value of association between gene expression and SNP; q-value, q-value calculated with a false discovery rate approach; B, Bonferroni corrected P-values; δ, effect size of the marker and its standard error (SE); ­A1, minor allele; MAF, minor allele frequency The estimated proportion of RFI variance explained Discussion by these 1078 SNPs (Table 4) did not increase sig- In this study, we used a systems-based approach that nificantly (0.66 ± 0.06) compared to the variance combines SNP co-association and muscle gene expres- explained by the 829 SNPs (0.61 ± 0.06) identified by sion analyses to identify candidate genes, biological AWM. The results from the subsequent cross-vali- pathways, and potential regulators of FE in pigs. An addi- dation analysis showed that the correlation between tional goal of our work was to develop a panel of SNPs observed and predicted RFI was slightly lower when that could predict FE phenotypes. The results obtained the 269 eSNPs were added (dropped from 0.65 to 0.60; from both sources of information allowed us to recapit- Table 4), which led us to conclude that there was no ulate the known biological processes that affect FE (i.e. improvement in RFI prediction accuracy derived from supported by current literature) and to identify novel bio- including these putative regulatory markers in the pre- logical pathways and candidate genes associated with this dictive SNP panel. phenotype. Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 13 of 17 The AWM co-association network that was found to be Combining the genes that were identified in both the associated with the nine evaluated phenotypes included genome and muscle transcriptome analyses resulted in 829 SNPs. The list of candidate loci that directly or indi- an extensive list with more than 1000 candidate genes rectly affected FE included 879 genes, which are involved that were shown to have either a direct or an indirect in a wide variety of biological processes. Among these, effect on FE traits. According to their functional annota- genes that displayed associations with a high proportion tion, these genes are involved in a broad set of biologi- of the analysed traits were functionally related to either cal processes, indicating that the molecular mechanisms the nervous system function or development, as well as to that control FE are highly interconnected and, in turn, the immune system. In agreement with this finding, and are not only associated with energy metabolism, but supporting their potential pleiotropic roles (since they also with immunity, behaviour, and the nervous system. were associated with several traits), it is worth mention- Interestingly, the nervous system pathways included the ing that several of these genes were previously reported melatonin signalling, glutamate receptor, and gustation to be associated with FE traits (SDK1, NGEF, CTNND2 pathways, which seem to play an important role in regu- and JAK1) [17, 48–50], fatness (SDK1, JAK1) [51, 52], lation of feed intake and FE [64–68]. Other pathways that growth (CACHD1, NGEF) [17, 53], or bone mineral den- were enriched in the list of candidate genes have also sity (DNAJC6) [54]. More interestingly, the AWM proce- been previously reported as associated with FE traits [16, dure allowed us to identify a group of TF genes that acted 28, 69], including aldosterone signalling in epithelial cells, as hubs in the topology of the network, which suggests ephrin receptor signalling, relaxin signalling, glycogen that they have a cooperative role in mediating a highly degradation, protein kinase A signalling, axonal guid- inter-connected regulatory cascade that seems pivotal ance signalling, semaphorin signalling in neurons, and for FE in pigs. This group of TF genes included ESR1, RhoGDI signalling pathways. SMARCA2, COPS5, GTF2H5, RUNX1, USP16, TCF7L2, One of the main advantages of our methodological MAML2 and LHX4, which have already been reported to approach lies in the joint interpretation of results from be associated with FE in pigs and other livestock species both structural and functional genomic studies. Indeed, [15, 18, 21, 26, 55–58]. they corroborate the fundamental role of the nervous The association between the skeletal muscle transcrip- system in the regulation of FE and, specifically, of the tome and the traits that were most markedly correlated pituitary gland jointly with the hypothalamus and thy- with FE was also explored, using different approaches. roid axis. First, several genes that are functionally related Results from DE between two groups of animals with to the function and development of the nervous system extreme FE pointed out that less efficient animals have a were included in the AWM network associated to FE higher protein turnover and increased energy expendi- traits, including SDK1, NPTX1, NGEF, and CTNND2. ture compared to their highly efficient counterparts. Moreover, LIM homeobox protein 4 (LHX4), a TF gene Discriminant analysis showed that the gene expression involved in the control of the differentiation and devel- profile of 200 genes (included in PC1) explained 24% of opment of the pituitary gland [70], was included in the the global phenotypic variance in FE-related traits and set of potential regulators detected in the co-association allowed clear discrimination between animals with high network. In addition, in the muscle transcriptome analy- and low RFI. This result allowed us to hypothesize that ses other relevant regulator genes were part of the list of muscle gene expression data can have a predictive ability DE genes, such as RXRG, which is expressed in pituitary for classifying pigs into those having high or low FE, as cells. RXRG-deficient mice display high metabolic rates observed by Piles et al. [59], who analysed liver and gut and resistance to weight gain when fed a high-fat diet gene expression data using machine-learning algorithms. [71], two features which may be explained by interference Based on results from the different transcriptomic analy- with the thyrotrope axis and/or from effects specific to ses, we identified other potential regulators of FE among skeletal muscle [71]. Other genes associated with regu- the genes that were functionally associated with FE. This lation of FE through the hypothalamus and thyroid axis group of potential regulator genes included ESRRG​, were also functionally associated with FE traits in our ZNF473, NFATC3, RXRG, PPARGC1A, NFKBIZ, TCEA1, transcriptomic analyses, such as core-binding factor beta CDCA7L, ZFP64, LZTFL1, RBL2, and CBFB, which were (CBFB) or cells inhibitor zeta (NFKBIZ). previously reported to be related to FE [15, 16, 20, 21, Another relevant observation from the joint inter- 58, 60–63]. Among these genes, it is worth mentioning pretation of the genomic and transcriptomic analyses that ESRRG​ showed the largest differences in expression is the key role of estrogen signalling in the regulation between animals with divergent FE and was among the of FE in pigs based on several findings that support top most relevant genes to explain the global phenotypic this observation. On the one hand, the estrogen recep- variance in FE-related traits in the multivariate analysis. tor 1 (ESR1) gene was among the AWM top regulators Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 14 of 17 of the co-association network. On the other hand, the possibility of an upward bias in the estimated genomic ESRRG​ gene, which encodes an estrogen-related recep- heritability due to the fact that the entire dataset was tor, showed a muscle expression pattern that was clearly used to both selecting the SNPs and performing the vali- associated with FE variance. These results are consistent dation processes. In any case, this bias would not invali- with previous studies at the transcriptomic [15, 20] and date the relevant increase in the explained variance and proteomic [55] levels that reported a link of the ESR1 and prediction accuracy when compared with those obtained ESRRG​ genes with FE in pigs. Estrogens control many with the whole SNP dataset (31K). These results support cellular processes and modulate growth and maintenance the relevance of preselecting SNPs to improve the pre- of the skeleton, as well as the cardiovascular and nervous diction accuracy and that there is little benefit in simply systems [72, 73]. In addition, increased estrogen recep- increasing SNP density, in agreement with previous stud- tor (ERα and ERβ) signalling suppresses energy intake ies [78, 79]. It is also worth noting that prediction accu- and increases energy expenditure [74]. Several studies racy when the whole SNP panel was used was worse than have shed light on the interaction between the estrogen that achieved with only pedigree information, revealing signalling pathway and the nervous system, and provide the uselessness of the whole SNP panel when the training elements to understand the regulatory role of estrogens data is of small size. on FE-related traits. At the level of the central nervous Inclusion of eSNPs in the panel of markers did not system, the hypothalamus is responsible for control- improve the prediction accuracy of RFI, in spite of the ling feed intake (appetite), energy metabolism, and body slight increase in explained RFI variance. Among other weight [74], and it has been reported that estrogens dis- reasons, this may be because of the modest size of the play a nucleus-specific action within the hypothalamus to sample that was used for our mRNA expression analy- modulate energy balance [75]. Noteworthy, the nuclear ses, i.e. 104 pigs, which limited the power of the study. In estrogen receptor ESR1, which is expressed on the pro- addition, most of the genes identified by the AWM pro- opiomelanocortin neurons, acts by suppressing feed cedure were related to nervous system functions, which intake and increasing energy expenditure [74] when it is may be poorly expressed in muscle. This agrees with the activated by the estrogen hormone. Furthermore, cross- limited number of genes that were detected using the talk has been reported between estrogen and the regula- AWM approach and that were expressed in GM, i.e. 290 tion of thyroid hormone-releasing receptor [76, 77]. out of 879. The study of other relevant tissues, such as Finally, from a more practical perspective, our aim was the pituitary gland and hypothalamus, may contribute to to evaluate potential applications of our results in ani- identifying candidate genes for FE and functional vari- mal breeding. Improving FE is one of the most relevant ants that are associated with their activity that could be objectives of the pig industry. Several investigations have useful to predict RFI phenotype. In any case, our results proven the possibility of applying successful selection allow us to conclude that the molecular markers that processes for FE [1–3], but selecting for FE is challenging were identified using the AWM approach present some because obtaining reliable measures of individual feed predictive ability for FE, thus opening the possibility of intake is difficult and costly. Thus, identifying accurate applying more accurate selection processes to improve predictive markers of FE could be of capital relevance FE in pigs with a shorter generation interval, since pre- for genetic selection for FE. The 829 SNPs identified dictions of FE could be obtained early in the life of the by the AWM approach explained a significant propor- animal. tion (61%) of the RFI phenotypic variance, which sug- gests that they have a promising potential for obtaining Conclusions a molecular-based prediction of FE records. This panel An integrative approach that combines different sources of SNPs offered the most favourable scenario regarding of information has made it possible to identify candidate RFI prediction, based on a genomic prediction accu- genes, pathways, and predictors that may be important racy of about 0.65, which was higher than the prediction in the determinism of FE in pigs. A list of genes that may accuracies for FE reported in previous studies on pigs have direct or indirect effects on FE was elaborated by (0.40–0.53) [11, 12], which used larger population sizes combining outputs from genomic and muscle transcrip- and a larger number of SNPs. The estimated proportion tome analyses. These genes are involved in a broad set of of variance explained by our panel of selected SNPs was biological processes that are mainly related to immunity, also greater than the estimated RFI heritability. Previ- behaviour, and nervous system, along with energy metab- ously, a simulation study suggested that selecting causal olism. The set of TF genes that were identified either as SNPs in the panel of predictor markers could increase the regulators of the co-association network or as functionally estimated prediction accuracy above what is expected for related to FE traits are potentially responsible of the mod- a given heritability [78]. However, we cannot rule out the ulation of a highly inter-connected regulatory cascade that Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 15 of 17 seems pivotal for FE in pigs. The list of these main regu- Union H2020 Research and Innovation programme under Marie Skłodowska- Curie Grant (P-Sphere) Agreement No. 6655919. E. Mármol-Sánchez is funded lator genes includes ESR1, ESRRG​, RXRG, PPARGC1A, with a Ph.D. Grant from the Spanish Ministry of Education (FPU15/01733). SMARCA2, COPS5, GTF2H5, RUNX1, USP16, TCF7L2, M. Ballester is recipient of a Ramon y Cajal post-doctoral fellowship (RYC- MAML2, LHX4, ZNF473, NFATC3, NFKBIZ, TCEA1, 2013–12573) from the Spanish Ministry of Economy and Competitiveness. Thanks also to Grant 2017 SGR 1719 from the Agency for Management of CDCA7L, ZFP64, LZTFL1, RBL2, and CBFB. Joint University and Research Grants and the CERCA Programme of the Generalitat interpretation of results from both the structural and de Catalunya. functional genomic studies indicated that the estrogen Availability of data and materials signalling pathway, the pituitary gland, the hypothalamus, Microarray data belonging to selected samples were deposited in the Gene and the thyroid axis play fundamental roles in the regula- Expression Omnibus (GEO) public repository, and are accessible through GEO tion of FE in pigs. The additive effects of a panel of 829 Series Accession Number GSE115484. The phenotypic and genotypic datasets used during the current study are available from the corresponding author on selected SNPs explained 61% of the phenotypic variance reasonable request. of RFI, whereas the prediction accuracy of RFI phenotype based on cross-validation was 0.65. These results offer a Ethics approval and consent to participate The research protocol was approved by the animal care and use committee of promising perspective about the usefulness of molecular the Institut de Recerca i Tecnologia Agroalimentàries (IRTA). approaches to predict FE and open the possibility of more accurate selection processes for improving FE in pigs Consent for publication Not applicable. with a shorter generation interval, since predictions of FE could be obtained early in the life of the animal. Competing interests The authors declare that they have no competing interests. Supplementary information Author details Supplementary information accompanies this paper at https​://doi. 1 Animal Breeding and Genetics Program, Institute for Research and Technol‑ org/10.1186/s1271​1-019-0490-6. ogy in Food and Agriculture (IRTA), Torre Marimon, 08140 Caldes de Montbui, Spain. 2 Department of Animal Genetics, Centre for Research in Agricultural Genomics (CRAG), CSCIC‑IRTA‑UAB‑UB, Campus de LA Universitat Autònoma Additional file 1: Figure S1. Representation of the system-based de Barcelona, 08193 Bellaterra, Spain. 3 Departament de Ciència Animal i dels approach employed to identify candidate genes and predictors of feed Aliments, Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain. efficiency. Figure S2. Classification error rates for each component in the sparse partial least squares discriminant analysis (sPLS-DA); the optimal Received: 13 December 2018 Accepted: 19 August 2019 number of variables to select in each component is indicated as a diamond. Figure S3. Classification performance obtained with the partial least squares discriminant analysis (sPLS-DA). Additional file 2: Table S1. Annotation of the 829 SNPs that were selected in the association weight matrix (AWM) procedure. References Additional file 3: Table S2. Description of genes identified as differen‑ 1. Gilbert H, Bidanel JP, Gruand J, Caritez JC, Billon Y, Guillouet P, et al. tially expressed between the high and low feed efficiency groups (HFE Genetic parameters for residual feed intake in growing pigs, with empha‑ and LFE). sis on genetic relationships with carcass and meat quality traits. J Anim Sci. 2007;85:3182–8. Additional file 4: Table S3. Description of the genes identified in the 2. Cai W, Casey DS, Dekkers JCM. Selection response and genetic parame‑ partial least squares discriminant analysis (sPLS-DA). ters for residual feed intake in Yorkshire swine. J Anim Sci. 2008;86:287–98. Additional file 5: Table S4. Description of the genes identified in the 3. Hoque MA, Kadowaki H, Shibata T, Oikawa T, Suzuki K. Genetic param‑ regularized canonical correlation analysis (rCCA). eters for measures of the efficiency of gain of boars and the genetic relationships with its component traits in Duroc pigs. J Anim Sci. Additional file 6: Table S5. Functional annotation of the identified 2007;85:1873–9. candidate genes. 4. Dritz SS. Influence of health on feed efficiency. In: Patience JF, editor. Additional file 7: Table S6. Description of the SNPs significantly associ‑ Feed efficiency in swine. Wageningen: Wageningen Academic Publishers; ated with the gene expression according to the expression genome-wide 2012. p. 225–37. association studies (eGWAS). 5. Patience JF, Rossoni-Serão MC, Gutiérrez NA. A review of feed efficiency in swine: biology and application. J Anim Sci Biotechnol. 2015;6:33. 6. Beaulieu AD, Williams NH, Patience JF. Response to dietary digestible Authors’ contributions energy concentration in growing pigs fed cereal grain-based diets. J YRC designed the integrative analysis approach. RQ and MA conceived and Anim Sci. 2009;87:965–76. supervised the experiment. JPS and RQ performed the phenotypic analyses 7. McCormack UM, Curião T, Buzoianu SG, Prieto ML, Ryan T, Varley P, et al. from farm records. YRC performed genomic and bioinformatic analyses. EM Exploring a possible link between the intestinal microbiota and feed and RGP performed gene expression and eSNP analyses, respectively. MB and efficiency in pigs. Appl Environ Microbiol. 2017;83:e00380. MA participated in the interpretation of results. YRC and RQ wrote the initial 8. Yang H, Huang X, Fang S, He M, Zhao Y, Wu Z, et al. Unraveling the fecal draft of the manuscript. All authors contributed to discussion of results and microbiota and metagenomic functional capacity associated with feed edition of the manuscript. All authors read and approved the final manuscript. efficiency in pigs. Front Microbiol. 2017;8:1555. 9. Camarinha-Silva A, Maushammer M, Wellmann R, Vital M, Preuss S, Benne‑ Funding witz J. Host genome influence on gut microbial composition and micro‑ Part of the research presented in this publication was funded by Grants bial prediction of complex traits in pigs. Genetics. 2017;206:1637–44. AGL2013-48742-C2-2-R, AGL2013-48742-C2-1-R, AGL2010-22208-C02-01 and 10. Koch RM, Swiger LA, Chambers D, Gregory KE. Efficiency of feed use in AGL2010-22208-C02-02 awarded by the Spanish Ministry of Economy and beef cattle. J Anim Sci. 1963;22:486–94. Competitiveness. Y. Ramayo-Caldas is financially supported by the European 11. Zhang C, Kemp RA, Stothard P, Wang Z, Boddicker N, Krivushin K, et al. Genomic evaluation of feed efficiency component traits in Duroc pigs Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 16 of 17 using 80 K, 650 K and whole-genome sequence variants. Genet Sel Evol. 33. Fortes MRS, Reverter A, Zhang Y, Collis E, Nagaraj SH, Jonsson NN, et al. 2018;50:14. Association weight matrix for the genetic dissection of puberty in beef 12. Do DN, Janss LLG, Jensen J, Kadarmideen HN. SNP annotation-based cattle. Proc Natl Acad Sci USA. 2010;107:13642–7. whole genomic prediction and selection: an application to feed effi‑ 34. Reverter A, Fortes MRS. Association weight matrix: a network-based ciency and its component traits in pigs. J Anim Sci. 2015;93:2056–63. approach towards functional genome-wide association studies. Methods 13. Gilbert H, Bidanel JP, Billon Y, Lagant H, Guillouet P, Sellier P, et al. Corre‑ Mol Biol. 2013;1019:437–47. lated responses in sow appetite, residual feed intake, body composition, 35. Ramayo-Caldas Y, Renand G, Ballester M, Saintilan R, Rocha D. Multi-breed and reproduction after divergent selection for residual feed intake in the and multi-trait co-association analysis of meat tenderness and other growing pig. J Anim Sci. 2012;90:1097–108. meat quality traits in three French beef cattle breeds. Genet Sel Evol. 14. Do DN, Strathe AB, Ostersen T, Pant SD, Kadarmideen HN. Genome- 2016;48:37. wide association and pathway analysis of feed efficiency in pigs reveal 36. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide candidate genes and pathways for residual feed intake. Front Genet. complex trait analysis. Am J Hum Genet. 2011;88:76–82. 2014;5:307. 37. Reverter A, Chan EKF. Combining partial correlation and an information 15. Reyer H, Oster M, Magowan E, Dannenberger D, Ponsuksili S, Wimmers K. theory approach to the reversed engineering of gene co-expression Strategies towards improved feed efficiency in pigs comprise molecu‑ networks. Bioinformatics. 2008;24:2491–7. lar shifts in hepatic lipid and carbohydrate metabolism. Int J Mol Sci. 38. Chomczynski P, Sacchi N. Single-step method of RNA isolation by acid 2017;18:1674. guanidinium thiocyanate-phenol-chloroform extraction. Anal Biochem. 16. Gondret F, Vincent A, Houée-Bigot M, Siegel A, Lagarrigue S, Causeur D, 1987;162:156–9. et al. A transcriptome multi-tissue analysis identifies biological pathways 39. Gautier L, Cope L, Bolstad BM, Irizarry RA. affy—analysis of Affymetrix and genes associated with variations in feed efficiency of growing pigs. GeneChip data at the probe level. Bioinformatics. 2004;20:307–15. BMC Genomics. 2017;18:244. 40. Sartor MA, Tomlinson CR, Wesselkamper SC, Sivaganesan S, Leikauf GD, 17. Horodyska J, Hamill RM, Varley PF, Reyer H, Wimmers K. Genome-wide Medvedovic M. Intensity-based hierarchical Bayes method improves association analysis and functional annotation of positional candidate testing for differentially expressed genes in microarray experiments. BMC genes for feed conversion efficiency and growth rate in pigs. PLoS ONE. Bioinform. 2006;7:538. 2017;12:e0173482. 41. Law CW, Chen Y, Shi W, Smyth GK. voom: precision weights unlock linear 18. Lkhagvadorj S, Qu L, Cai W, Couture OP, Barb CR, Hausman GJ, et al. Gene model analysis tools for RNA-seq read counts. Genome Biol. 2014;15:R29. expression profiling of the short-term adaptive response to acute caloric 42. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical restriction in liver and adipose tissues of pigs differing in feed efficiency. and powerful approach to multiple testing. J R Stat Soc Series B Meth‑ Am J Physiol-Regul Integr Comp Physiol. 2009;298:R494–507. odol. 1995;57:289–300. 19. Horodyska J, Wimmers K, Reyer H, Trakooljul N, Mullen AM, Lawlor PG, 43. Lê Cao KA, González I, Déjean S. integrOmics: an R package to et al. RNA-seq of muscle from pigs divergent in feed efficiency and unravel relationships between two omics datasets. Bioinformatics. product quality identifies differences in immune response, growth, 2009;25:2855–6. and macronutrient and connective tissue metabolism. BMC Genomics. 44. Aguilar I, Tsuruta S, Masuda Y, Lourenco D, Legarra A, Misztal I. BLUPF90 2018;19:791. suite of programs for animal breeding with focus on genomics. In: 20. Jing L, Hou Y, Wu H, Miao Y, Li X, Cao J, et al. Transcriptome analysis of Proceedings of the 11th World Congress on genetics applied to livestock mRNA and miRNA in skeletal muscle indicates an important network for production: 11–16 February 2018; Auckland; 2018. differential residual feed intake in pigs. Sci Rep. 2015;5:11953. 45. Legarra A, Christensen OF, Aguilar I, Misztal I. Single step, a general 21. Vincent A, Louveau I, Gondret F, Tréfeu C, Gilbert H, Lefaucheur L. Diver‑ approach for genomic selection. Livest Sci. 2014;166:54–65. gent selection for residual feed intake affects the transcriptomic and 46. VanRaden PM. Efficient methods to compute genomic predictions. J proteomic profiles of pig skeletal muscle. J Anim Sci. 2015;93:2745–58. Dairy Sci. 2008;91:4414–23. 22. Morales PE, Bucarey JL, Espinosa A. Muscle lipid metabolism: role of lipid 47. Wood JD, Nute GR, Richardson RI, Whittington FM, Southwood O, Plastow droplets and perilipins. J Diabetes Res. 2017;2017:1789395. G, et al. Effects of breed, diet and muscle on fat deposition and eating 23. Pedersen BK. Muscle as a secretory organ. Compr Physiol. quality in pigs. Meat Sci. 2004;67:651–67. 2013;3:1337–62. 48. Balasubramanian MN, Panserat S, Dupont-Nivet M, Quillet E, Montfort 24. Civelek M, Lusis AJ. Systems genetics approaches to understand complex J, Le Cam A, et al. Molecular pathways associated with the nutritional traits. Nat Rev Genet. 2014;15:34–48. programming of plant-based diet acceptance in rainbow trout following 25. Williams EG, Auwerx J. The convergence of systems and reductionist aAp‑ an early feeding exposure. BMC Genomics. 2016;17:449. proaches in complex trait analysis. Cell. 2015;162:23–32. 49. Zhou N. Using RNA-seq to characterize the biological basis of variation in 26. Do DN, Ostersen T, Strathe AB, Mark T, Jensen J, Kadarmideen HN. feed efficiency in broiler chickens. Master thesis, University of Delaware; Genome-wide association and systems genetic analyses of residual feed 2015. intake, daily feed consumption, backfat and weight gain in pigs. BMC 50. Xi YM, Wu F, Zhao DQ, Yang Z, Li L, Han ZY, et al. Biological mechanisms Genet. 2014;15:27. related to differences in residual feed intake in dairy cows. Animal. 27. Xu Z, Ji C, Zhang Y, Zhang Z, Nie Q, Xu J, et al. Combination analysis of 2016;10:1311–8. genome-wide association and transcriptome sequencing of residual feed 51. Cardoso TF, Cánovas A, Canela-Xandri O, González-Prendes R, Amills M, intake in quality chickens. BMC Genomics. 2016;17:594. Quintanilla R. RNA-seq based detection of differentially expressed genes 28. Ramayo-Caldas Y, Ballester M, Sánchez JP, González-Rodríguez O, Revilla in the skeletal muscle of Duroc pigs with distinct lipid profiles. Sci Rep. M, Reyer H, et al. Integrative approach using liver and duodenum RNA- 2017;7:40005. Seq data identifies candidate genes and pathways associated with feed 52. Davoli R, Luise D, Mingazzini V, Zambonelli P, Braglia S, Serra A, et al. efficiency in pigs. Sci Rep. 2018;8:558. Genome-wide study on intramuscular fat in Italian Large White pig breed 29. Gallardo D, Pena RN, Amills M, Varona L, Ramírez O, Reixach J, et al. Map‑ using the PorcineSNP60 BeadChip. J Anim Breed Genet. 2016;133:277–82. ping of quantitative trait loci for cholesterol, LDL, HDL, and triglyceride 53. Okumura N, Matsumoto T, Hayashi T, Hirose K, Fukawa K, Itou T, et al. serum concentrations in pigs. Physiol Genomics. 2008;35:199–209. Genomic regions affecting backfat thickness and cannon bone circum‑ 30. Quintanilla R, Pena RN, Gallardo D, Cánovas A, Ramírez O, Díaz I, et al. Por‑ ference identified by genome-wide association study in a Duroc pig cine intramuscular fat content and composition are regulated by quanti‑ population. Anim Genet. 2013;44:454–7. tative trait loci with muscle-specific effects. J Anim Sci. 2011;89:2963–71. 54. Rothammer S, Kremer PV, Bernau M, Fernandez-Figares I, Pfister-Schär J, 31. Cánovas A, Quintanilla R, Amills M, Pena RN. Muscle transcriptomic pro‑ Medugorac I, et al. Genome-wide QTL mapping of nine body composi‑ files in pigs with divergent phenotypes for fatness traits. BMC Genomics. tion and bone mineral density traits in pigs. Genet Sel Evol. 2014;46:68. 2010;11:372. 55. Fu L, Xu Y, Hou Y, Qi X, Zhou L, Liu H, et al. Proteomic analysis indicates 32. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, et al. that mitochondrial energy metabolism in skeletal muscle tissue is nega‑ PLINK: a tool Set for whole-genome association and population-based tively correlated with feed efficiency in pigs. Sci Rep. 2017;7:45291. linkage analyses. Am J Hum Genet. 2007;81:559–75. Ramayo‑Caldas et al. Genet Sel Evol (2019) 51:48 Page 17 of 17 56. Graczyk M, Reyer H, Wimmers K, Szwaczkowski T. Detection of the impor‑ 69. Grubbs JK, Fritchen AN, Huff-Lonergan E, Dekkers JCM, Gabler NK, Lon‑ tant chromosomal regions determining production traits in meat-type ergan SM. Divergent genetic selection for residual feed intake impacts chicken using entropy analysis. Br Poult Sci. 2017;58:358–65. mitochondria reactive oxygen species production in pigs. J Anim Sci. 57. Connor EE, Kahl S, Elsasser TH, Parker JS, Li RW, Van Tassell CP, et al. 2013;91:2133–40. Enhanced mitochondrial complex gene function and reduced liver size 70. Pfaeffle RW, Hunter CS, Savage JJ, Duran-Prado M, Mullen RD, Neeb ZP, may mediate improved feed efficiency of beef cattle during compensa‑ et al. Three novel missense mutations within the LHX4 gene are associ‑ tory growth. Funct Integr Genomics. 2010;10:39–51. ated with variable pituitary hormone deficiencies. J Clin Endocrinol 58. Lkhagvadorj S. Effects of selection for low residual feed intake and feed Metab. 2008;93:1062–71. restriction on gene expression profiles and thyroid axis in pigs. PhD 71. Haugen BR, Jensen DR, Sharma V, Pulawa LK, Hays WR, Krezel W, et al. thesis, Iowa State University; 2010. Retinoid X receptor γ-deficient mice have increased skeletal muscle 59. Piles M, Fernandez-Lozano C, Velasco-Galilea M, González-Rodríguez O, lipoprotein lipase activity and less weightgain when fed a high-fat diet. Sánchez JP, Torrallardona D, et al. Machine learning applied to transcrip‑ Endocrinology. 2004;145:3679–85. tomic data to identify genes associated with feed efficiency in pigs. 72. Lee HR, Kim TH, Choi KC. Functions and physiological roles of two types Genet Sel Evol. 2019;51:10. of estrogen receptors, ERα and ERβ, identified by estrogen receptor 60. Hardie LC, VandeHaar MJ, Tempelman RJ, Weigel KA, Armentano LE, knockout mouse. Lab Anim Res. 2012;28:71–6. Wiggans GR, et al. The genetic and biological basis of feed efficiency in 73. Marino M, Galluzzo P, Ascenzi P. Estrogen signaling multiple pathways to mid-lactation Holstein dairy cows. J Dairy Sci. 2017;100:9061–75. impact gene transcription. Curr Genomics. 2006;7:497–508. 61. Mukiibi R, Vinsky M, Keogh KA, Fitzsimmons C, Stothard P, Waters SM, 74. Mauvais-Jarvis F, Clegg DJ, Hevener AL. The role of estrogens in control of et al. Transcriptome analyses reveal reduced hepatic lipid synthesis and energy balance and glucose homeostasis. Endocr Rev. 2013;34:309–38. accumulation in more feed efficient beef cattle. Sci Rep. 2018;8:7303. 75. Morentin PB, González-García I, Martins L, Lage R, Fernández-Mallo D, 62. Li Z, Liu X, Zhang P, Han R, Sun G, Jiang R, et al. Comparative transcrip‑ Martínez-Sánchez N, et al. Estradiol regulates brown adipose tissue tome analysis of hypothalamus-regulated feed intake induced by thermogenesis via hypothalamic AMPK. Cell Metab. 2014;20:41–53. exogenous visfatin in chicks. BMC Genomics. 2018;19:249. 76. Kimura N, Arai K, Sahara Y, Suzuki H, Kimura N. Estradiol transcriptionally 63. Yuan J, Wang K, Yi G, Ma M, Dou T, Sun C, et al. Genome-wide association and posttranscriptionally up-regulates thyrotropin-releasing hormone studies for feed intake and efficiency in two laying periods of chickens. receptor messenger ribonucleic acid in rat pituitary cells. Endocrinology. Genet Sel Evol. 2015;47:82. 1994;134:432–40. 64. Wideman CH, Murphy HM. Constant light induces alterations in mela‑ 77. Vasudevan N, Ogawa S, Pfaff D. Estrogen and thyroid hormone receptor tonin levels, food intake, feed efficiency, visceral adiposity, and circadian interactions: physiological flexibility by molecular specificity. Physiol Rev. rhythms in rats. Nutr Neurosci. 2009;12:233–40. 2002;82:923–44. 65. Prunet-Marcassus B, Desbazeille M, Bros A, Louche K, Delagrange P, 78. Pérez-Enciso M, Rincón JC, Legarra A. Sequence-vs. chip-assisted Renard P, et al. Melatonin reduces body weight gain in Sprague Dawley genomic selection: accurate biological information is advised. Genet Sel rats with diet-induced obesity. Endocrinology. 2003;144:5347–52. Evol. 2015;47:43. 66. Campos CA, Ritter RC. NMDA-type glutamate receptors participate in 79. Lippert C, Quon G, Kang EY, Kadie CM, Listgarten J, Heckerman D. The reduction of food intake following hindbrain melanocortin receptor benefits of selecting phenotype-specific variants for applications of activation. Am J Physiol Regul Integr Comp Physiol. 2015;308:R1–9. mixed models in genomics. Sci Rep. 2013;3:1815. 67. Thoma V, Kobayashi K, Tanimoto H. The role of the gustatory system in the coordination of feeding. eNeuro. 2017. https​://doi.org/10.1523/ ENEUR​O.0324-17.2017. Publisher’s Note 68. Bannai M, Torii K. Detection of dietary glutamate via gut–brain axis. J Springer Nature remains neutral with regard to jurisdictional claims in pub‑ Anim Sci. 2013;91:1974–81. lished maps and institutional affiliations. Ready to submit your research ? Choose BMC and benefit from: • fast, convenient online submission • thorough peer review by experienced researchers in your field • rapid publication on acceptance • support for research data, including large and complex data types • gold Open Access which fosters wider collaboration and increased citations • maximum visibility for your research: over 100M website views per year At BMC, research is always in progress. Learn more biomedcentral.com/submissions