Genome‐wide association study for feed efficiency in collective cage‐raised rabbits under full and restricted feeding

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 contribution 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.


Introduction
Before the availability of feeding devices for individual recording of feed intake (FI) of animals raised in groups, breeding programs to improve the feed efficiency (FE) of monogastric species achieved important genetic responses by using traits that could be measured individually in animals housed in groups and were genetically correlated with feed efficiency as selection criteria. In this context, FE should be understood as a general concept that reflects the degree of efficacy in the use of feed resources for performance. In the case of rabbits (Estany et al. 1992) and poultry (Emmerson 1997), the selection criterion is traditionally growth rate or body weight at slaughter, whereas in the case of pigs (Sather & Fredeen 1978) it is an index based on growth rate and backfat thickness. Genetic improvement of FE via indirect selection for these criteria has been possible given that they show high heritabilities and moderate correlations with direct measures of FE such as feed conversion ratio (FCR) or residual feed intake. In all of the aforementioned species, and in particular in rabbits, FE traits, jointly with prolificacy, are the most economically relevant ones (Cartuche et al. 2014). In addition, the 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. 2018).
Owing to the non-availability of electronic feeders for individual recording of FI in rabbits housed in groups, Drouilhet et al. (2016) performed a selection experiment to improve FE in which animals were housed individually.
Despite this strategy offering an interesting framework for understanding FE from a metabolic perspective, it overlooks social interactions between cage mates, which are crucial when animals are raised in groups and especially under restricted feeding (Piles et al. 2017). This is a common practice in commercial rabbit farms to control digestive disorders during fattening (Gidenne et al. 2012). Therefore, it could certainly be argued that studies with individually housed rabbits do not reflect the reality of commercial farms where animals are reared in groups. However, this experiment provides compelling evidence of favorable genetic responses even when evaluated on animals raised in collective cages (Garreau et al. 2019).
The present study uses data collected from an experiment designed to estimate genetic parameters of FE in animals raised in groups. Therefore, the available information consists of weekly group records of FI and individual records of body weight (BW). Piles & Sánchez (2019) studied the data collected in this experiment from a quantitative genetic perspective, estimating heritabilities and genetic correlations of growth and FI on animals raised in groups and under either full or restricted feeding. They also proposed breeding value predictions for FE measures derived from cage-recorded FI and individual growth and metabolic weight.
After the initial rabbit genome assembly (Lindblad-toh et al. 2011), Carneiro et al. (2014) released an improved version with the aim of identifying domestication sweeps in rabbits. From the SNPs detected in this study, 200 000 SNPs were included in one array commercialized by Affymetrix, opening up possibilities to conduct genomic studies based on dense panels in this species.
The objective of this study was to identify genomic regions and potential candidate genes associated with traits involved in the growth and FE of meat rabbits raised in collective cages under different feeding regimes using a high-density SNP array for this species. To achieve this objective, it was necessary first to create a model for handling collective cage performance records in the framework of GWAS studies, this being a partial objective in our study.

Animals
The animals used in this study belonged to the Caldes line (Gómez et al., 2002), and the experiment was conducted at the rabbit farm of the Institute of Agriculture and Food Research and Technology. They were randomly sampled from four batches during the first semester of 2014 and from an additional batch in spring 2016. The animals from 2014 were raised in a semi-open-air facility and the fattening period was from 30 to 66 days of age. Eight animals were kept in each cage. The batch in 2016 was produced on a different farm under controlled environmental conditions, which produced a better growth rate and a shorter fattening period (30-60 days of age); on this other farm six animals were kept in each cage. Beyond the farms' environmental differences and the number of animals per cage, the recorded data and management protocols were the same in both facilities.
After weaning, kits were randomly assigned to one of two feeding regimen (FR) treatments: ad libitum (F) or restricted to 75% of the ad libitum intake (R). In order to obtain homogeneous groups regarding animal size, the kits under each FR were assigned to one of two groups based on their BW at weaning: large size (LS, i.e. kits with BW >700 g) and small size (SS, i.e. kits with BW ≤700 g). Animals from the same litter were distributed to both FRs. To obtain feed restriction to 75% of the ad libitum FI, the amount of feed supplied during week 1 was computed as 0.75 times the average feed intake of kits on F in a specific group j (j = LS or SS) during the previous week (i.e. i-1), plus 10% corresponding to the estimated increase in FI as the animals grew, i.e. FI R;ji ¼ 0:75 Â ð1 þ 0:1Þ Â FI F;j iÀ1 ð Þ À Á for i = 1-5 and j = LS or SS.
This amount of feed was multiplied by the number of animals present in the cage to determine group feed requirements. The amount of feed for week 1 was computed from data that were recorded in previous experiments on the same production line with animals raised in the same season. The actual amounts of feed provided to the restricted animals were, on average, 75 and 74% the ad libitum intake in LS and SS kits respectively. A maximum of two kits per litter were allocated to the same cage in order to minimize collinearity between maternal and pre-weaning environmental effects and cage effects.
In both experimental groups (F and R), the recorded raw data consisted of weekly individual BW, and for the case of the F weekly cage, FI. In both groups, kits were fed the same standard pellet diet, supplied once per day in a feeder with three places, and containing prescribed antibiotics to control gut disorders. In both experimental facilities, feed was changed to a standard feed without drugs during the 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 and the derived FE measures. Therefore, in both farms, the growing period controlled was from 30 to 56 days of age; thus, a total of four weekly individual BW records were retained per animal and three weekly group FI measurements were considered per cage.
From these raw records, individual average daily gain (ADG) was computed as the regression coefficient of the within-animal BW records on their ages. This was done for each FR, obtaining individual ADG on ad libitum (ADG F ) or restricted (ADG R ) FR. For the 99 cages on F, individual average daily feed intake ( ADFI F ) was computed as the total feed intake of the cage during the whole fattening period divided by the number of days and number of rabbits present in the cage during the whole fattening period. Also, individual average daily feed conversion ratio ( ADFCR F ) and individual average daily residual feed intake ( ADRFI F ) were computed. The first was the ratio between ADFI F and ADG F cage average, and the second was the residual of a batch-nested multiple regression of ADFI F on the ADG F cage average and the cage average of the mid fattening period day metabolic weight.
Two datasets were employed in the analyses, one containing individual growth (ADG F and ADG R ) of genotyped animals (438) and another including growth from genotyped animals as well as from all their non-genotyped cage mates (438 + 1032). This second dataset also included cage average traits ( ADFI F , ADFCR F and ADRFI F ). Table 1 shows the number of animals per feeding regime in cages containing genotyped animals from the five batches. In Table 2, raw statistics of the traits under study are shown. They refer to the animals mentioned before (i.e. genotyped animals and non-genotyped cage mates).

DNA extraction and SNP genotyping
The DNA extraction was carried out using the NucleoSpin Tissue (250prep; Macherey-Nagel) commercial kit from liver samples of 438 rabbits collected immediately after slaughter (66 days of age). DNA extracts were sent to an Affymetrix platform to conduct genotyping using the Axiom Rabbit Genotyping Array 'Axiom_OrCunSNP' (Thermo Fisher Scientific), which includes 199 692 variants. Only 161 830 were segregating in our population and, after retaining the SNPs mapped in autosomes in the OryCun2.0 assembly and applying standard quality control criteria, 114 604 SNPs were kept for further analysis. The applied quality control criteria comprised retaining animals having at least 90% of SNPs correctly genotyped, SNPs with less than 5% missing genotype data and SNPs with a MAF higher than 5%. The LD (r 2 ) decay pattern from our population was assessed using PLINK 1.9 (Chang et al. 2015).
Prior to a pairwise LD computation, in order to reduce the computational effort, the genotype file was pruned to retrieve just one SNP every 20 kb; thus, the resolution of the obtained LD decay pattern was as low as 0.02 Mb.

Statistical analysis
Two different modeling approaches were adopted to conduct the GWASs: Regression analysis on the allele content of each SNP This model was applied to individually recorded traits (ADG F and ADG R ) and was implemented using QXPAK (Pérez-Enciso & Misztal 2011). The procedure implemented with this model is frequently called 'EMMAX' (Kang et al. 2010).
The general equation of this model fitting the alternative hypothesis is as follows: where a particular record of a given trait under study (y ijklmno ) -ADG F or ADG R -(one at the time) is explained by the effect of the allele content (SNP ip : 0, 1, 2 depending on the number of copies of the reference allele) in the pth genomic position of the ith animal, reflected by the regression coefficient at that particular position (α p Þ which represents the allele substitution effect, the effect of the jth batch level (B j , five levels), the effect of the kth level of the order of parity in which the animal was born (P k , four levels), the effect of the lth level of size of the litter in which  the animal was born (L l , seven levels), the effect of the mth level of the type of cage [S m , two levels, cages with animals of large body weight at weaning (>700 g) or of low body weight at weaning (≤700 g)], the effect of the nth cage (C n ), the oth litter (l o ) and the ith breeding value (a i ), the last three being random effects. Thus, each random factor had associated with it a variance component to be estimated. For cage and litter effects, a diagonal structure was assumed between the different levels, whereas for the additive genetic effect the numerator relationship matrix was used to define the covariance between the individuals. Finally, a diagonal normal distribution was assumed for the residual term, e ijklmno .
At each genomic position, two models were fitted by maximum likelihood, including or not (null model) the regression on the SNP allele content (1). Then, likelihood values at their maximum were compared using likelihood ratio tests. This ratio follows, under the null hypothesis, a chi-squared distribution with 1 degree of freedom; P-values were computed from this theoretical distribution.
Bivariate analysis considering each recorded trait, individual growth or cage records, jointly with the allele content of each SNP This statistical model was considered as a way to perform GWASs on group mean records. With regard to individual traits (ADG F and ADG R ), the model was the same as that fitting the null hypothesis in the case of regression analysis. For explaining cage records ( ADFI F , FCR F and ADRFI F ), a model similar to that considered by Piles & Sánchez (2019) was adopted. The bivariate model was defined by jointly considering, as correlated traits, one performance trait at a time (either individual or cage average) and the allele content at each SNP. The equation for explaining this allele content considered only an overall mean and the additive genetic effect in addition to a residual term . The effect of the marker on the trait under study was estimated through the genetic covariance of both traits.  proved that this approach is equivalent, for individual records and complete observations, to the EMMAX model that is commonly used, the main advantage being the possibility of including missing genotypes.
The model equations for the bivariate analysis fitting individually recorded traits were the following, Note that this model was applied to a different set of individual records from that employed with model (1). In this second case, we considered individual records from both genotyped animals and their non-genotyped cage mates.
In the case of the analysis of group records, the equations involved in the bivariate model were the following: Group means, y jmn , i.e. traits of interests ( ADFI F , FCR F and ADRFI F ), are explained by the effect of the jth batch level (B j , five levels), the effect of the mth level of the type of cage (S m , two levels) and the averages of litters (l nk ) and additive genetic effects (a 1,nk ) associated with the N n individuals in the nth cage. Litter, additive genetic and residual effects are random factors, assumed to follow normal distributions, indexed by their respective (co)variances to be estimated using an EM-REML procedure.
Breeding values for the two traits analyzed at a time were assumed to follow a joint multivariate normal distribution of the following form: Similarly, for the residual term, the assumed distribution was the following: In the case of the residual effects, a null covariance was considered between SNP allele content and the performance trait. For the case of the additive genetic effects, this covariance (σ a1;a 2 ) under the alternative model was assumed to be non-null, representing in this case the association, at a genetic level, between breeding values for the trait of interest and the SNP genotypes. Under the null model, σ a1,a2 was set to zero. The REML likelihood values at their maximum were used to construct likelihood ratio tests allowing exploration of the significance of σ a1,a2 estimates. This was done by computing P-values from the theoretical distribution of the ratio under the null hypothesis, a chisquared distribution with one degree of freedom. From this model, the estimated effect for each SNP position was calculated as a function of the estimated additive genetic covariance ( σ a1,a2 ) and the SNP frequency (f p ) : In the two statistical methods, multiple test correction was performed following the procedure by Storey (2002) to adjust raw P-values to a positive false discovery rate of 0.05; this was done using the R package 'qvalue' (Storey et al. 2019). The adjustment was done at two different levels: first, at genome-wide level considering all of the tests conducted; and second, within each chromosome. In the second case, thresholds for declaring significance varied across chromosomes, and they were much less strict than those applied at genomic level. To define the genomic regions associated with the analyzed traits, those significant SNPs that were less than 1 Mb apart were grouped in the same QTL region. This distance threshold was defined based on a preliminary assessment of LD as a function of the physical distance between SNPs. For functional categorization of the annotated genes, GOs were determined using ClueGO version 2.5.0 plug-in of Cytoscape (Bindea et al. 2009). The functions assigned to the proposed candidate genes include metabolic, behavioural or immunological pathways. Orthologous human gene names were retrieved from the Ensembl Genes 98 Database for functional categorization when a rabbit gene name was not assigned to the gene stable id. Furthermore, information from the Mouse Genome Database (Eppig et al. 2017) and Genecards (Safran et al. 2002) was used to identify gene functions affecting the analyzed phenotypes.

Results
Two modeling approaches were used to conduct a GWAS on five phenotypic traits related to individual growth and group FE using 438 rabbits genotyped with AxiomOr-CunSNP (114 604 SNPs after quality control).
At genome-wide level, after multiple testing correction, neither of the methods returned significant associations. However, when multiple test correction was applied within each chromosome, 189 SNPs (Table S1) located in nine Oryctolagus cuniculus chromosome (OCC) regions (3, 5, 6, 9, 12, 13, 16, 17 and 21) were declared as significantly associated with any of the five studied traits, i.e. ADG R , ADG F , ADFI F , FCR F and ADRFI F . It is important to describe the LD pattern decay (Fig. 1) to properly determine that the QTL intervals to be defined cover regions in relatively high LD. The LD between regions with a distance of 1 Mb was nearly 0.2. Thus, we assumed that significantly associated SNPs within a 1 Mb region pertain to the same QTL. Table 3 summarizes the significantly associated regions with the traits of interest at chromosome level. In addition, graphical representation of the results obtained for the different traits and methods is presented in Manhattan plots (Figs 2-4).
Eight chromosomal regions located at OCCs 3, 5 and 21, were declared to be associated with ADG F (Table 3 and Fig. 2). Two of the five regions on the distal segment of OCC 3 (102.22-102.37 and 109.07-109.08 Mb) were significantly associated with the trait using both modeling approaches. The estimated effects of the SNPs with the 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 SNP at 107.99 Mb of OCC 3) detected with model 1. The effects of the other OCC 3-associated regions were estimated to be close to 4 g/day. For ADG F , 78 ensembl_gene_ids were annotated on the declared QTL regions of OCC 3 (Table S2). One candidate gene, carbonic anhydrase 2 (CA2), was identified in the region 100.99-101.11 Mb, whereas two candidate genes, NADH:ubiquinone oxidoreductase complex assembly factor 6 (NDUFAF6) and tumor protein p53 inducible nuclear protein 1 (TP53INP1), were proposed for ADG F in the region 109.07-110.88 Mb of OCC 3 (Table 4). In OCC 5, two significantly associated regions were detected with model 1; one region comprised a single SNP in position 9.07 Mb and the other comprised two SNPs in the region 18.95-18.97 Mb. The magnitude of the strongest estimated effects for these regions was similar to those estimated on OCC 3 (between 3.5 and 5.5 g/day). In these regions, 19 ensembl_gene_ids were annotated (Table S2). Furthermore, one promising candidate gene for ADG F alpha-ketoglutarate dependent dioxygenase (FTO) was identified at 9.07 Mb in OCC 5 (Table 4). Finally, one region in OCC 21 compressing 1.29 Mb (7.17-8.46 Mb) was also associated with ADG F . In this region, 26 SNPs were found to be significantly associated with the trait. Within this region, AX-147049623, the SNP with the strongest association had an effect of 3.51 g/day. Remarkably, this SNP was located inside an intron of the Ataxin 2 (ATXN2) gene (Table S1)one of the four candidate genes (ATXN2, ACAD10, TRAFD1 and PTPN11) identified among the 71 ensembl_gene_ids annotated in this region (Table 4 and Table S2). These candidate genes contained another 10 SNPs significantly associated with ADG F (Table S1).
ADG R showed significant associations with SNPs on OCCs 9, 12, 13 and 17 (Table 3 and Fig. 3). For this trait, however, the two models declared different chromosomal regions as significantly associated with the trait. Model 1 declared a QTL region at the proximal region of OCC 13 (0.40-2.09 Mb) containing 50 significant SNPs. The estimated SNP effect having the strongest association (minimum q-value within the region) was 3.41 g/day. Two candidate genes (RC3H1 and TNFSF18), out of 37 annotated ensembl_gene_ids, were found in this region (Table 4  and Table S2). For the same trait, model 2 declared significant signals on OCCs 9 (29.66-31.00 Mb), 12 (99.88 Mb) and 17 (73.57 Mb-74.16 Mb). The SNP effects of the QTL regions on OCCs 9 and 17 were lower (approximately between 1 and 1.5 g/day) than that detected on OCC 12, for which an effect of 3.73 g/day was estimated, which was a similar magnitude to the effects of the SNPs associated with ADG F . In these regions, 74 ensembl_gene_ids were annotated (Table S2) and four candidate genes were proposed (FEZF2 and PTPRG on OCC 9, and LGALS3 and TMEM260 on OCC 17; Table 4).
The GWAS conducted with model 2 for the cage performance traits, ADFI F , FCR F and ADRFI F , declared six significantly associated regions on OCCs 5, 6, 16 and 21 (Table 3 and Fig. 4). Region 3.70-3.85 Mb on OCC 5 was declared to be associated with ADFI F and comprised 12 significant SNPs; that with the strongest association had a MAF of 0.37 and an estimated effect equal to 0.85 g/day. In this region, 20 ensembl_gene_ids were annotated (Table  S2) and two candidate genes for ADFI F were identified (CEBPA and KCTD15) ( Table 4). For FCR F , two significant signals were detected on OCCs 6 (26. . The estimated effects of the SNPs with the strongest statistical association within these regions were large -0.47 and 0.52 feed conversion units ((grams of feed/day)/(grams of growth/day)) respectively. The most significant SNP on OCC 6 had a low MAF (0.06), whereas that of the most significant SNPs on the region of OCC 16 was much higher (0.42). For FCR F , 63 ensembl_gene_ids were annotated (Table S2) and two candidate genessalt inducible kinase 1B (putative) (SIK1B) on OCC 6 and phospholipase A2 group IVA (PLA2G4A) on OCC 16 (Table 4)were retained.
Finally, for the last studied FE trait, ADRFI F , significant associations were detected in three regions of OCC 21: 3.89-4.33 and 7.16-7.70 Mb, and one single SNP at 9.21 Mb. The second region was particularly relevant as it also contained SNPs significantly associated with ADG F . The MAF of the most significant SNP within this region was 0.37, and its estimated effect was 2.16 g/day. A total of 146 ensembl_gene_ids were annotated on OCC 21 (Table S2), and the same candidate genes as those previously proposed for ADG F in this region were retrieved for ADRFI F (Table 4).

Discussion
To our knowledge, this is the first GWAS for growth and feed efficiency traits performed in a rabbit population using a dense SNP chip panel. Our study, in addition, also introduces a new modeling approach allowing study of the association of traits recorded as group averages, when not all of the individuals in the group have been genotyped. This methodology was originally proposed for gene-assisted selection when a certain percentage of the candidates have not been genotyped for the major gene of interest . Modelling the SNP allele content using animal models has also been proposed as a tool to detect low-quality SNPs within the panels (Forneris et al. 2015), an SNP being declared as erroneously genotyped when its heritability estimate is significantly different from 1. With this work, we extend the scope of application of such models to GWASs, in particular, to GWASs on group average performance traits. Previous studies (Zhang et al. 2018) have addressed the problem in the context of experiments where the limiting factor is the capability to generate individual phenotypes, but all of the individuals in the design were genotyped. In this case, it has been shown that pooling individual records to produce pool phenotypes and then explaining these pooled data by the mean genotype of the group produced considerable gains in the power of statistical tests over simple random sampling, i.e. random selection of as many individual phenotypes as pools were defined. This result could be expected as in the analyses of the pooled phenotypes all of the available genotypes are included, whereas in the study of a random sample of individual records only a subset of them are considered, and this sampling is particularly sensitive to low-frequency markers. Our study, although related to the aforementioned problem, has a completely different motivation: on the one hand, there is no individual alternative to the group average phenotype recorded, and on the other hand, the experimental limitations constrain the number of genotyped animals to only a few of those responsible for the group average phenotype. In this situation, a much smaller  amount of information is available for the analysis, and thus, lower power would be expected. We do not formally assess the efficiency of our proposed model, but its limited power seems obvious. On this regard, Sánchez et al. (2018) reported an important reduction in the capability to detect simulated QTL regions for one trait (out of three) that is considered as a group average with respect to the situation in which all traits are studied as individual phenotypes. They implemented a multitrait Bayesian procedure similar to the single-step association methods (Wang et al. 2012), which relies on derivation of SNP effects from genomic predictions using a multiple regression in which the SNP effects are treated as random factors. A possible means of validation of the results from the proposed bivariate model is to analyze those traits individually recorded with the two approaches. For the case of ADG F the results obtained are partially the same, for example regions in OCC 3 are detected with both methods. However, other regions detected for this trait and all those declared for ADG R , which is a trait with lower heritability (Piles & Sanchez 2019), are not the same across the models. One reason for this is that the datasets used by each method are different. In the case of model 1, only records from genotyped animals are considered, but with the bivariate model, records from both genotyped and non-genotyped animals are jointly considered, and the pedigree is used to predict genotypes of non-genotyped animals with records.
As stated, given the available information on the cage performance records, the expected statistical power was low. Thus, in order to allow for a certain degree of signal detection, the threshold for significance declaration was deliberately reduced to a chromosome-wide level. In this situation, we have successfully identified 17 chromosomal regions associated with the analyzed traits. To allow comparison between traits, the estimated SNP effects within the regions can be expressed relative to their estimated phenotypic variance (Table 2). To this end, we approximated the additive genetic variance associated with each Name of the most significant SNP within the region.  Table 3, by considering the SNP effect with the strongest association (minimum q-value) within the region and its frequency. The additive variances of the QTL regions in Table 3 represent 5-8, 0.5-8 and 0.5-2% of the phenotypic variance of ADG F , ADG R , and both ADFI F and ADRFI F respectively. The percentage of phenotypic variance explained by the additive genetic effect for one of the QTL regions declared for FCR F is particularly high -65% for the region on OCC 16. It could be difficult to propose a validation method for these results free from the assumptions in the model for the analysis, because FCR F is recorded as the group average. Nevertheless, a simple regression of group FCR F on group average genotype for the SNP AX-147107945 at bp 82858725 in OCC 16, the SNP with the strongest association within the region (Table 3), also showed a strong magnitude -0.20 (0.05) FCR units per unit of change on the cage average genotype (P = 3.38 × 10 −5 ). This means that the expected FCR in a cage with all of the animals heterozygous for this SNP will be 0.20 units larger than that in a cage with all of the animals homozygous of one type and 0.20 units lower than that in a cage with all of the animals homozygous of the other type. Twenty candidate genes, located in 12 QTL regions, have been proposed to explain the phenotypic variation of the traits under study; this proposition was done based on their biological functions. It is worth mentioning the FTO gene, annotated on OCC 5 for ADG F , which has been previously associated with growth traits in rabbits (Zhang et al. 2013;Zhang et al. 2014). Furthermore, NDUFAF6, which was also annotated for ADG F in a region of OCC 3, has recently been described as a candidate gene for growth-related traits in pigs (Ji et al. 2019).
It is also relevant to highlight ATXN2, acyl-CoA dehydrogenase family member 10 (ACAD10), TRAF-type zinc finger domain containing 1 (TRAFD1) and protein tyrosine phosphatase non-receptor type 11 (PTPN11) genes as they map in a region of OCC 21 with a pleiotropic effect for both ADG F and ADRFI F and have 11 significant SNPs located in their introns. Identifying pleiotropic regions for both ADG F and ADRFI F could be considered an unexpected result as ADRFI F is a trait obtained after the phenotypic correction of FI by growth and metabolic weight. However, this phenotypic correction does not grant a null genetic correlation, and in fact, in the population under study it has been reported that the genetic correlation between ADG F and ADRFI F was 0.58 (Piles & Sanchez 2019). Piles & Sanchez (2019) showed that growth recorded in animals fed under restriction (ADG R ) is a trait genetically different from growth recorded in animals fed ad libitum. Our results could be said to support this as for ADG R we have declared chromosomal regions different from those declared for ADG F . Nonetheless, this could be also a simple consequence of our reduced statistical power. In these regions, candidate genes associated with behavioral patterns (FEZF2, PTPRG and LGALS3) or involved in immunity and/or lipid metabolism (RC3H1, TNFSF18 and TMEM260) were identified. Finally, it is worth highlighting the CCAAT enhancer-binding protein alpha (CEBPA) gene, annotated on OCC 5 for ADFI F , which has recently been identified as an upstream regulator of several differentially expressed genes down-regulated in adipose tissue of high-feed-efficiency pigs (Horodyska et al. 2019).
In spite of our loose significance threshold setting, we feel relatively confident of having adequately controlled the rate of false-positive signals that we have declared. In support of our results, we have identified some candidate genes that

Conclusions
We have proposed a number of QTL regions linked to the observed variation of the studied traits using a complex statistical model for fitting cage FE and feed intake, jointly with individual genotypes. To our knowledge, this is the first time this type of statistical model has been used within the framework of GWAS studies. The information content on cage average performances is quite limited, thus we have reduced the threshold for significance declaration to a chromosome-wide level. In spite of this loose significance threshold definition, the declared QTL seem to harbor genes that can clearly be regarded as functional candidates for the traits of interest. Our results seem to support the idea that the growth of animals fed on restriction is under a different genetic control that that of animals fed ad libitum as we have identified different QTL regions for both traits. It is remarkable that genes related to behavioral patterns have been proposed as candidates for ADG R . Regarding FE, some of the QTL regions that we declared to harbor candidate genes which are involved in lipid and energy metabolism have a pleiotropic effect for both ADG F and ADRFI F : In spite of these promising results, further functional research is warranted to validate these genes. Overall, our results lay an important foundation for future studies to unravel the underlying genetic bases driving growth and FE regulation in rabbits.