FULL PAPER doi: 10.1111/age.12988 Genome-wide association study for feed efficiency in collective cage-raised rabbits under full and restricted feeding J. P. Sánchez* , A. Legarra†, M. Velasco-Galilea*, M. Piles*, A. Sánchez‡, O. Rafel*, O. González-Rodrı́guez* and M. Ballester* *Animal Breeding and Genetic Program, Institute of Agriculture and Food Research and Technology, Caldes de Montbui 08140, Spain. † GenPhySE, National Institute for Agronomic Research, Castanet-Tolosan 31326, France. ‡Centre for Research in Agricultural Genomics, Campus Universitat Autònoma de Barcelona, Cerdanyola 08193, Spain. Summary Feed efficiency (FE) is one of the most economically and environmentally relevant traits in the animal production sector. The objective of this study was to gain knowledge about the genetic control of FE in rabbits. To this end, GWASs were conducted for individual growth under two feeding regimes (full feeding and restricted) and FE traits collected from cage groups, using 114 604 autosome SNPs segregating in 438 rabbits. Two different models were implemented: (1) an animal model with a linear regression on each SNP allele for growth trait; and (2) a two-trait animal model, jointly fitting the performance trait and each SNP allele content, for FE traits. This last modeling strategy is a new tool applied to GWAS and allows information to be considered from non-genotyped individuals whose contribu- tion is relevant in the group average traits. A total of 189 SNPs in 17 chromosomal regions were declared to be significantly associated with any of the five analyzed traits at a chromosome-wide level. In 12 of these regions, 20 candidate genes were proposed to explain the variation of the analyzed traits, including genes such as FTO, NDUFAF6 and CEBPA previously associated with growth and FE traits in monogastric species. Candidate genes associated with behavioral patterns were also identified. Overall, our results can be considered as the foundation for future functional research to unravel the actual causal mutations regulating growth and FE in rabbits. Keywords candidate gene, feed efficiency, genome-wide association study, growth, pooled records, rabbit, restricted feeding poultry (Emmerson 1997), the selection criterion is tradi- Introduction tionally growth rate or body weight at slaughter, whereas Before the availability of feeding devices for individual in the case of pigs (Sather & Fredeen 1978) it is an index recording of feed intake (FI) of animals raised in groups, based on growth rate and backfat thickness. Genetic breeding programs to improve the feed efficiency (FE) of improvement of FE via indirect selection for these criteria monogastric species achieved important genetic responses has been possible given that they show high heritabilities by using traits that could be measured individually in and moderate correlations with direct measures of FE such animals housed in groups and were genetically correlated as feed conversion ratio (FCR) or residual feed intake. In all with feed efficiency as selection criteria. In this context, FE of the aforementioned species, and in particular in rabbits, should be understood as a general concept that reflects the FE traits, jointly with prolificacy, are the most economically degree of efficacy in the use of feed resources for perfor- relevant ones (Cartuche et al. 2014). In addition, the mance. In the case of rabbits (Estany et al. 1992) and improvement of FE is expected to have positive effects for decreasing the environmental footprint of the rabbit production industry (Gidenne et al. 2017; Cesari et al. Address for correspondence 2018). J.P. Sánchez, Animal Breeding and Genetic Program, Institute of Owing to the non-availability of electronic feeders for Agriculture and Food Research and Technology, 08140 Caldes de individual recording of FI in rabbits housed in groups, Montbui, Spain. Drouilhet et al. (2016) performed a selection experiment to E-mail: juanpablo.sanchez@irta.es improve FE in which animals were housed individually. Accepted for publication 26 June 2020 © 2020 The Authors. Animal Genetics published by 799 John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. 800 Sánchez et al. Despite this strategy offering an interesting framework for produced on a different farm under controlled environmen- understanding FE from a metabolic perspective, it overlooks tal conditions, which produced a better growth rate and a social interactions between cage mates, which are crucial shorter fattening period (30–60 days of age); on this other when animals are raised in groups and especially under farm six animals were kept in each cage. Beyond the farms’ restricted feeding (Piles et al. 2017). This is a common environmental differences and the number of animals per practice in commercial rabbit farms to control digestive cage, the recorded data and management protocols were disorders during fattening (Gidenne et al. 2012). Therefore, the same in both facilities. it could certainly be argued that studies with individually After weaning, kits were randomly assigned to one of two housed rabbits do not reflect the reality of commercial farms feeding regimen (FR) treatments: ad libitum (F) or restricted where animals are reared in groups. However, this exper- to 75% of the ad libitum intake (R). In order to obtain iment provides compelling evidence of favorable genetic homogeneous groups regarding animal size, the kits under responses even when evaluated on animals raised in each FR were assigned to one of two groups based on their collective cages (Garreau et al. 2019). BW at weaning: large size (LS, i.e. kits with BW >700 g) The present study uses data collected from an experiment and small size (SS, i.e. kits with BW ≤700 g). Animals from designed to estimate genetic parameters of FE in animals the same litter were distributed to both FRs. To obtain feed raised in groups. Therefore, the available information restriction to 75% of the ad libitum FI, the amount of feed consists of weekly group records of FI and individual supplied during week 1 was computed as 0.75 times the records of body weight (BW). Piles & Sánchez (2019) average feed intake of kits on F in a specific group j (j = LS studied the data collected in this experiment from a or SS) during the previous week (i.e. i–1), plus 10% quantitative genetic perspective, estimating heritabilities corresponding to the estimated increase in FI as the animals   and genetic correlations of growth and FI on animals raised grew, i.e. FIR;ji ¼ 0:75  ð1 þ 0:1Þ  FIF;jði1Þ for i = 1–5 in groups and under either full or restricted feeding. They and j = LS or SS. also proposed breeding value predictions for FE measures This amount of feed was multiplied by the number of derived from cage-recorded FI and individual growth and animals present in the cage to determine group feed metabolic weight. requirements. The amount of feed for week 1 was computed After the initial rabbit genome assembly (Lindblad-toh from data that were recorded in previous experiments on et al. 2011), Carneiro et al. (2014) released an improved the same production line with animals raised in the same version with the aim of identifying domestication sweeps in season. The actual amounts of feed provided to the rabbits. From the SNPs detected in this study, 200 000 restricted animals were, on average, 75 and 74% the SNPs were included in one array commercialized by ad libitum intake in LS and SS kits respectively. A maximum Affymetrix, opening up possibilities to conduct genomic of two kits per litter were allocated to the same cage in order studies based on dense panels in this species. to minimize collinearity between maternal and pre-weaning The objective of this study was to identify genomic environmental effects and cage effects. regions and potential candidate genes associated with traits In both experimental groups (F and R), the recorded raw involved in the growth and FE of meat rabbits raised in data consisted of weekly individual BW, and for the case of collective cages under different feeding regimes using a the F weekly cage, FI. In both groups, kits were fed the same high-density SNP array for this species. To achieve this standard pellet diet, supplied once per day in a feeder with objective, it was necessary first to create a model for three places, and containing prescribed antibiotics to handling collective cage performance records in the frame- control gut disorders. In both experimental facilities, feed work of GWAS studies, this being a partial objective in our was changed to a standard feed without drugs during the study. last week of fattening. Thus, records from the last week were discarded for the analysis because of the effect that the lack of antibiotics in the feed might have on growth rate, FI Material and methods and the derived FE measures. Therefore, in both farms, the growing period controlled was from 30 to 56 days of age; Animals thus, a total of four weekly individual BW records were The animals used in this study belonged to the Caldes line retained per animal and three weekly group FI measure- (Gómez et al., 2002), and the experiment was conducted at ments were considered per cage. the rabbit farm of the Institute of Agriculture and Food From these raw records, individual average daily gain Research and Technology. They were randomly sampled (ADG) was computed as the regression coefficient of the from four batches during the first semester of 2014 and within-animal BW records on their ages. This was done for from an additional batch in spring 2016. The animals from each FR, obtaining individual ADG on ad libitum (ADGF) or 2014 were raised in a semi-open-air facility and the restricted (ADGR) FR. For the 99 cages on F, individual fattening period was from 30 to 66 days of age. Eight average daily feed intake (ADFI F ) was computed as the total animals were kept in each cage. The batch in 2016 was feed intake of the cage during the whole fattening period © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 GWAS for group feed efficiency traits in rabbits 801 divided by the number of days and number of rabbits Table 2 Basic statistics for the studied traits present in the cage during the whole fattening period. Also, First Third Phenotypic individual average daily feed conversion ratio (ADFCR  F) Trait N Mean SD quartile quartile variance2,2 and individual average daily residual feed intake (ADRFI F)  were computed. The first was the ratio between ADFIF and ADGF (g/ 758 53.2 9.4 50.3 58.8 77.6 day) ADGF cage average, and the second was the residual of a ADGR (g/ 712 35.4 8.0 32.2 40.4 55.0  F on the ADGF cage batch-nested multiple regression of ADFI day) average and the cage average of the mid fattening period  F (g/ ADFI 99 151.4 17.0 141.8 162.7 289.3 day metabolic weight. day)1,1  F [(g/ FCR 99 2.8 0.2 2.7 3.0 0.2 Two datasets were employed in the analyses, one day)/(g/ containing individual growth (ADGF and ADGR) of geno- day)]1,1 typed animals (438) and another including growth from  F (g/ ADRFI 99 0.0 5.9 −3.3 3.4 143.8 genotyped animals as well as from all their non-genotyped day) 1,1 cage mates (438 + 1032). This second dataset also ADGF, Average daily gain in rabbits fed ad libitum; ADGR, average included cage average traits (ADFI  F , ADFCR F and  F average daily feed daily gain in rabbits fed under restriction; ADFI  ADRFIF ). Table 1 shows the number of animals per feeding  F , cage average daily feed intake in rabbits fed ad libitum; FCR regime in cages containing genotyped animals from the five  F , cage average daily conversion ratio in rabbits fed ad libitum; ADRFI batches. In Table 2, raw statistics of the traits under study residual feed intake in rabbits fed ad libitum. 1 Refers to cage traits. are shown. They refer to the animals mentioned before (i.e. 2 Estimated using model 2. For cage average records, the residual genotyped animals and non-genotyped cage mates). variance of the model accounts for the number of animals involved in the mean; thus, these quantities actually represent individual variation. DNA extraction and SNP genotyping The DNA extraction was carried out using the NucleoSpin Prior to a pairwise LD computation, in order to reduce the Tissue (250prep; Macherey-Nagel) commercial kit from computational effort, the genotype file was pruned to liver samples of 438 rabbits collected immediately after retrieve just one SNP every 20 kb; thus, the resolution of slaughter (66 days of age). DNA extracts were sent to an the obtained LD decay pattern was as low as 0.02 Mb. Affymetrix platform to conduct genotyping using the Axiom Rabbit Genotyping Array ‘Axiom_OrCunSNP’ (Thermo Statistical analysis Fisher Scientific), which includes 199 692 variants. Only 161 830 were segregating in our population and, after Two different modeling approaches were adopted to conduct retaining the SNPs mapped in autosomes in the OryCun2.0 the GWASs: assembly and applying standard quality control criteria, 114 604 SNPs were kept for further analysis. The applied Regression analysis on the allele content of each SNP quality control criteria comprised retaining animals having at least 90% of SNPs correctly genotyped, SNPs with less This model was applied to individually recorded traits than 5% missing genotype data and SNPs with a MAF (ADGF and ADGR) and was implemented using QXPAK higher than 5%. The LD (r2) decay pattern from our (Pérez-Enciso & Misztal 2011). The procedure implemented population was assessed using PLINK 1.9 (Chang et al. 2015). with this model is frequently called ‘EMMAX’ (Kang et al. 2010). The general equation of this model fitting the alternative Table 1 Number of individual and cage records per batch and feeding hypothesis is as follows: regime. Genotyped and non-genotyped animals are distinguished for the individual records yijklmno ¼ SNPip  αp þ B j þ Pk þ Ll þ Sm þ cn þ lo þ ai þ eijklmno (1) Individuals Non- where a particular record of a given trait under study Genotyped genotyped Cages (yijklmno ) – ADGF or ADGR – (one at the time) is explained by the effect of the allele content (SNPip : 0, 1, 2 depending on Batch R F R F R F the number of copies of the reference allele) in the pth 1 28 26 68 62 12 11 genomic position of the ith animal, reflected by the 2 41 35 103 84 18 15 regression coefficient at that particular position (αp Þ which 3 58 63 190 209 31 34 represents the allele substitution effect, the effect of the jth 4 35 59 93 124 16 23 batch level (Bj, five levels), the effect of the kth level of the 5 46 47 50 49 16 16 order of parity in which the animal was born (Pk, four F, Animals fed ad libitum; R, animals fed under restriction. levels), the effect of the lth level of size of the litter in which © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 802 Sánchez et al. the animal was born (Ll, seven levels), the effect of the mth In the case of the analysis of group records, the equations level of the type of cage [Sm, two levels, cages with animals involved in the bivariate model were the following: of large body weight at weaning (>700 g) or of low body N  N  n 1 n 1 weight at weaning (≤700 g)], the effect of the nth cage (Cn),  ¼ B j þ Sm þ ∑ yjmn lnk þ ∑ a1,nk þ e1,jmn (3) the oth litter (lo) and the ith breeding value (ai), the last k¼1 N n n k¼1 N n n three being random effects. Thus, each random factor had SNPip ¼ μ þ a2,ip þ e2,ip associated with it a variance component to be estimated. For cage and litter effects, a diagonal structure was assumed Group means, yjmn  F , FCR  , i.e. traits of interests (ADFI  F between the different levels, whereas for the additive genetic  and ADRFIF ), are explained by the effect of the jth batch effect the numerator relationship matrix was used to define level (Bj, five levels), the effect of the mth level of the type of the covariance between the individuals. Finally, a diagonal cage (Sm, two levels) and the averages of litters (lnk) and normal distribution was assumed for the residual term, additive genetic effects (a1,nk ) associated with the Nn eijklmno . individuals in the nth cage. Litter, additive genetic and At each genomic position, two models were fitted by residual effects are random factors, assumed to follow maximum likelihood, including or not (null model) the normal distributions, indexed by their respective (co)vari- regression on the SNP allele content (1). Then, likelihood ances to be estimated using an EM-REML procedure. values at their maximum were compared using likelihood Breeding values for the two traits analyzed at a time were ratio tests. This ratio follows, under the null hypothesis, a assumed to follow a joint multivariate normal distribution chi-squared distribution with 1 degree of freedom; P-values of the following form: were computed from this theoretical distribution. " # ! σ 2a1 σ a1, a2 a≈N 0, A : σ a1, a2 σ 2a2 Bivariate analysis considering each recorded trait, individual growth or cage records, jointly with the allele content of each Similarly, for the residual term, the assumed distribution SNP was the following: " # ! This statistical model was considered as a way to perform σ 2e1 0 GWASs on group mean records. With regard to individual e≈N 0, I : 0 σ 2e2 traits (ADGF and ADGR), the model was the same as that fitting the null hypothesis in the case of regression analysis. In the case of the residual effects, a null covariance was For explaining cage records (ADFI  F , FCR  F and ADRFI F ), a considered between SNP allele content and the performance model similar to that considered by Piles & Sánchez (2019) trait. For the case of the additive genetic effects, this was adopted. The bivariate model was defined by jointly covariance (σ a1;a2 ) under the alternative model was assumed considering, as correlated traits, one performance trait at a to be non-null, representing in this case the association, at a time (either individual or cage average) and the allele genetic level, between breeding values for the trait of content at each SNP. The equation for explaining this allele interest and the SNP genotypes. Under the null model, σ a1, a2 content considered only an overall mean and the additive was set to zero. The REML likelihood values at their genetic effect in addition to a residual term (Legarra & maximum were used to construct likelihood ratio tests Vitezica 2015). The effect of the marker on the trait under allowing exploration of the significance of σ a1, a2 estimates. study was estimated through the genetic covariance of both This was done by computing P-values from the theoretical traits. Legarra & Vitezica (2015) proved that this approach distribution of the ratio under the null hypothesis, a chi- is equivalent, for individual records and complete observa- squared distribution with one degree of freedom. From this tions, to the EMMAX model that is commonly used, the model, the estimated effect for each SNP position was main advantage being the possibility of including missing calculated as a function of the estimated additive genetic genotypes. covariance ( σ ) and the SNP frequency (f p ) (Legarra & a1, a2 The model equations for the bivariate analysis fitting Vitezica, 2015): individually recorded traits were the following, σ a1, a2 yijklmno ¼ B j þ Pk þ Ll þ Sm þ cn þ lo þ a1,i þ e1,ijklmno (2) α¼ p 2∗f p ∗ð1  f p Þ SNPip ¼ μ þ a2,ip þ e2,ip In the two statistical methods, multiple test correction Note that this model was applied to a different set of was performed following the procedure by Storey (2002) to individual records from that employed with model (1). In adjust raw P-values to a positive false discovery rate of 0.05; this second case, we considered individual records from both this was done using the R package ‘qvalue’ (Storey et al. genotyped animals and their non-genotyped cage mates. 2019). The adjustment was done at two different levels: first, at genome-wide level considering all of the tests © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 GWAS for group feed efficiency traits in rabbits 803 conducted; and second, within each chromosome. In the different traits and methods is presented in Manhattan plots second case, thresholds for declaring significance varied (Figs 2–4). across chromosomes, and they were much less strict than Eight chromosomal regions located at OCCs 3, 5 and 21, those applied at genomic level. To define the genomic were declared to be associated with ADGF (Table 3 and regions associated with the analyzed traits, those significant Fig. 2). Two of the five regions on the distal segment of OCC SNPs that were less than 1 Mb apart were grouped in the 3 (102.22–102.37 and 109.07–109.08 Mb) were signifi- same QTL region. This distance threshold was defined based cantly associated with the trait using both modeling on a preliminary assessment of LD as a function of the approaches. The estimated effects of the SNPs with the physical distance between SNPs. strongest association within each region ranged from 3.34 g/day (for a SNP on the region 100.90 Mb–101.11- Mb of OCC 3) when model 2 was used to 6.55 g/day (for an Gene annotation and functional analysis SNP at 107.99 Mb of OCC 3) detected with model 1. The Associated regions were annotated considering 1 Mb effects of the other OCC 3-associated regions were estimated around the previously defined intervals in the rabbit to be close to 4 g/day. For ADGF, 78 ensembl_gene_ids were genome. Gene annotations were retrieved from the Ensembl annotated on the declared QTL regions of OCC 3 (Table S2). Genes 98 Database with the BIOMART software (Smedley et al. One candidate gene, carbonic anhydrase 2 (CA2), was 2015) using the OryCun2.0 reference assembly. Further- identified in the region 100.99–101.11 Mb, whereas two more, functional predictions of the significantly associated candidate genes, NADH:ubiquinone oxidoreductase com- SNPs were performed with VEP software (McLaren et al. plex assembly factor 6 (NDUFAF6) and tumor protein p53 2016). inducible nuclear protein 1 (TP53INP1), were proposed for For functional categorization of the annotated genes, GOs ADGF in the region 109.07–110.88 Mb of OCC 3 (Table 4). were determined using ClueGO version 2.5.0 plug-in of In OCC 5, two significantly associated regions were detected Cytoscape (Bindea et al. 2009). The functions assigned to with model 1; one region comprised a single SNP in position the proposed candidate genes include metabolic, beha- 9.07 Mb and the other comprised two SNPs in the region vioural or immunological pathways. Orthologous human 18.95–18.97 Mb. The magnitude of the strongest estimated gene names were retrieved from the Ensembl Genes 98 effects for these regions was similar to those estimated on Database for functional categorization when a rabbit gene OCC 3 (between 3.5 and 5.5 g/day). In these regions, 19 name was not assigned to the gene stable id. Furthermore, ensembl_gene_ids were annotated (Table S2). Furthermore, information from the Mouse Genome Database (Eppig et al. one promising candidate gene for ADGF alpha-ketoglutarate 2017) and Genecards (Safran et al. 2002) was used to dependent dioxygenase (FTO) was identified at 9.07 Mb in identify gene functions affecting the analyzed phenotypes. OCC 5 (Table 4). Finally, one region in OCC 21 compressing 1.29 Mb (7.17–8.46 Mb) was also associated with ADGF. In this region, 26 SNPs were found to be significantly Results associated with the trait. Within this region, AX- Two modeling approaches were used to conduct a GWAS 147049623, the SNP with the strongest association had on five phenotypic traits related to individual growth and an effect of 3.51 g/day. Remarkably, this SNP was located group FE using 438 rabbits genotyped with AxiomOr- inside an intron of the Ataxin 2 (ATXN2) gene (Table S1) – CunSNP (114 604 SNPs after quality control). one of the four candidate genes (ATXN2, ACAD10, TRAFD1 At genome-wide level, after multiple testing correction, and PTPN11) identified among the 71 ensembl_gene_ids neither of the methods returned significant associations. annotated in this region (Table 4 and Table S2). These However, when multiple test correction was applied within candidate genes contained another 10 SNPs significantly each chromosome, 189 SNPs (Table S1) located in nine associated with ADGF (Table S1). Oryctolagus cuniculus chromosome (OCC) regions (3, 5, 6, 9, ADGR showed significant associations with SNPs on OCCs 12, 13, 16, 17 and 21) were declared as significantly 9, 12, 13 and 17 (Table 3 and Fig. 3). For this trait, associated with any of the five studied traits, i.e. ADGR, however, the two models declared different chromosomal ADGF, ADFI F , FCR  F and ADRFI F. regions as significantly associated with the trait. Model 1 It is important to describe the LD pattern decay (Fig. 1) to declared a QTL region at the proximal region of OCC 13 properly determine that the QTL intervals to be defined (0.40–2.09 Mb) containing 50 significant SNPs. The esti- cover regions in relatively high LD. The LD between regions mated SNP effect having the strongest association (mini- with a distance of 1 Mb was nearly 0.2. Thus, we assumed mum q-value within the region) was 3.41 g/day. Two that significantly associated SNPs within a 1 Mb region candidate genes (RC3H1 and TNFSF18), out of 37 anno- pertain to the same QTL. tated ensembl_gene_ids, were found in this region (Table 4 Table 3 summarizes the significantly associated regions and Table S2). For the same trait, model 2 declared with the traits of interest at chromosome level. In addition, significant signals on OCCs 9 (29.66–31.00 Mb), 12 graphical representation of the results obtained for the (99.88 Mb) and 17 (73.57 Mb–74.16 Mb). The SNP effects © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 804 Sánchez et al. Figure 1 LD (r2) decay pattern. (a) Up to 8 Mb; (b) up to 1 Mb. of the QTL regions on OCCs 9 and 17 were lower Discussion (approximately between 1 and 1.5 g/day) than that To our knowledge, this is the first GWAS for growth and detected on OCC 12, for which an effect of 3.73 g/day was feed efficiency traits performed in a rabbit population using estimated, which was a similar magnitude to the effects of a dense SNP chip panel. Our study, in addition, also the SNPs associated with ADGF. In these regions, 74 introduces a new modeling approach allowing study of the ensembl_gene_ids were annotated (Table S2) and four association of traits recorded as group averages, when not candidate genes were proposed (FEZF2 and PTPRG on all of the individuals in the group have been genotyped. This OCC 9, and LGALS3 and TMEM260 on OCC 17; Table 4). methodology was originally proposed for gene-assisted The GWAS conducted with model 2 for the cage  F , FCR  F and ADRFI F , declared six selection when a certain percentage of the candidates have performance traits, ADFI not been genotyped for the major gene of interest (Legarra & significantly associated regions on OCCs 5, 6, 16 and 21 Vitezica, 2015). Modelling the SNP allele content using (Table 3 and Fig. 4). Region 3.70–3.85 Mb on OCC 5 was  F and comprised 12 animal models has also been proposed as a tool to detect declared to be associated with ADFI low-quality SNPs within the panels (Forneris et al. 2015), significant SNPs; that with the strongest association had a an SNP being declared as erroneously genotyped when its MAF of 0.37 and an estimated effect equal to 0.85 g/day. In heritability estimate is significantly different from 1. With this region, 20 ensembl_gene_ids were annotated (Table  F were identified this work, we extend the scope of application of such models S2) and two candidate genes for ADFI  F , two significant to GWASs, in particular, to GWASs on group average (CEBPA and KCTD15) (Table 4). For FCR performance traits. Previous studies (Zhang et al. 2018) signals were detected on OCCs 6 (26.28–26.44 Mb) and 16 have addressed the problem in the context of experiments (82.86–83.26 Mb). The estimated effects of the SNPs with where the limiting factor is the capability to generate the strongest statistical association within these regions individual phenotypes, but all of the individuals in the were large – 0.47 and 0.52 feed conversion units ((grams of design were genotyped. In this case, it has been shown that feed/day)/(grams of growth/day)) respectively. The most pooling individual records to produce pool phenotypes and significant SNP on OCC 6 had a low MAF (0.06), whereas then explaining these pooled data by the mean genotype of that of the most significant SNPs on the region of OCC 16  F , 63 ensembl_gene_ids the group produced considerable gains in the power of was much higher (0.42). For FCR statistical tests over simple random sampling, i.e. random were annotated (Table S2) and two candidate genes – salt selection of as many individual phenotypes as pools were inducible kinase 1B (putative) (SIK1B) on OCC 6 and defined. This result could be expected as in the analyses of phospholipase A2 group IVA (PLA2G4A) on OCC 16 the pooled phenotypes all of the available genotypes are (Table 4) – were retained.  F , significant included, whereas in the study of a random sample of Finally, for the last studied FE trait, ADRFI individual records only a subset of them are considered, and associations were detected in three regions of OCC 21: this sampling is particularly sensitive to low-frequency 3.89–4.33 and 7.16–7.70 Mb, and one single SNP at markers. Our study, although related to the aforementioned 9.21 Mb. The second region was particularly relevant as it problem, has a completely different motivation: on the one also contained SNPs significantly associated with ADGF. The hand, there is no individual alternative to the group MAF of the most significant SNP within this region was average phenotype recorded, and on the other hand, the 0.37, and its estimated effect was 2.16 g/day. A total of 146 experimental limitations constrain the number of genotyped ensembl_gene_ids were annotated on OCC 21 (Table S2), animals to only a few of those responsible for the group and the same candidate genes as those previously proposed  F (Table 4). average phenotype. In this situation, a much smaller for ADGF in this region were retrieved for ADRFI © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 GWAS for group feed efficiency traits in rabbits 805 Table 3 QTL regions associated with the studied traits according to the two employed methods. Effect estimates and MAF are reported Initial Final position position SNPs in Significant SNPs q- Method Trait Region1,1 OCC2,2 (Mb) (Mb) the region in the region3 SNP name4 Value5 Effect6 MAF7 QXPAK ADGR 11 13 0.40 2.09 90 50 AX- 146990063 3.41 0.0056 (g/ 0.24 day) QXPAK ADGF 2 3 102.22 102.37 11 9 AX- 147016699 3.68 0.0350 (g/ 0.39 day) QXPAK ADGF 3 3 107.99 107.99 1 1 AX- 146983203 6.55 0.0350 (g/ 0.06 day) QXPAK ADGF 4 3 109.07 110.88 111 16 AX- 146982129 4.05 0.0211 (g/ 0.24 day) QXPAK ADGF 5 3 113.46 113.46 1 1 AX- 147140896 3.65 0.0410 (g/ 0.20 day) QXPAK ADGF 6 5 9.07 9.07 1 1 AX- 147010974 3.57 0.0416 (g/ 0.29 day) QXPAK ADGF 7 5 18.95 18.97 2 2 AX- 147049894 5.46 0.0416 (g/ 0.11 day) QXPAK ADGF 8 21 7.17 8.46 67 26 AX- 147102744 3.51 0.0135 (g/ 0.23 day) BI ADGR 9 9 29.66 31.00 66 29 AX- 147167857 1.67 0.0039 (g/ 0.07 day) BI ADGR 10 12 99.88 99.88 0 1 AX- 146984543 3.73 0.0222 (g/ 0.05 day) BI ADGR 12 17 73.57 74.16 29 7 AX- 147012391 0.95 0.0183 (g/ 0.16 day) BI ADGF 1 3 100.99 101.11 3 4 AX- 147009110 3.34 0.0399 (g/ 0.49 day) BI ADGF 2 3 102.22 102.37 10 11 AX- 147016699 3.71 0.0302 (g/ 0.39 day) BI ADGF 4 3 109.07 109.88 58 11 AX- 147097036 3.85 0.0313 (g/ 0.24 day) BI  F FCR 14 6 26.28 26.44 16 10 AX- © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 806 Sánchez et al. Table 3 (Continued) Initial Final position position SNPs in Significant SNPs q- Method Trait Region1,1 OCC2,2 (Mb) (Mb) the region in the region3 SNP name4 Value5 Effect6 MAF7 147140966 0.47 0.0015 [(g/ 0.06 day)/ (g/ day)] BI  F FCR 15 16 82.86 83.26 26 7 AX- 147107945 0.52 0.0482 [(g/ 0.42 day/ (g/ day)] BI  F ADFI 13 5 3.70 3.85 13 12 AX- 147126724 0.85 0.0278 (g/ 0.37 day) BI  F ADRFI 16 21 3.89 4.33 26 8 AX- 147145784 1.14 0.0175 (g/ 0.25 day) BI  F ADRFI 8 21 7.16 7.70 34 15 AX- 147081855 2.16 0.0030 (g/ 0.37 day) BI  F ADRFI 17 21 9.21 9.21 0 1 AX- 147132637 1.34 0.0321 (g/ 0.35 day) 1 Annotated region, match to Table S2. 2 Oryctolagus cuniculus chromosome. 3 Chromosome-wise q-value <0.05. 4 Name of the most significant SNP within the region. 5 Within-region minimum chromosome-wide q-value. 6 Absolute value of the effect of the most significant SNP. 7 MAF of most significant SNP. amount of information is available for the analysis, and declared for ADGR, which is a trait with lower heritability thus, lower power would be expected. We do not formally (Piles & Sanchez 2019), are not the same across the models. assess the efficiency of our proposed model, but its limited One reason for this is that the datasets used by each method power seems obvious. On this regard, Sánchez et al. (2018) are different. In the case of model 1, only records from reported an important reduction in the capability to detect genotyped animals are considered, but with the bivariate simulated QTL regions for one trait (out of three) that is model, records from both genotyped and non-genotyped considered as a group average with respect to the situation animals are jointly considered, and the pedigree is used to in which all traits are studied as individual phenotypes. predict genotypes of non-genotyped animals with records. They implemented a multitrait Bayesian procedure similar As stated, given the available information on the cage to the single-step association methods (Wang et al. 2012), performance records, the expected statistical power was which relies on derivation of SNP effects from genomic low. Thus, in order to allow for a certain degree of signal predictions using a multiple regression in which the SNP detection, the threshold for significance declaration was effects are treated as random factors. deliberately reduced to a chromosome-wide level. In this A possible means of validation of the results from the situation, we have successfully identified 17 chromosomal proposed bivariate model is to analyze those traits individ- regions associated with the analyzed traits. To allow ually recorded with the two approaches. For the case of comparison between traits, the estimated SNP effects within ADGF the results obtained are partially the same, for the regions can be expressed relative to their estimated example regions in OCC 3 are detected with both methods. phenotypic variance (Table 2). To this end, we approxi- However, other regions detected for this trait and all those mated the additive genetic variance associated with each © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 GWAS for group feed efficiency traits in rabbits 807 Figure 2 Manhattan plots for average daily gain recorded in animals fed ad libitum (ADGF) for models (a) QXPAK and (b) BI Figure 3 Manhattan plots for average daily gain recorded in animals fed under restriction (ADGR) obtained for models (a) QXPAK and (b) BI. QTL region, defined in Table 3, by considering the SNP recorded as the group average. Nevertheless, a simple effect with the strongest association (minimum q-value)  F on group average genotype for the regression of group FCR within the region and its frequency. The additive variances SNP AX-147107945 at bp 82858725 in OCC 16, the SNP of the QTL regions in Table 3 represent 5–8, 0.5–8 and with the strongest association within the region (Table 3), 0.5–2% of the phenotypic variance of ADGF, ADGR, and also showed a strong magnitude – 0.20 (0.05) FCR units  F and ADRFI both ADFI  F respectively. The percentage of per unit of change on the cage average genotype (P = 3.38 phenotypic variance explained by the additive genetic effect × 10−5). This means that the expected FCR in a cage with  F is particularly for one of the QTL regions declared for FCR all of the animals heterozygous for this SNP will be 0.20 high – 65% for the region on OCC 16. It could be difficult to units larger than that in a cage with all of the animals propose a validation method for these results free from the homozygous of one type and 0.20 units lower than that in a assumptions in the model for the analysis, because FCR F is cage with all of the animals homozygous of the other type. © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 808 Sánchez et al.  F ), (b) average daily residual feed intake in Figure 4 Manhattan plots using model BI for (a) average daily feed intake in rabbits fed ad libitum (ADFI  F ).  F ) and (c) average daily feed conversion ratio in rabbits fed ad libitum (FCR rabbits fed ad libitum (ADRFI Twenty candidate genes, located in 12 QTL regions, have reported that the genetic correlation between ADGF and been proposed to explain the phenotypic variation of the  F was 0.58 (Piles & Sanchez 2019). ADRFI traits under study; this proposition was done based on their Piles & Sanchez (2019) showed that growth recorded in biological functions. It is worth mentioning the FTO gene, animals fed under restriction (ADGR) is a trait genetically annotated on OCC 5 for ADGF, which has been previously different from growth recorded in animals fed ad libitum. associated with growth traits in rabbits (Zhang et al. 2013; Our results could be said to support this as for ADGR we Zhang et al. 2014). Furthermore, NDUFAF6, which was have declared chromosomal regions different from those also annotated for ADGF in a region of OCC 3, has recently declared for ADGF. Nonetheless, this could be also a simple been described as a candidate gene for growth-related traits consequence of our reduced statistical power. In these in pigs (Ji et al. 2019). regions, candidate genes associated with behavioral pat- It is also relevant to highlight ATXN2, acyl-CoA dehy- terns (FEZF2, PTPRG and LGALS3) or involved in immunity drogenase family member 10 (ACAD10), TRAF-type zinc and/or lipid metabolism (RC3H1, TNFSF18 and TMEM260) finger domain containing 1 (TRAFD1) and protein tyrosine were identified. Finally, it is worth highlighting the CCAAT phosphatase non-receptor type 11 (PTPN11) genes as they enhancer-binding protein alpha (CEBPA) gene, annotated map in a region of OCC 21 with a pleiotropic effect for both on OCC 5 for ADFI  F , which has recently been identified as ADGF and ADRFI  F and have 11 significant SNPs located in an upstream regulator of several differentially expressed their introns. Identifying pleiotropic regions for both ADGF genes down-regulated in adipose tissue of high-feed-effi- and ADRFI F could be considered an unexpected result as ciency pigs (Horodyska et al. 2019).  F is a trait obtained after the phenotypic correction of ADRFI In spite of our loose significance threshold setting, we feel FI by growth and metabolic weight. However, this pheno- relatively confident of having adequately controlled the rate typic correction does not grant a null genetic correlation, of false-positive signals that we have declared. In support of and in fact, in the population under study it has been our results, we have identified some candidate genes that © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 GWAS for group feed efficiency traits in rabbits 809 Table 4 Candidate genes and their associated functions for the QTL regions declared Initial Final Method OCC1 Trait Mb Mb Gene Function Region2 BI 3 ADGF 100.99 101.11 CA2 Respiration and transport of CO2/bicarbonate 1 QXPAK, 3 ADGF 102.22 102.37 — — 2 BI QXPAK 3 ADGF 107.99 107.99 — — 3 QXPAK, 3 ADGF 109.07 110.88 NDUFAF6, TP53INP1 Mitochondrial respiration, oxidative stress 4 BI response QXPAK 3 ADGF 113.46 113.46 — — 5 QXPAK 5 ADGF 9.07 9.07 FTO, AKTIP Growth 6 QXPAK 5 ADGF 18.95 18.97 — — 7 QXPAK, 21  F ADGF,ADRFI 7.16 9.21 ATXN2, ACAD10, TRAFD1, Energy homeostasis, immunity 8,17 BI PTPN11 BI 9 ADGR 29.66 31 FEZF2, PTPRG Behaviour 9 BI 12 ADGR 99.88 99.88 — — 10 QXPAK 13 ADGR 0.4 2.09 RC3H1, TNFSF18 Immunity 11 BI 17 ADGR 73.57 74.16 LGALS3, TMEM260 Circadian rhythm, immunity, lipid metabolism 12 BI 5  F ADFI 3.7 3.85 CEBPA, KCTD15 Energy homeostasis, adipogenesis, feed 13 Behaviour/food intake BI 6  F FCR 26.28 26.44 SIK1B Hormone signaling 14 BI 16  F FCR 82.86 83.26 PLA2G4A Lipid metabolism, inflammatory response 15 BI 21  F ADRFI 3.89 4.33 SELENOM Energy metabolism 16 1 Oryctolagus cuniculus chromosome. 2 Annotated region, match to Table S2. have already been associated with similar traits in other Acknowledgements rabbit and pigs populations. The work conducted by the staff of the Institute of Agriculture and Food Research and Technology’s rabbit Conclusions farm is acknowledged. This study was financed by Feed-a- We have proposed a number of QTL regions linked to the Gene, which has received funding from the European observed variation of the studied traits using a complex Union’s H2020 program under grant agreement no. 63353 statistical model for fitting cage FE and feed intake, jointly and the Spanish projects RTA2014-00015-C2-01, with individual genotypes. To our knowledge, this is the RTA2011-00064-00-00 and RTI2018-097610-R-I00. M. first time this type of statistical model has been used within Ballester is financially supported by a Ramon y Cajal the framework of GWAS studies. The information content contract (RYC-2013-12573) from the Spanish Ministry of on cage average performances is quite limited, thus we have Economy and Competitiveness. The authors declare that reduced the threshold for significance declaration to a there is no conflict of interest regarding the publication of chromosome-wide level. In spite of this loose significance this article. threshold definition, the declared QTL seem to harbor genes that can clearly be regarded as functional candidates for the Availability of data traits of interest. Our results seem to support the idea that the growth of animals fed on restriction is under a different Genotypes can be found here: https://doi.org/10.5281/ze genetic control that that of animals fed ad libitum as we nodo.3611097 have identified different QTL regions for both traits. It is remarkable that genes related to behavioral patterns have References been proposed as candidates for ADGR. Regarding FE, some of the QTL regions that we declared to harbor candidate Bindea G., Mlecnik B., Hackl H. et al. (2009) ClueGO: A Cytoscape genes which are involved in lipid and energy metabolism plug-in to decipher functionally grouped gene ontology and have a pleiotropic effect for both ADGF and ADRFI F : In spite pathway annotation networks. Bioinformatics 25, 1091–3. of these promising results, further functional research is Carneiro M., Rubin C.J., Di Palma F. et al. (2014) Rabbit genome warranted to validate these genes. Overall, our results lay analysis reveals a polygenic basis for phenotypic change during domestication. Science 345, 1074–9. an important foundation for future studies to unravel the Cartuche L., Pascual M., Gómez E.A. et al. (2014) Economic weights underlying genetic bases driving growth and FE regulation in rabbit meat production. World Rabbit Science 22, 165–77. in rabbits. © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810 810 Sánchez et al. Cesari V., Zucali M., Bava L., Gislon G., Tamburini A., Toschi I. using gene content multiple-trait BLUP. Genetics Selection Evolu- (2018) Environmental impact of rabbit meat: The effect of tion 47, 89. production efficiency. Meat Science 145, 447–54. Lindblad-Toh K., Garber M., Zuk O. et al. (2011) A high-resolution Chang C.C., Chow C.C., Tellier L.C., Vattikuti S., Purcell S.M., Lee map of human evolutionary constraint using 29 mammals. J.J. (2015) Second-generation PLINK: rising to the challenge of Nature 478, 476. larger and richer datasets. Gigascience 4(1), 7. http://dx.doi.org/ McLaren W., Gil L., Hunt S.E. et al. (2016) The ensembl variant 10.1186/s13742-015-0047-8 effect predictor. Genome Biology 17, 122. Drouilhet L., Achard C.S., Zemb O. et al. (2016) Direct and Pérez-Enciso M., Misztal I. (2011) Qxpak5: Old mixed model correlated responses to selection in two lines of rabbits selected solutions for new genomics problems. BMC Bioinformatics 12, for feed efficiency under ad libitum and restricted feeding: I. 202. Production traits and gut microbiota characteristics. Journal of Piles M., Sánchez J.P. (2019) Use of group records of feed intake to animal science 94, 38–48. select for feed efficiency in rabbit. Journal of Animal Breeding and Emmerson D.A. (1997) Commercial approaches to genetic selection Genetics 136, 474–83. for growth and feed conversion in domestic poultry. Poultry Piles M., David I., Ramon J., Canario L., Rafel O., Pascual M., Ragab Science 76, 1121–5. M., Sánchez J.P. (2017) Interaction of direct and social genetic Eppig J.T., Smith C.L., Blake J.A., Ringwald M., Kadin J.A., effects with feeding regime in growing rabbits. Genetics Selection Richardson J.E., Bult C.J. (2017) Mouse Genome Informatics Evolution 49, 58. (MGI): resources for mining mouse genetic, genomic, and Safran M., Solomon I., Shmueli O. et al. (2002) GeneCards™ 2002: biological data in support of primary and translational research. towards a complete, object-oriented, human gene compendium. Systems Genetics 2017, 47–73. Bioinformatics 18, 1542–3. Estany J., Camacho J., Baselga M., Blasco A. (1992) Selection Sánchez J.P., Legarra A., Piles M. (2018) Multiple trait single step response of growth rate in rabbits for meat production. Genetics Bayesian GWAS on pooled data. Proceedings of the World Congress Selection Evolution 24, 527. on Genetics Applied to Livestock Production, 11.275 (http:// Forneris N.S., Legarra A., Vitezica Z., Tsuruta S., Aguilar I., Misztal www.wcgalp.org/proceedings/2018/multiple-trait-single-step-ba I., Cantet R.J.C. (2015) Quality control of genotypes using yesian-gwas-pooled-data) heritability estimates of gene content at the marker. Genetics 199, Sather A.P., Fredeen H.T. (1978) Effect of selection for lean growth 675–81. rate upon feed utilization by the market hog. Canadian Journal of Garreau H., Ruesche J., Gilbert H. et al. (2019) Estimating direct Animal Science 58, 285–9. genetic and maternal effects affecting rabbit growth and feed Smedley D., Haider S., Durinck S. et al. (2015) The BioMart efficiency with a factorial design. Journal of Animal Breeding and community portal: an innovative alternative to large, centralized Genetics 136, 168–73. data repositories. Nucleic Acids Research 43, 589–98. Gidenne T., Combes S., Fortun-Lamothe L. (2012) Feed intake Storey J.D. (2002) A direct aproach to false discovery rates. Journal limitation strategies for the growing rabbit: effect on feeding of the Royal Statistical Society: Series B (Statistical Methodology) behaviour, welfare, performance, digestive physiology and 64, 479–98. health: a review. Animal 6, 1407–19. Storey JD, Bass AJ, Dabney A, Robinson D (2019) qvalue: Q-value Gidenne T., Fortun-Lamothe L., Bannelier C. et al. (2017) Direct estimation for false discovery rate control. R package version and correlated responses to selection in two lines of rabbits 2.18.0, http://github.com/jdstorey/qvalue selected for feed efficiency under ad libitum and restricted feeding: Wang H, Misztal I, Aguilar I, Legarra A, Muir W (2012) Genome- III. Digestion and excretion of nitrogen and minerals. Journal of wide association mapping including phenotypes from relatives Animal Science 95, 1301–12. without genotypes. Genetics Research 94, 73–83. Gómez E.A., Rafel O., Ramon J. (2002) The Caldes strain Rabbit Zhang GW, Gao L, Chen SY et al. (2013) Single nucleotide genetic resources in mediterranean countries. Options Méditer- polymorphisms in the FTO gene and their association with ranéennes: Série B. Etudes et Recherches 38, 189–98. growth and meat quality traits in rabbits. Gene 527, 553–7. Horodyska J., Reyer H., Wimmers K., Trakooljul N., Lawlor P.G., Zhang G.W., Jia W., Chen S.Y., Jia X.B., Wang J., Lai S.J. (2014) Hamill R.M. (2019) Transcriptome analysis of adipose tissue Association between the IRS1 and FTO genes regulates body from pigs divergent in feed efficiency reveals alteration in gene weight in rabbits. Gene 548, 75–80. networks related to adipose growth, lipid metabolism, extracel- Zhang W., Liu A., Albert P.S., Ashmead R.D., Schisterman E.F., lular matrix, and immune response. Molecular Genetics and Mills J.L. (2018) A pooling strategy to effectively use genotype Genomics 294, 395–408. data in quantitative traits genome-wide association studies. Ji J., Yan G., Chen D., Xiao S., Gao J., Zhang Z. (2019) An Statistics in Medicine 37, 4083–95. association study using imputed whole-genome sequence data identifies novel significant loci for growth-related traits in a Duroc × Erhualian F2 population. Journal of Animal Breeding and Supporting information Genetics 136, 217–28. Additional supporting information may be found online in Kang H.M., Sul J.H., Zaitlen N.A. et al. (2010) Variance component the Supporting Information section at the end of the article. model to account for sample structure in genome-wide associ- Table S1. Description of the 189 chromosome-wise signif- ation studies. Nature Genetics 42, 348–54. Legarra A., Vitezica Z. (2015) Genetic evaluation with major genes icant SNPs and features annotated to them. and polygenic inheritance when some animals are not genotyped Table S2. List of genes annotated in the QTL regions. © 2020 The Authors. Animal Genetics published by John Wiley & Sons Ltd on behalf of Stichting International Foundation for Animal Genetics, 51, 799–810