**Impact of Goji Berries (***Lycium barbarum***) Supplementation on the Energy Homeostasis of Rabbit Does: Uni- and Multivariate Approach**

**Laura Menchetti 1,2, Giulio Curone 3, Egon Andoni 4, Olimpia Barbato 2,\*, Alessandro Troisi 5, Bernard Fioretti 6, Angela Polisca 2, Michela Codini 7, Claudio Canali 2, Daniele Vigo <sup>3</sup> and Gabriele Brecchia 3,\***


Received: 12 October 2020; Accepted: 28 October 2020; Published: 30 October 2020

**Simple Summary:** The energy balance during the reproduction cycle is a problematic issue for livestock species because it has consequences not only on animal welfare but also on the profitability of the farm. The adoption of new nutritional strategies could improve both of these aspects. In the present study, the supplementation with goji berries was proposed and evaluated on the rabbit, which is both a livestock animal and a useful animal model. Goji berry is the fruit of *Lycium barbarum* that is a natural resource made up of several compounds with biological activities and their consumption could be beneficial for the health and the general well-being of humans and animals. Its effect on several hormones and metabolites involved on energy balance of rabbit doe were evaluated by using both uni- and multivariate approach. Our finding, in addition to describing the intricate relationships between body conditions, hormones and metabolites during pregnancy and lactation, suggested that the supplementation with goji berry in the rabbit diet at low percentage could improve some aspects of energy metabolism and, in particular, doe's insulin sensitivity. Conversely, the intake of high doses of goji raises concerns due to the risk of excessive fattening and worsening of insulin resistance.

**Abstract:** This study examined the effects of goji berries dietary supplementation on the energetic metabolism of doe. Thirty days before artificial insemination, 75 New Zealand White does were assigned to three different diets: commercial standard diet (C) and supplemented with 1% (LG) and 3% (HG) of goji berries, respectively. Body conditions, hormones and metabolites were monitored until weaning. Body weight and BCS were higher in HG than C (*p* < 0.05). LG showed lower T3/T4 ratio and cortisol concentrations (*p* < 0.05) and tended to have lower indices of insulin resistances (*p* < 0.1) than HG. Compared to control, leptin was higher in HG at AI (*p* < 0.01) and in LG during lactation (*p* < 0.05). Two principal components were extracted by multivariate analysis describing the relationships between (1) non-esterified fatty acids, insulin and glucose levels, and (2) body conditions and leptin metabolism. The first component highlighted the energy deficit and the insulin resistance of the does during pregnancy and lactation. The second one showed that leptin, body weight and Body Condition Score (BCS) enhance as levels of goji berries in the diet increase. Thus, the effects of goji supplementation are dose-dependent: an improvement on energy metabolism was achieved with a low-dose while the highest dose could determine excessive fattening and insulin resistance in does.

**Keywords:** goji berries; rabbit; insulin resistance; leptin; non-esterified fatty acids; pregnancy; lactation; body condition score; principal component analysis

#### **1. Introduction**

Recently, research has turned towards natural products, rich in compounds with high biological activity that have favourable effects on health with low production expenses and reduced side effects for the prevention and treatment of various pathologies. In this contest, goji berry, the fruit of *Lycium barbarum* L., consumed in several Asian countries for a long time as a traditional tonic food and natural herbal remedy, has attracted a lot of attention also in Western countries in the last decades [1,2]. In fact, there are a lot of evidences showing that the berry has numerous potential beneficial effects for the general well-being of the individual and for the prevention and treatment of numerous pathologies [1–4]. Beside several biological active compounds such as carotenoids, vitamins (riboflavin, thiamin and ascorbic acid) and flavonoids, goji berry is primarily rich in polysaccharides [1] which are responsible for the main beneficial pharmacological effects of the fruit both in vitro [5–8] and in vivo in various laboratory animal species [9,10] and in clinical trials in humans [11,12]. The effects of the goji berries are mainly studied in laboratory animals such as mice and rats [9,10,13] and only a few trials were conducted using the rabbit [14–16] although it is considered a useful experimental animal model [17–22]. Moreover, only a limited number of researches evaluated the effects of goji berry on the reproductive and productive performance, other than on the quality of meat, in livestock animals, rabbits included [15,16,23–25].

During pregnancy, the energy demands increase to favor the growth of fetuses and prepare the mother to the subsequent lactation. For these reasons, some changes in the energy homeostasis are necessary. In fact, during this physiological status, the pregnant animal enhances the feed intake and the mobilization of the body reserves as an adaption mechanisms to the modification in energetic request [26–28]. Different hormones (insulin, leptin, T3, T4 and cortisol) and metabolites such as glucose and non-esterified fatty acids (NEFA) are implicated in maintaining energy homeostasis during pregnancy in different species [29–32]. To our knowledge, no studies have evaluated the impact of the goji berries enriched diet on the hormonal control of the energetic homeostasis during pregnancy of the rabbit does.

Therefore, the aims of this research were to evaluate the effects of the integration of the rabbit diet with goji berries at two different concentrations, 1% and 3%, on the body conditions as well as on the levels of several metabolic hormones and metabolites involved in energy homeostasis in pregnant and lactating rabbits does. First, the patterns of these parameters in the three groups was assessed individually using a univariate approach. Subsequently, a multivariate approach was used to identify the main hormonal and metabolic profiles and the effects of goji berries inclusion in the feed on these frameworks.

#### **2. Materials and Methods**

The study was performed at the facilities of the Faculty of Veterinary Medicine, of the Agricultural University of Tirana, Albania. The experimental protocol was in agreement with the national rules on the use of livestock animal for experimental and other scientific purposes, their handling, protection, and welfare. Every effort has been made to reduce animal discomfit and to use only the number of animals sufficient to produce valid results.

#### *2.1. Animals and Experimental Design*

Seventy-five New Zealand White nulliparous does were maintained in single cage in controlled environmental conditions: temperature ranged from +18 to +23 ◦C, relative humidity being from 60% to 75%, and the lighting schedule of 16 L:8 D. Rabbit does, on the base of a random chose, were assigned to three groups depending by the diet administered (*n* = 25/group): Control (C), 1% goji (LG) and 3% goji (HG). Does were fed with commercial feed (control) or the same feed integrated with 1% and 3% of goji berries in LG and HG, respectively (Table 1 [16]). The animals were undergoing artificial insemination (AI) after 30 days of nutritional adaptation to the experimental dietary regime. During this period, they received 150 g/d of feed, while after AI, the rabbits had ad libitum access to feed. Water was always administered ad libitum. An injection of 0.8 μg of synthetic GnRH (Receptal, Hoechst-Roussel Vet, Milan, Italy) just before AI was used to induce the ovulation [33]. AI was executed with 0.5 mL of diluted fresh semen. Pregnancy was determined by abdominal palpation 10 days after AI. Eleven does for each group were monitored during pregnancy and lactation until weaning. Lactation was controlled, the does had access to the nest once a day for ten minutes until day 18 post-partum. The nest was open on the 19th day of lactation. Weaning occurred at 35 days of age separating the little rabbits from their mothers.

**Table 1.** Formulation and chemical composition (as fed) of control (C) and experimental diets supplemented with goji berries. LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries.


\* Per kg diet: vitamin A 11,000 IU; vitamin D3 2000 IU; vitamin B1 2.5 mg; vitamin B2 4 mg; vitamin B6 1.25 mg; vitamin B12 0.01 mg; alpha-tocopherol acetate 50 mg; biotine 0.06 mg; vitamin K 2.5 mg; niacin 15 mg; folic acid 0.30 mg; D-pantothenic acid 10 mg; choline 600 mg; Mn 60 mg; Fe 50 mg; Zn 15 mg; I 0.5 mg; Co 0.5 mg. \*\* Estimated by Maertens et al. [34].

#### *2.2. Body Conditions*

During the experimental period, the feed intake was registered daily while body weight (BW) and body condition score (BCS) were recorded weekly, between 7:30 and 9:00 AM, from time 0 (before the administration of the three experimental diets) until day 35 post-partum (weaning). Weights and BCS were also measured on the day of the IA. The aggregated BCS (from 0 to 4) was obtained by summing the score (0–2) estimated for the rump and the loin [28].

#### *2.3. Hormone and Metabolite Assays*

Blood samples were collected to evaluate the metabolic hormones (insulin, leptin, T3, T4 and cortisol) and metabolites (glucose and NEFA) concentrations. These were withdrawn at the time 0 (basal), at AI, post-partum, 20th day of lactation (top lactation), and at the post-weaning. The samples were extracted from the central ear vein into tubes containing EDTA, and instantly centrifuged at 3000× *g* for 20 min; subsequently, plasma was frozen and stored until evaluated for hormones and metabolites. Radioimmunological procedures (RIA) were used to evaluate plasma cortisol, insulin, leptin, triiodothyronine (T3) and thyroxine (T4) concentrations as previously reported [26,27,33]. Leptin concentrations were determined by using the multispecies leptin kit based on a double antibody RIA method (Linco Research Inc., St. Charles, MO, USA). The intra- and inter-assay coefficients of variations were 3.4% and 8.7%, respectively while the limit of sensitivity was 1.0 ng/mL.

A porcine insulin RIA kit was used to quantify the plasma levels of insulin making use of double antibody/PEG technique (Linco Research Inc., Saint Charles, MO, USA). A purified recombinant human insulin was used to produce the standard curve as well as labelled antigen while the antibody against the porcine insulin (antiserum) was produced in the guinea pig. The minimum level of insulin detected by the kit was 2 μU/mL. The coefficients of variations were 6.8% for the intra-assay and 9.2% for the interassay.

Thyroid hormones (T3 and T4) were assayed using the procedure reported by the manufacturer on the brochure of the kit (Immunotech, Prague, Czech Republic). The sensitiveness of the analysis was 0.26 nmol/L for T3 and 10.63 nmol/L for T4. The coefficients of variations of T3 were 6.3% for the intraassay and 7.7%, inter-assay; whereas, for T4 were 3.29% for the intraassay and 7.53% for the inter-assay.

CORT kit was used to determine the cortisol plasma concentrations (Immunotech, Prague, Czech Republic). The lower level of cortisol detected by the kit was 2.5 nM and intra-assay and inter-assay coefficients of variations were 5.8% and 9.2%, respectively.

The metabolites glucose and NEFA were assayed by colorimetric methods according to Menchetti et al. [31], and García-García et al. [35] respectively. The NEFA plasma concentrations were evaluated using a colorimetric assay from Wako that use two enzymatic reactions (NEFA-C, Wako Chemicals GmbH, Neuss, Germany), and based on the ability of NEFA to acylate coenzyme A in the presence of CoA synthetase.

A glucose Infinity kit from Sigma (Sigma Diagnostic Inc., St. Louis, MO, USA), that take advantage of the glucose oxidase method, was used to assess the glucose plasma levels.

Finally, the homeostasis model assessment for insulin resistance (HOMA) was used as index of insulin-resistance and it was calculated using a formula previously reported [27,31]: [insulin concentration × (glucose concentration/18)]/22.5. Insulin resistance (low insulin sensibility) is evident when the HOMA-IR is high while low HOMA-IR values are expression of high insulin sensitivity.

#### *2.4. Statistical Considerations and Data Analysis*

Statistical analyses were performed with SPSS Statistics version 25 (IBM, SPSS Inc., Chicago, IL, USA). We defined *p* ≤ 0.05 as significant and a *p* value between 0.1 and 0.05 as a trend.

First, data were analysed by a univariate approach using the Linear Mixed Model. The models evaluated the effects of the Group (three levels: C, LG or HG), Time (as days, depending on the dependent variable) and Group × Time interaction. The Time was included as a repeated factor. The Sidak method was used for multiple comparisons. We verified assumptions and outliers by using diagnostic graphics. Results were expressed as estimated marginal means ± standard error (SE) but raw data were presented in figures.

Then, the principal component analysis (PCA) was used to describe relationships between body condition, hormones, and metabolites and convert the numerous variables into a few profiles describing the hormonal and metabolic framework. The very low or very high correlations between variables were checked using correlation matrices in order to identify suitable variables for the PCA [36–38]. Moreover, the Determinant of the correlation matrix and the Bartlett's test of Sphericity were also used. The sampling adequacy was tested using a Kaiser–Meyer–Olkin (KMO) value >0.6 as index of adequate factorability. The number of components (PCs) to retain was chosen according to the Kaiser's criterion (eigenvalues > 1) and Varimax rotation was applied. The interpretation of PCs was mainly based on items having loadings greater than |0.4|. Then, the scores were calculated for each PC using regression techniques to create two new variables (PC1 and PC2) as linear combinations of the indices of body conditions, hormones, and metabolites.

Finally, the effect of percentage of goji berries inclusion on these PCs was evaluated using Generalized Linear Models (GLMs). Normal and Identity were set as the probability distribution and the link function, respectively, while the Time was included as Within-Subject Effects. In the GLMs, the two PCs were evaluated separately as dependent variables and the percentage of goji berries inclusion was included as predictor [16]. In addition, the number of days of supplementation was included as covariate. Results were expressed as b coefficient with SE and *p* value from Wald Chi-square statistics [16,39].

The G\*Power® software was used to calculate the achieved power [40]. An F test was chosen for this power analysis by setting α = 0.05 and the effect size = 0.25 (medium effect size) [41]. For the 3 experimental groups (C, LG or HG) and the five repeated measures (basal, AI, post-partum, 20th day of lactation, and post-weaning), a power of 90.5% was reached using a total sample size of 33 animals.

#### **3. Results**

#### *3.1. Univariable Approach*

#### 3.1.1. Body Conditions

BW was influenced by time (F = 9.59; *p* < 0.001) and group (F = 5.43; *p* = 0.005). In supplemented rabbits, the BW increase was early, starting from the first week of pregnancy in LG (4044 ± 128 g; *p* < 0.05) and from AI in HG (4054 ± 88 g; *p* < 0.01) groups. After pregnancy, it returned at baseline values after one week of lactation in C (3969 ± 151 g) and LG (4000 ± 128 g) while, in the HG group, it remained elevated up to 20 d of lactation (4220 ± 139 g) (Figure 1).

The aggregate BCS showed significant effects for time (F = 9.51; *p* < 0.001), group (F = 4.19; *p* = 0.016) and their interaction (F = 1.82; *p* = 0.007). Marginal means were 1.1 ± 0.1, 1.3 ± 0.1, and 1.5 ± 0.1 in C, LG, and HG, respectively. HG group showed higher BCS values than C starting from week 3 of supplementation, (*p* < 0.001; Figure 2). From mid-pregnancy onwards, the differences were no longer significant, perhaps also due to the high variability of the data.

Food intake did not differ between groups during pregnancy (144 ± 2 g/d, 146 ± 2 g/d, and 148 ± 3 g/d for C, LG, and HG, respectively; F = 0.78; *p* = 0.459) while it was greater in the control group than supplemented rabbits during lactation (320 ± 2 g/d, 301 ± 2 g/d, and 287 ± 3 g/d for C, LG, and HG, respectively; F = 42.00; *p* < 0.001).

**Figure 1.** Changes in the body weight (BW) during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

**Figure 2.** Changes in the body condition score (BCS) during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

#### 3.1.2. Hormone and Metabolite

Overall, insulin concentrations reduced from baseline (6.4 ± 0.4 μU/mL) to post-partum (3.7 ± 0.6 μU/mL; F = 6.71; *p* < 0.001). They increased during lactation (6.5 ± 0.6) but, at weaning, lower concentrations than baseline were found (3.7 ± 0.8 μU/mL; *p* = 0.003). With regards the group effect, only a trend was found at post-partum, when values of LG tended to be lower than HG (*p* = 0.081; Figure 3).

**Figure 3.** Changes in insulin concentrations during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

The HG group showed the highest leptin concentrations at AI (2.9 ± 0.1 ng/mL; *p* < 0.01). Subsequently the concentrations decreased (*p* < 0.05) and returned to baseline values after 20 d of lactation (2.5 ± 0.3 ng/mL). Constant values over time were found in C group (marginal mean: 2.3 ± 0.1 ng/mL) while LG showed greater leptin concentrations during lactation (2.7 ± 0.2 ng/mL) compared to baseline (2.2 ± 0.1 ng/mL; *p* = 0.013; Figure 4).

**Figure 4.** Changes in leptin concentrations during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

A progressive reduction until postpartum was observed in T3 concentrations (F = 9.17; *p* < 0.001). This reduction was more marked in the LG group which showed lower marginal means (1.9 ± 0.1 nmol/L and 1.6 ± 0.1 nmol/L in C and LG groups; respectively; *p* < 0.05) and lower postpartum values (1.7 ± 0.2 nmol/L and 1.0 ± 0.2 in C and LG groups, respectively; *p* < 0.05) than the control. After weaning, baseline values were found for both C and LG while HG showed lower T3 concentration (2.4 ± 0.2 nmol/L and 1.3 ± 0.3 nmol/L at time 0 and post weaning, respectively; *p* < 0.01; Figure 5).

**Figure 5.** Changes in T3 concentrations during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

Only changes over time (F = 5.97; *p* < 0.001) were found for T4 which decreased until post-partum (from 51.8 ± 1.9 nmol/L at T0 to 38.8 ± 2.6 nmol/L at post-partum; *p* < 0.001); subsequently, it increased (46.5 ± 0.2 nmol/L at day 20 of lactation) but its values on the last time, (42.9 ± 3.8 nmol/L) were lower than in the first sampling time (*p* < 0.05). No group-related differences were significant (F = 0.75; *p* = 0.475; Figure 6).

**Figure 6.** Changes in T4 concentrations during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

Changes over time in T3/T4 ratio were only found for the Control group (*p* = 0.043) which showed significantly higher values than the LG group at the postpartum (0.034 ± 0.003 and 0.029 ± 0.005 for C and LG, respectively; *p* = 0.010; Figure 7).

**Figure 7.** Changes in T3/T4 ratio during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

Cortisol concentrations progressively decreased in all groups, from 334.3 ± 7.4 nmol/L at baseline to 258 ± 14.4 nmol/L at post-weaning (F = 12.63; *p* < 0.001). The LG group (272.5 ± 7.0 nmol/L; *p* < 0.01) had, on average, the lowest concentrations of cortisol but no multiple comparisons reached statistical significance (Figure 8).

**Figure 8.** Changes in cortisol concentrations during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

No change over time was observed in the control group while a reduction in glucose concentrations compared to baseline (5.7 ± 0.4 mmol/L) was observed post-partum (4.2 ± 0.4 mmol/L; *p* < 0.001) and during lactation (4.3 ± 0.4 mmol/L; *p* < 0.001) in the LG group. Moreover, marginal means of LG group (5.5 ± 0.2 mmol/L) was lower than C (5.8 ± 0.2 mmol/L; *p* < 0.05; Figure 9).

**Figure 9.** Changes in glucose concentrations during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

Conversely, NEFA showed significant changes over time (F = 30.67; *p* < 0.001) with no difference between groups (F = 1.81; *p* = 0.167). In particular, greater values on their concentration were found at post-partum (0.33 ± 0.02 mmol/L; *p* < 0.001) and post-weaning (0.37 ± 0.03; *p* < 0.001) compared to baseline (0.14 ± 0.02 mmol/L; Figure 10).

**Figure 10.** Changes in non-esterified fatty acids (NEFA) concentrations during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

Overall, a reduction in HOMA was found at post-partum (0.05 ± 0.01; *p* < 0.01) and post-weaning (0.05 ± 0.01; *p* < 0.01) compared to baseline values (0.09 ± 0.01). This index tended to be lower in LG (0.06 ± 0.01) than HG (0.08 ± 0.01; *p* = 0.063) and, in particular, lower values were found in the post-partum (*p* = 0.072) and at day 20 of lactation (*p* = 0.030; Figure 11).

**Figure 11.** Changes in homeostasis model assessment for insulin resistance index (HOMA) during nutrition adaptation and productive cycle of rabbit does. C = control diet; LG = diet supplemented with 1% of goji berries; HG = diet supplemented with 3% of goji berries. AI = artificial insemination.

#### *3.2. Multivariable Approach*

According to the correlation matrix, all the variables indicating body condition and hormonal and metabolic profiles were included in the PCA except for HOMA (Table 2). The Bartlett's Test of Sphericity (*p* < 0.001) indicated that this dataset was suitable for a data reduction technique and the KMO (=0.70) confirmed the sampling adequacy of the PCA. The PCA extracted two PCs, which together account for 46.1% of the variance (Table 2). The PC1 was a bipolar component including NEFA with negative loading and insulin, cortisol, glucose and thyroid hormones with positive loadings. In the PC2, the highest loadings were found for leptin and body condition indices, BW and BCS (Table 2 and Figure 12).


**Table 2.** Loadings of factors extracted with the principal component analysis.

Factor loadings with an absolute value greater than 0.4 are in bold. NEFA = non-esterified fatty acids; T3 = triiodothyronine; T4 = thyroxine; BW = body weight; BCS = Body Condition Score.

**Figure 12.** Component plot in rotate space (Varimax rotation). The position of the variables in the Cartesian graph indicates that PC1 (x-axis) was positively associated with insulin, cortisol, glucose, and thyroid hormones while negatively with NEFA levels. The PC2 (y-axis) was positively associated with body weight (BW), leptin, and BCS.

The two PCs were then analysed by regression techniques to evaluate how the two hormonal-metabolic profiles change as the concentration of supplementation increases (Table 3). Independently from the duration of supplementation, for every 1% increase in the goji concentration, the PC2 score increased by 0.275 (*p* < 0.001). This result indicated that the effects of goji on body conditions and leptin levels are dose-dependent: as the percentage of goji supplementation increases, the BW, the leptin, and the BCS of does tend to increase (Figure 13). In contrast, there was no linear relationship between goji berries inclusion in the diet and PC1 (b = 0.019, *p* > 0.1), which included NEFA, insulin, cortisol, glucose, and thyroid hormones.


**Table 3.** Results of regression analyses including percentage of goji berries inclusion in feed as independent variables and principal component as outcome. Days of supplementation were included as covariate.

**Figure 13.** Scatter plot showing the relationship between the percentage of goji berries inclusion in feed (independent variable, x-axis) and PC2 (dependent variable, y-axis). The PC2 was positively associated with body weight, leptin, and BCS. The parameters of regression line were estimated by a generalized linear model including the days of supplementation as covariate.

Both components changed linearly over time (*p* < 0.01). The signs of b indicated a progressive reduction of PC1 (b = −0.017, *p* < 0.001) while, with lower slope and significance, an increase over time for PC2 (b = 0.007, *p* < 0.01).

#### **4. Discussion**

Energy homeostasis is finely regulated by many factors which, in a synergistic way, have the task of maintaining a constant body weight. Maintaining this balance is never easy and, in fact, there are frequent alterations that lead to pathological conditions such as obesity and diabetes. During pregnancy, energy homeostasis is even more complicated because the body must also ensure the growth of the foetus and prepare for the next lactation [27,31,32]. The energy balance during a reproduction cycle is one of the most problematic issues for rabbit breeding because it has consequences not only on the animal welfare but also on the profitability of the farm [28,29]. The adoption of new nutritional strategies for the rabbit doe could, therefore, improve both of these aspects. In the present study, a product which is receiving growing interest as a nutraceutical [1,2] and that has shown promising results for the rabbit's growth performance [15,25] was proposed: the goji berry. Two percentages of supplementation were evaluated in order to find the optimal concentration for goji inclusion in the food of rabbit doe.

The effects of goji supplementation were evaluated on different hormones and metabolites responsible for the maintenance of the energy homeostasis: leptin, insulin, cortisol, thyroid hormones, NEFA, and glucose. All of them are involved in the changes of energy intake, efficiency, and expenditure but the success is only assured if their actions are integrated and strictly coordinated. For this reason, in addition to the description of the individual patterns of these hormones and metabolites, an assessment of the overall picture during the doe's reproduction cycle was proposed by using multivariate analysis techniques.

In the group fed with standard diet, lower insulin, T3, and cortisol values were found at the end of the production cycle compared to the pre-pregnancy. The T4 remarkably decreased at the post-partum meanwhile the T3/T4 ratio increased. Later, an increase in T4 concentrations was observed but, after weaning, its values were lower than in the first sampling time. Conversely, significant increases were observed for NEFA concentrations with peaks at the post-partum and post-weaning. The concentrations of leptin and glucose did not undergo great changes in the timing of collection. Finally, HOMA insulin sensitivity index showed a fluctuating trend but, at the end of the observation period, its values had decreased.

These hormones and metabolites have recently been evaluated by Menchetti et al. [31] in pregnant and pseudopregnant rabbits. However, this previous study focused on changes during the 32 days post insemination, the sampling was weekly, and no observations were done after the birth. By using weekly sampling, these authors highlighted short-term changes during pregnancy, especially in insulin and leptin. These changes suggested a condition of insulin and leptin resistance that characterizes the pregnancy of many mammals including women [30,32,42–45]. In the present study, the changes were assessed at longer intervals but the concentrations at early and late pregnancy were comparable to those of Menchetti et al. [31]. Moreover, the longer-term evaluations of the present study allowed to highlight the physiological conditions of the lactating rabbits. In particular, the post-partum increases in NEFA and insulin resistance confirmed the state of negative energy balance during early lactation. In this phase, in fact, there is the big gap between the energy the energy required for maintenance and milk production and that taken with the diet. Several factors can contribute to this severe energy deficit of the does, such as its prolificacy, the quantity and characteristics of the milk, and the intense reproductive rhythm [29,46].

This state is characterized by an intense mobilization of the body reserves which increase the NEFA concentrations. At the same time, the increase in insulin was associated with moderate concentrations of glucose triggering a peripheral insulin resistance status [47] which could contribute to raise NEFA levels. Insulin is, indeed, an anti-lipolytic hormone and the insulin resistance can favor lipolysis and fatty acid availability.

The condition of insulin resistance was exacerbated when the highest goji supplementation was included in the feed (3%, HG group). A similar condition has been found in lactating dairy cows consuming excessive concentrate [47,48]. The insulin resistance develops as a homeorhetic adaptation to ensure a sufficient glucose supply for the lactating mammary gland [49] but it has been negatively associated with milk production, reproduction performance, and animal health [50,51]. The insulin resistance observed in the present study could therefore explain the low milk yields reported by Menchetti et al. [15] in rabbits receiving 3% of goji supplementation. Goji berries contain high amounts of polysaccharides and polyphenols [1]. Although most of the studies reported that these compounds exert beneficial effects on metabolic diseases [4,8,52], their excessive consumption during pregnancy and lactation is debated. For example, Caimari et al. [53] demonstrated that supplementations with a grape seed extract (rich in polyphenols) in pregnant and lactating rats induced insulin and adiponectin resistance. On the other hand, other authors report several positive effects of the administration of polyphenols in lactation: they seem to increase the immunoglobulins content in colostrum of sows and improve the preweaning survivability [54], to protect bovine mammary epithelial cells from oxidative stress [55] and also to exert positive effects on the adult offspring's metabolism [56,57].

Thus, our findings suggest that the quantity and source of energy during lactation affect milk yield in rabbit doe by modulating body mobilization and insulin sensitivity. Moreover, they raise concerns about the excessive consumption of polyphenol-enriched foods during lactation although dose, source, species, and duration of polyphenols administration should be evaluated from time to time.

Menchetti et al. [15] also reported greater BWs in rabbits fed with the 3% of goji. In the present study, a greater BCS was found in the same group of does. This could suggest that high-doses of Goji supplementation also leads to excessive fattening which, in rabbit does, is negatively associated with reproductive performance [58].

The effects of a lower supplementation diverged: the does supplemented with 1% goji showed lower T3 at post-partum as well as lower cortisol levels during lactation; their leptin concentrations were lower at AI and higher during lactation while insulin, glucose and HOMA values suggested a higher insulin sensitivity during lactation. This hormonal framework would seem to favour the milk yield as high productions were recorded by Menchetti et al. [15] with the low-dose goji supplementation. For example, the low levels of T3 at the beginning of lactation could decrease thermogenesis diverting

nutrients to the mammary glands. As mentioned above, the positive association between insulin sensitivity and lactation has been already showed in ruminants and humans [59]. These findings suggest that the effects of goji berries on insulin resistance are dose-dependent and a moderate consumption could improve insulin sensitivity during the reproductive cycle.

The increase in leptin concentrations during lactation could explain the differences in food intake between groups as this is the main function of leptin, as well as the most studied [60]. Some authors reported hypoleptinemia during lactation and have hypothesized that leptin resistance may also continue after giving birth [61,62]. However, our results do not seem to support these theses. On the other hand, maternal leptin circulating during lactation may be important for the neurological development of the new-borns because it is secreted into milk and it is taken up during suckling [63]. As observed in several species including humans, the leptin ingested trough suckling could have a role not only in the changes during the perinatal period but also in the metabolic programming of adulthood phenotype [64,65]. Finally, Koch et al. [66] demonstrated that rabbit mammary tissue expresses leptin mRNA and protein during lactation suggesting that it could be also involved in the mammary gland development.

In the present study, all the variables described so far were then included in the multivariate analysis to identify the hormonal and metabolic profiles characterizing rabbit pregnancy and lactation. Two dimensions were extracted.

The conditions of insulin resistance and energy deficit were well described by the first component extracted by the principal component analysis. In this component, insulin and glucose had the same sign (both positive loadings) and this could be explained by the reduction in insulin sensitivity: insulin rises as the glucose increases but, under conditions of insulin resistance, cells fail to respond normally and high blood glucose levels are still maintained. In the first dimension, instead, the insulin was opposed to NEFA levels. The anti-lipolytic effect of insulin could explain the opposite signs of these items. Indeed, the known effect of systemic insulinization is a decline in circulating NEFA concentrations [67]. However, rabbits' pregnancy and lactation are characterized by elevated NEFA levels. To explain the negative association between NEFA and insulin levels in this context, it can also be hypothesized that chronically elevated NEFA levels reduce the insulin-secretory capacity in rabbit does, as already proposed for high-yielding dairy cows [51].

Overall, the first dimension progressively reduced, suggesting that insulin, cortisol, glucose, and thyroid hormones decreased during the productive cycle of the doe while, in the meantime, the NEFAs increased. Similar profiles are described in lactating cows [47].

The second dimension extracted by the principal component analysis included leptin, a hormone secreted mainly by adipose tissue [60]. Therefore, its positive associations with body weight and BCS found in this component are not surprising. The positive loadings of glucose and insulin in this dimension also confirm that they could be factors implicated in leptin regulation [60]. Of major interest is the positive relationship between this second component and goji supplementation indicating that leptin, BW and BCS enhance as the percentage of goji berries inclusion in the diet increases. In agreement with the univariate results, this finding suggests that there is a risk of excessive fattening in rabbit does supplemented with goji and that, probably, a high percentage of inclusion ultimately determines detrimental reproductive effects.

It is worthwhile noting that 3% of goji supplementation had a different impact in growing rabbits: it guaranteed good growth performance and an improvement in the meat quality [15,25]. The use of goji in rabbit diet was also positively perceived by consumers who probably appreciate the use of nutraceutical products also in animal nutrition [25]. Overall, the goji inclusion in the feed could be a promising nutritional strategy from both a marketing and animal welfare point of view. However, the dosage of goji berries should be differentiated according to the category of animals and physiological status.

#### **5. Conclusions**

The rabbit doe is subjected to a severe energy deficit during reproductive cycle so that several physiological phenomena, such as the increase in NEFA and insulin resistance, are comparable to those of lactating cows. The quantity and source of energy during pregnancy and lactation modulate body mobilization, insulin sensitivity and leptin levels that, in turn, could affect milk yield. This balance is so delicate that the effects of goji supplementation diverged according to the dose: a moderate consumption seems improve insulin resistance while a high-dose could make animals excessively fat and reduce insulin sensitivity. This work has proved that goji berry could play a role in rabbit nutrition, but the formulation should be carefully modulated according to the physiological state.

**Author Contributions:** Conceptualization, L.M. and G.B.; methodology, L.M., G.C., D.V. and G.B.; software, L.M.; formal analysis, L.M.; investigation, M.C., O.B., C.C. and E.A.; resources, E.A. and O.B.; data curation, L.M.; writing—original draft preparation, L.M. and G.B.; writing—review and editing, L.M., G.C., D.V. and G.B.; visualization, L.M.; supervision, G.B.; project administration, D.V. and G.B.; funding acquisition, O.B., B.F., A.P. and A.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study has been financially supported by Regione Umbria (Italy) (Grant PSR 2007/2013; no. 44750050831).

**Acknowledgments:** The authors wish to thank Giovanni Migni for his excellent technical assistance and Margherita Angelucci for the English revision.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Relationship between Body Chemical Composition and Reproductive Traits in Rabbit Does**

**Meriem Taghouti 1, Javier García 2, Miguel A. Ibáñez 3, Raúl E. Macchiavelli <sup>4</sup> and Nuria Nicodemus 2,\***


**Simple Summary:** At the beginning of the productive life of rabbit does, there must be a balance between ensuring at least a minimal degree of bodily development to guarantee a successful reproductive life, and the minimization of the unproductive rearing period, but nowadays there is no clear recommendation about the optimal moment for the first artificial insemination (AI). A better body condition at the first AI (higher body protein, fat and energy), that indicates a higher degree of maturity of the rabbit doe, did not influence fertility at the first AI (that is usually very high), but improved it at the second AI (that is usually lower than the first one). The percentage of kits born alive at the first and at the second AI also were positively influenced by the body protein content at the first AI. We can conclude that the degree of maturity at the first AI is a key point to optimize the does reproductive success, with body fat and body protein content being relevant factors.

**Abstract:** The relationship among live weight, chemical body composition and energy content (at artificial insemination (AI) and three days before parturition), estimated by bioelectrical impedance with fertility rates and the percentage of kits born alive, was studied during the first three AI. The first AI was conducted at 16 weeks of age in 137 rabbit does that weighted 3.91 ± 0.46 kg. Their body chemical composition was 17.4 ± 0.50%, 16.1 ± 2.6%, 1067 ± 219 kJ/100 g body weight, for protein, fat and energy, respectively. An increase in body protein, fat and energy content at the first AI did not affect fertility at the first AI but improved it at the second AI (*p* ≤ 0.030). Moreover, an increase in body fat and energy content at the second AI improved fertility at the second AI (*p* ≤ 0.001). Fertility at the third AI was positively influenced by body protein at the third AI and the increase in body protein and fat between the second parturition and the third AI (*p* ≤ 0.030). The percentage of kits born alive at the first and at the second AI improved with the increase in body protein at the first AI (*p* ≤ 0.040). In conclusion, a minimal body protein and fat content is required at the first AI to optimize the reproductive performance in young does.

**Keywords:** body composition; fertility; kits born alive; rabbit does

#### **1. Introduction**

In the last four decades, rabbit production underwent a noticeable change from a traditional and familiar organization to industrial and intensive systems. Consequently, genetic selection programs and new breeding management systems were established, improving the production of the new hybrid lines used. This has led to an increase in

**Citation:** Taghouti, M.; García, J.; Ibáñez, M.A.; Macchiavelli, R.E.; Nicodemus, N. Relationship between Body Chemical Composition and Reproductive Traits in Rabbit Does. *Animals* **2021**, *11*, 2299. https:// doi.org/10.3390/ani11082299

Academic Editors: Rosa María García-García and Maria Arias Alvarez

Received: 29 May 2021 Accepted: 29 July 2021 Published: 4 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

female nutritional requirements [1,2], health problems [3–7] and welfare necessity [8–11]. Nutritional strategies in reproductive females need to be global and consider both the shortterm productive factors (litter size, milk production, or fertility) as well as the long-term factors (body condition or health). Therefore, the expected improvement in nutritional management should be based on an accurate analysis of the requirements of the doe, its evolution during successive reproductive cycles and the identification of crucial moments in the life of rabbit does to optimize productivity and longevity.

Two of the key points of the reproductive success of rabbit does are their birth weight and maturity at the first artificial insemination (AI). There is an optimal threshold for birth weight (>57 g) to optimize the initial reproductive performance, which is associated with an increase in the live weight and fat reserves at the first AI [12–14]. The latter are the traits used to define the maturity of rabbit does at the first AI, and both are also related to the nutritional rearing strategy and the time of AI [13,15]. However, there is no clear recommendation indicating the weight and/or body condition of the nulliparous rabbit doe at the onset of its reproductive life [16]. In the current management systems, does are inseminated at a fixed age with minor or no consideration of their weight and chemical body composition. As in other species, there must exist a balance between ensuring at least the minimal degree of body development needed to guarantee a successful reproductive life, and minimizing the unproductive rearing period. Accordingly, the study of the body chemical composition of rabbit does seems to be a useful tool not only to improve feeding, but also for general rabbit doe management [17]. The final aim is to extend the lifespan of rabbit does, which is limited by the relatively high early mortality and culling rate in intensive production systems [5].

A new non-destructive technique to estimate rabbit doe body chemical composition (moisture, ash, protein, fat, and energy content) based on the bioelectrical impedance measurement, live weight and physiological status of rabbit does, was developed by Pereda [18]. This method is easier and cheaper than TOBEC [19], and in both methods, the variations of gut contents are included in the error term. Furthermore, it allows the prediction of total fat and energy content (not only the perirenal fat content, as it is with the case with the ultra-sound technique [20]), as well as the body protein, which is not usually estimated with other methodologies to evaluate body condition.

The aim of this work was to establish the relationship between chemical body composition (at AI and parturition), determined using bioelectrical impedance, and both fertility and the percentage of kits born alive during the first three inseminations, and to identify the most important moments to record the body chemical composition.

#### **2. Materials and Methods**

#### *2.1. Animal Husbandry and Management*

Data of the estimated body chemical composition, the fertility and percentage of kits born alive recorded in two different farms in 2010 were used in order to obtain a wide variation in chemical body composition. In farm A, 106 crossbred (New Zealand White × Californian) rabbit does from the UPV hybrid line (genetic line selected for prolificacy in Polytechnic University of Valencia in Spain using the crossline A × Line V) were used. After the first parturition, does were submitted to three different breeding systems, defined by the parturition–AI and parturition–weaning intervals (4/32, 11/35 and 14/42, Table 1) to obtain a wide variability of body condition. In farm B, 37 crossbred (New Zealand White × Californian) rabbit does from the Hyplus hybrid line (prolific hybrid maternal line selected by Hypharm in France) were used and they were inseminated 11 d after parturition and litter weaned at 35 d of age. Rabbit does were inseminated for the first time between 16 and 18 weeks of age. Rearing period of rabbit does was not controlled and they were fed ad libitum after the first insemination. Non pregnant-non lactating rabbit does were restricted to around 150 g/d. Two commercial diets were used, one in farm A containing 17.5% crude protein, 32.0% neutral detergent fibre and 9.95 MJ digestible energy/kg (based on dehydrated alfalfa, wheat bran and barley), and another in

farm B containing 16.5% crude protein, 30% neutral detergent fibre, and 10.8 MJ digestible energy/kg (based on dehydrated alfalfa, sugar beet pulp, sunflower meal, rye and wheat). All rabbit does were submitted to a cycle of 16 h light/8 h dark, and heating, cooling and forced ventilation systems allowed the building's temperature to be maintained between 18 and 24 ◦C. Animals were handled according to the principles for the care of animals in experimentation published by the Spanish Royal Decree 53/2013 [21] and favourably assessed retrospectively by the Ethics Committee of the Polytechnic University of Madrid.


**Table 1.** Description of the animals and breeding systems used.

The relationship between body chemical composition and fertility was studied during the first three cycles as the number of replicates decreased from the third parturition onwards (Table 2). The percentage of kits born alive was only studied at the first and second parturition (N = 121 and 83, respectively) due to the reduction in replicates in the third parturition. The range of variation in chemical body composition of these does at the moments of insemination and parturition is shown in Table 2.

**Table 2.** Estimated chemical body composition <sup>1</sup> of nulliparous and primiparous rabbit does during the first three reproductive cycles at artificial insemination and parturition.


<sup>1</sup> LW: Live body weight (g); Moist. (moisture); Prot. (protein); fat and ash: % LW; Ener. (body energy): kJ/100 g body weight. <sup>2</sup> SD: Standard deviation; Min.: Minimal; Max.: Maximal.

> Seminal doses with more than 20 million spermatozoa in 0.5 mL of a commercial diluent (Magapor S.L.) were made using a pool of fresh heterospermic semen from bucks selected for growth performance. In order to synchronize the oestrus in the second, third and fourth insemination, 48 h before insemination, the does were injected with 25 IU of eCG (Equine Chorionic Gonadotropin, Segiran, Lab. Ovejero, León) [22]. On the day of insemination, the does received an intramuscular injection of 10 μg of buserelin Suprefact®

(Hoechst Marison Roussel, S.A., Madrid). Buserelin is a Gonadotropin-releasing hormone agonist (GnRH agonist), which is used to induce ovulation in rabbit does [23].

#### *2.2. Experimental Procedures and Data Recording*

In order to determine the chemical body composition of the does, a bioelectrical impedance analysis (BIA) technique was used [18,24]. Measurements were taken with the four-terminal analyzer Quantum II (Model BIA-101, RJL Systems, Detroit, MI, USA). This device generates an alternating current of 425 μA of intensity at a frequency of 50 kHz. It is provided with 2 black electrodes to conduct the electrical current through the doe's body, and 2 red electrodes to register the resistance and reactance resulting from the passing current. When it is used with rabbit does, a needle (Terumo, 21G × 1 1/2, 0.8 × 40 mm, nr 2) is inserted in the extremity of each electrode. The needles must pass through the skin of the rabbit doe at four reference points along the loin (two near the scapula and other in the distal part of the loin). Impedance is defined by the equation: Impedance = (Resistance<sup>2</sup> + Reactance2) 1/2. These three parameters, in addition to the physiological status, live body weight and parturition order were used further to estimate the body composition of the rabbit does. Rabbit does were weighed and their body composition estimated on the days of artificial insemination and parturition (three days before parturition, both pregnant and non-pregnant does) during the first three reproductive cycles. The regression equations developed and validated by Nicodemus et al. [24] and Pereda [18] enabled the prediction of the doe's body content in moisture, protein, fat and ash expressed as percentages or as g/100 g of body weight. Furthermore, energetic body content expressed in MJ or in kJ/100 g of body weight could be also estimated.

#### *2.3. Data Treatment and Statistical Analysis*

In order to establish correlations among parameters of chemical body composition (moisture, fat, protein, ash and energy) at the moments of insemination and parturition during the first three parities, we used Pearson's correlation coefficient (CORR procedure of the SAS system) (SAS Inst., Cary, NC, USA). The GENMOD procedure of the SAS system [25] was used to study the relationship between chemical body composition, fertility during the first three parturitions and the percentage of kits born alive in the two first parturitions, as the logistic regression [26,27] is an adequate tool for modeling proportions [28]. In the generalized linear models used, the link function employed was 'logit' (Equation (1)) and we considered that both fertility and the percentage of kits born alive were proportions arising from a binomial distribution. In equation 1, 'p' is the mean of the proportion of fertility and kits born alive.

$$\text{Logit}\,(\mathbf{p}) = \ln\left(\mathbf{p}/1-\mathbf{p}\right); \mathbf{p}\in(0,1)\tag{1}$$

Our first aim was to determine whether the variations in fertility and percentage of kits born alive (dependent variables) were significantly affected by the body chemical composition and its variation between insemination and parturitions (independent covariates). The model also included, as fixed effects, the farm and the breeding system (except for nulliparous does). Afterwards, only significant covariates were retained for the interpretation of coefficient estimates (β) and the linear predictor (η = x'β). In order to display the covariates' effects, three levels (low, medium and high), representative of the range of variation of each covariate, were fixed for each significant covariate (body constituent) based on the data used in this investigation. To predict the means of fertility and percentage of kits born alive expected for each level (the mean or predicted values, 'p') the inverse logit function was used (Equation (2)).

$$\mathbf{p} = \mathbf{e}^{\eta/} (\mathbf{1} + \mathbf{e}^{\eta}) \tag{2}$$

Among the criteria of the goodness of fit, we used the likelihood ratio chi-square, 'QL', also called 'deviance', to establish whether models fitted appropriately [29]. The inverse

function of the logit transformation was also used to calculate a confidence interval for the expected percentages of fertility and kits born alive for each level [29]. In fact, in the SAS output and for each estimator calculated, there are the values of lower (Lη) and upper (Uη) limits for its 95% confidence interval. In this case, the method consists in transforming the linear values of the limits to the log-link scale using Equation (3).

$$\text{Upper limit } \mathbf{p} = \mathbf{e}^{\text{U}\eta} / (1 + \mathbf{e}^{\text{U}\eta}) ; \text{Lower limit } \mathbf{p} = \mathbf{e}^{\text{L}\eta} / (1 + \mathbf{e}^{\text{L}\eta}) \tag{3}$$

To calculate a standard error for the mean of predicted values, 'p', for each level of the covariates retained, we used the Delta method [30], which involves the standard error of the coefficient estimate (SEη) as well as the predicted value, as explained in Equation (4).

$$\text{SEp} \approx \text{p}(1-\text{p})\,\text{SEp} \tag{4}$$

The results were transformed from the logit scale.

#### **3. Results**

#### *3.1. Correlation among Does' Corporal Constituents during the First Three Parturitions*

The initial live weight and chemical composition of rabbit does had a wide range of variation, which was more important for the fat than for the protein content (16 vs. 3%, respectively, for coefficient of variation (Table 2)). Rabbit does showed an important degree of growth between the first and the second AI (3906 vs. 4187 g, weight of does at 1st and 2nd AI, *p* < 0.05 (Table 2)). Live body weight was positively correlated (*p* < 0.001) with fat (r = 0.43 to 0.83) and energy body content (r = 0.43 to 0.71), and negatively with moisture (r = −0.35 to −0.72) (expressed on % and MJ/100 g body weight, respectively) from the first AI to the third parturition (Table 3). In this period, body fat was inversely and closely correlated with moisture (r = −0.95 to −0.99), and positively with body energy (r = 0.79 to 0.98), while body protein was positively correlated with ash content (r = 0.30 to 0.83). Regarding moisture, a negative correlation was observed with energy (r = −0.43 to −0.71).

**Table 3.** Correlation among estimated chemical body constituents <sup>1</sup> during the first three reproductive cycles at insemination and parturition 2.


1. LW: Live body weight (g); Moist. (moisture); Prot. (protein); fat and ash are expressed in % LW; Ener. (energy): KJ/100 g body weight. <sup>2</sup> \*\*\*: *p* < 0.001; \*\*: *p* < 0.01; \*: *p* < 0.05; †: 0.05 < *p* < 0.10.

> Protein content at the first AI was not correlated with protein content at any other moment of the cycle (Table 4). However, protein content was positively correlated at parturition and the subsequent AI (both at the first and second parturition; *p* < 0.001). In contrast, the body fat content was mainly positively correlated in the period between the first AI and the third parturition, and the same was observed for body energy content (*p* < 0.05 (Table 4)).


**Table 4.** Correlation among estimated protein, fat and energy body contents at artificial insemination (AI) and parturition during the first three cycles 1.

<sup>1</sup> N = 137, 132, 96, 96 and 84 rabbit does for 1st partum, 2nd AI, 2nd partum, 3rd AI and 3rd partum, respectively. <sup>2</sup> NS: Not significant (*p* > 0.05). \*\*\*: *p* < 0.001; \*\*: *p* < 0.01; \*: *p* < 0.05; †: 0.05 < *p* < 0.10.

#### *3.2. Relationship between Body Condition and Fertility during the First Three Parturitions: Nulliparous, Primiparous and Multiparous Rabbit Does*

The fertility rate registered in the first parturition was 93.4%. There were no differences between data from the two farms used for the analysis (*p* > 0.05). Neither live body weight, nor its chemical composition at the first AI, affected the fertility rate at this moment. In fact, there were no differences in the initial body composition and live weight at the first AI between pregnant and non-pregnant does (Table 5), although these results have to be taken with caution due to the small number of non-pregnant rabbits.

**Table 5.** Estimated chemical body composition of nulliparous pregnant and non-pregnant rabbit does at the first artificial insemination (AI) 1.


<sup>1</sup> N = 128 and 9 for pregnant and non-pregnant does, respectively. <sup>2</sup> Measured 3 d before parturition. RMSE: root mean square error.

As expected, in the second parturition, fertility was lower than in the first parturition (56.2 vs. 93.4%, respectively; *p* < 0.05). Breeding system and farm did not have any effect on fertility in the second parturition. Recorded fertility means were 54.0%, 50.8% and 65.8% for the studied breeding systems R4, R11, and R14, respectively. Live body weight at the moments of first AI, first parturition, and second AI did not affect the fertility rate in the second parturition. Fertility in the second AI was related to chemical body composition at the first AI, in contrast to that observed at the first parturition. Body protein content (*p* = 0.007), fat (*p* = 0.030), and energy (*p* < 0.001) at the first AI positively affected the fertility rate in the second AI (Table 6). The relationship between these parameters and fertility

is lineal in the logit scale. Consequently, the model was fitted using three fixed levels of body protein, fat, energy, and protein/energy ratio. These levels were chosen to cover the data range used for the analysis. Fat and energy were positively correlated at the first AI (r = 0.79, *p* < 0.001), and with their values at the second AI (r = 0.33 and 0.74; *p* < 0.001 (Tables 3 and 4)), while at the first AI protein was negatively correlated with fat but did not have any correlation with energy, and was not related to protein at the second AI. The three levels determined for each constituent and used in the model are shown in Table 6. The increase in body protein from 16 to 18% and body fat from 10 to 20% increased fertility from 19.9 to 72.0% and from 31.3 to 69.5%, respectively. The results showed that does with higher energy content (1400 vs. 700 kJ/100 g) underwent improvements in fertility from 19.1 to 83.6%.


**Table 6.** Effect of estimated body composition at the first insemination on fertility in the second insemination 1.

<sup>1</sup> N = 132. <sup>2</sup> SE: Standard error. <sup>3</sup> LL and UL: lower and upper limits for a confidence interval (95%). <sup>4</sup> QL: Deviance of the models.

A higher body protein/energy ratio at the moment of the first AI negatively influenced the fertility rate in the second AI (*p* < 0.001 (Table 6)). This ratio reflects again the importance of the energetic content, since high ratios (30 g/MJ) corresponded to low corporal energy (800 kJ/100 g) rather than to high body protein content. In fact, the latter showed lower variability (2.9 vs. 20.5%, coefficient of variation of body protein and energy at the first AI, respectively (Table 2)).

Chemical body composition at the second AI was also related to fertility in the second AI, as detailed in Table 7. Fat and energy contents at the second AI were positively related to fertility at this moment (*p* < 0.001) in a similar way to their relation at the first AI. These two variables were positively correlated both at the first and at the second AI (*p* < 0.001 (Table 3)), and were also positively correlated with fat and energy content at the first AI (*p* < 0.001; Table 4), although their coefficient of variation increased at the second AI (Table 2).

Body protein at the second AI did not exert any effect on fertility at the second AI, in contrast with its observed effect at the first AI. It was positively correlated with fat and energy content at the second AI. Protein contents at the first and at the second AI were not correlated (Table 4), although their average values and variability were similar (Table 2), and were lower than body protein values observed from second parturition onwards. A higher body protein/energy ratio at the second AI again negatively affected the fertility rate at the second AI, reflecting the effect of energy content at the second AI.


**Table 7.** Effect of estimated body composition at the second AI on fertility in the second artificial insemination 1.

<sup>1</sup> N = 132. <sup>2</sup> SE: Standard error. <sup>3</sup> LL and UL: lower and upper limits for a confidence interval (95%). <sup>4</sup> QL: Deviance of the models.

Moreover, to study the factors that influenced the fertility at the second AI, multiple linear regressions were fitted considering two independent variables: live weight or body composition at the first AI combined with their respective variations between the first and second AI (Table 8). Live body weight at the first AI and weight gain during the interval between the first AI and the first parturition as well as during the interval between the first and second AI were positively related to fertility at the second parturition.

**Table 8.** Effect of live weight at the first artificial insemination (AI) and its increment from the first insemination to the first parturition (model 1) or from the first to the second inseminations (model 2) on fertility at the second parturition.


<sup>1</sup> QL: Deviance of the models. <sup>2</sup> Increment-I: increment of body weight between the first AI and the first parturition expressed as a percentage of body weight at the first AI. <sup>3</sup> Increment-II: increment of body weight between the first and second AI expressed as a percentage of body weight at the first AI. <sup>4</sup> SE: standard error. Model 1: *p*-weight = 0.010; *p*-increment1 = 0.006. Model 2: *p*-weight = 0.010; *p*-increment2 = 0.001.

The fertility rate recorded at the third AI was 73.9%. No effect of breeding system or farm was found. Recorded fertility means were 56.2.0%, 90.3% and 75.8% for the studied breeding systems R4, R11, and R14, respectively, showing differences in the first two values (*p* = 0.007). The fertility rate recorded at the third AI was not affected by body composition (fat and energy) at the moments of the first and second AI. Unlike the fertility at the first and second AI, we registered a positive effect of protein content on fertility at the third AI (*p* = 0.030 (Table 9)). Body protein content at the third AI and second parturition were positively correlated (Table 4), and none of them were correlated with body fat or energy content (Table 3). Moreover, we studied the effect of the changes in chemical components between consecutive AI and parturitions. A positive effect of body protein gain between the first and third AI (*p* = 0.020), and between the second parturition and third AI (*p* < 0.001 (Table 9)), was observed. Similarly, body fat gain between the second parturition and third AI was also positively related to fertility rate at the third AI (*p* = 0.020 (Table 9)).


**Table 9.** Effect of estimated body protein content at the third artificial insemination (AI), its increments and fat increment on fertility at the third artificial insemination 1.

<sup>1</sup> N = 96. <sup>2</sup> SE: standard error. <sup>3</sup> LL and UL: lower and upper limits for a confidence interval (95%). <sup>4</sup> Protein content at 3rd AI. <sup>5</sup> Protein increment between the moments of the first and third AI expressed as a percentage of the initial protein content.6 Protein increment between the moments of third AI and second parturition expressed as a percentage of protein content at the second parturition. <sup>7</sup> Fat increment between the moments of the third AI and second parturition expressed as a percentage of fat content at the second parturition.

The effect of the previous reproductive success (or not) is reflected in the strong effects of both protein and fat gain before the third AI on the fertility at the third AI. In fact, pregnant does at the third AI showed lower fertility at the second AI, which allowed them a better body reserve recovery (especially protein) at the third AI compared to non-pregnant does (that showed higher fertility in the previous AI (Table 10)). This is confirmed by the positive relationship between fertility rates at the third and first AI (*p* = 0.050), and the negative relationship between fertility at the third and second AI (*p* = 0.050). We observed that does that gave birth at the second parturition had lower chances of becoming pregnant at the third AI (*p* = 0.043; Table 10).

**Table 10.** Effect of the change of live body weight and body composition, between the second parturition and the third artificial insemination (AI), and fertility at the second AI on the reproductive success at the third artificial insemination 1.


<sup>1</sup> N = 96. <sup>2</sup> expressed as a percentage of the content at 2nd parturition. <sup>3</sup> RMSE: root mean square error. <sup>4</sup> NS: non-significant.

#### *3.3. Relationship between Body Condition and Percentage of Kits Born Alive during the First Two Parturitions: Nulliparous and Primiparous Does*

The percentage of kits born alive out of the total born in the first parturition was 93.4%, and the number of kits born alive per doe was 8.00 ± 2.97. Among chemical body constituents, body protein (*p* = 0.040) and energy contents (*p* = 0.010) at the first AI increased the percentage of kits born alive (Table 11). Both variables were not correlated at the first AI (Table 3). The relationship between these parameters and the percentage of kits born alive is linear in the logit scale. Consequently, the model with logit link was made using three fixed levels of body protein and energy contents. These levels were chosen to cover the data range used for the analysis. When protein content at the first AI increased from 16 to 18%, the percentage of kits born alive increased from 88.7% to 95.2%, respectively. Furthermore, when body energy at the first AI increased from 900 to 1300 kJ/100 g, the percentage of kits born alive increased from 92.4 to 95.0, respectively.

**Table 11.** Effect of body protein and energy content at the first insemination on the percentage of kits born alive at the first parturition 1.


<sup>1</sup> N = 121. <sup>2</sup> SE: standard error. <sup>3</sup> LL and UL: lower and upper limits for a confidence interval (95%). <sup>4</sup> Deviance of each model.

Body fat content was negatively correlated with body protein and positively correlated with energy content at the first AI (r = −0.21 and 0.79, respectively; *p* < 0.05 (Table 3)), and no effect on the percentage of kits born alive was observed. Live weight at the first AI was also correlated with body protein, energy and fat contents (r = −0.46, 0.65 and 0.83, respectively; *p* < 0.001 (Table 3)) but was not related to the percentage of kits born alive. Furthermore, a higher percentage of kits born alive was observed in farm B compared to farm A (98.1 and 88.7%, respectively; *p* < 0.001), but the number of kits born alive per doe was not different between the two farms.

The percentage of kits born alive recorded in the second parturition was 87.5% and the number of kits born alive per doe was 10.3 ± 4.46. This rate is lower than that observed in the first parturition (93.4%), but the number of kits born alive increased (8.00 ± 2.96 kits born alive/doe in the first parturition). No effect of the breeding system or farm on the percentage of kits born alive at the second parturition was reported. Weight and body composition at the first AI were not related to this trait. However, a higher body protein at the second AI (*p* < 0.001) was associated with an increase in the percentage of kits born at the second parturition (Table 12).

**Table 12.** Effect of body protein content at the second insemination on percentage of kits born alive at the second parturition 1.


<sup>1</sup> N = 83. <sup>2</sup> SE: standard error. <sup>3</sup> LL and UL: lower and upper limits for a confidence interval (95%). <sup>4</sup> Deviance of the model.

#### **4. Discussion**

A general problem observed in rabbit does is their low fertility rate in the second parturition [31]. It is usually explained by the energy deficiency observed before the second insemination [16]. In this period, it might be difficult for the rabbit doe to meet the requirements for both pregnancy and growth due to the limited feed intake [1,32–35]. However, when the energy supply was increased, it was dedicated to milk production with no limitation of reserve mobilization [36,37]. In this way, Pascual [13] hypothesized that this negative energetic balance would be a natural adaptation to optimize evolutionary success. In this context, it is interesting to study the relationship between the factors related to body chemical composition and their influence on fertility and kit survival at birth, in order to identify the threshold to be met by rabbit does at the beginning of their reproductive life.

The evolution of the live weight and chemical composition of rabbit does from their first AI onwards indicated that they were still growing when inseminated the first two times, which agreed with recent data [38]. Live body weight, body energy and fat were closely and positively correlated, which was similar to the correlation between perirenal fat thickness and body energy content reported previously [39,40]. In contrast, body protein had a minor or no correlation with the latter traits, but a negative one with live weight from the second parturition onwards. These results agree with rabbit doe maturation in this period, which would depend on maturity at the first AI and on reproductive success. Once maturity is reached (or nearly reached), the changes in body weight might be mainly associated with fat mobilization and/or deposition. This would agree with the moderated and positive correlation between live weight and body condition score [41].

The differences in chemical body composition and live weight between the first AI and the first parturition seemed to be related to the success in the first insemination, which was not influenced by body condition or live weight. This effect was reported by other authors [2,34,42]. It may be explained by the specific situation of non-pregnant does, which would use the entire intake for body protein and fat accretion, and accordingly, energy accretion, towards the completion of the final step to reach their maturity, where the fat deposition is much higher than the protein deposition (24 vs. 6% increment, respectively). Meanwhile, pregnant does have to supply gestation requirements that are especially important during the last 10 days of gestation, and which can impair not only fat content [17,43,44], but also protein balance [45], with respect to non-pregnant does. It must be taken into account that rabbit does inseminated later show a higher feed intake capacity [46] and lower growth requirements. Furthermore, rabbit does reduce feed intake in the days before parturition, contributing to the impairment of their nutrient balance [43].

The impairment of fertility in the second insemination reflects the specific situation of primiparous rabbit does, which suffer a negative energetic balance during their first pregnancy and lactation, compared to non-pregnant does, which seems to negatively influence reproductive performance and especially fertility [2,32,34], although this could be considered a natural adaptation, as commented before [13]. Live body weight at the moments of first AI, first parturition and second AI did not affect the fertility rate in the second parturition. This agrees with the findings of Rommers et al. [46], who did not observe any effect of body weight at the first AI on fertility rate in the first two parturitions. Nevertheless, the combination of a high live weight and high weight gain between the first two AI was also related to better fertility. Rabbit does that lose weight between the first two AI, regardless their initial body weight, were unable to present an acceptable reproductive performance in the second cycle. In this period, live weight was positively correlated with body fat and energy content, but not with protein, as does are finishing their protein accretion. The relationship between the fertility in the second AI and the chemical body composition at the first AI confirms the importance of the rearing management of rabbit does. Diets used [2,14,42,47] and the time of first insemination [15,48–50] influenced body composition, which consequently affected fertility. Therefore, at the end of the rearing period, reproductive does should reach an optimal body condition (minimal body protein, fat and energy content), assuring an adequate feed intake and body development, which

enable high fertility rates during first parturitions [16]. In this sense, Pascual et al. [13] stated that body data at the first AI are a sign of doe soma and might be related to its productive potential. In this context, the supplementation of reproductive sows with certain daily amounts of amino acids enabled an adequate retention of nitrogen that led to an acceptable reproductive performance [51]. In this study, the threshold in the body composition at the first AI to avoid a sharp reduction in fertility in the second insemination might be set at 18% protein and 20% fat, but few does met this condition (12 and 7%, respectively). Other studies where the time of first insemination was delayed (to 18.4 or 19.5 weeks of age) rendered an increase in the body protein at that moment, as well as higher fertility values (83–87%) that were also associated with a higher fat content [38], although this was not always the case [52]. When the latter two studies were considered together, the does that were successful in the first five AI were lighter and had less body fat than the average (although their mean and range were similar to the current study), and the same body protein (although it was in the upper threshold that was previously mentioned: 17.9%) [53], suggesting the potential relevance of body protein at the beginning of the productive life. Another difference between these extraordinary does and the average population was their higher fat mobilization between the second AI and the first weaning (this was recovered between weaning and the third AI) [52]. These results partially differed from those of Theilgaard et al. [54], who indicated that there was no positive effect of perirenal fat at the first AI on reproductive life. In fact, a higher risk of culling was associated with high fat mobilization, although does that were too lean also seemed to increase their risk of culling. Similarly, Castellini et al. [50], using perirenal fat, found that does that were too fat and too lean (at AI) showed the poorest fertility. Recent results confirmed the negative effect of fatness at the first AI on the risk of being culled and litter size [14]. The disagreement among these studies and the current one might be related to the different fatness range (probably the absence of does that are too fat at the first AI in our work: maximal fat content at the first AI: 22.1%; Table 2), which might be associated with the time of the first AI (in the latter studies, nulliparous does were inseminated later than in the current one). Body fat and energy at the second AI was also related to fertility in the second AI, which might also reflect the observed positive influence of initial body condition on fertility at the second AI.

Quevedo et al. [55] suggested that the success of AI 11 days after parturition was conditioned by the rabbit doe's body condition at parturition rather than at insemination. This fact was not observed in this work, probably due to the fact that in our work, we measured body condition at parturition three days before birth (to avoid disturbing the doe) instead of after birth. This prevented the recording of an important proportion of fat mobilization, which was described by Savietto et al. [17]. Once confirmed, the BIA measurement did not alter the doe immediately after parturition (unpublished results); subsequent studies recorded it just after parturition of the doe [38,52].

The reproductive success at the beginning of reproductive life also influenced the fertility of the third AI. In fact, rabbits that are reproductively successful during the first two parturitions were more vulnerable, and consequently, their reproductive performance was impaired in the third AI if they did not have the opportunity to recover. In this sense, Castellini et al. [50,56] proposed the delaying of the second AI after weaning in order to allow rabbit does to recover properly from the first gestation lactation and continue their growth. Otherwise, reproductive success at the beginning of the reproductive life of rabbit does, when they cannot recover their body reserves, might worsen the rabbit doe's productivity and shorten its life span [13].

The absence of effects of the breeding system and farm may be explained by the synchronization of rabbit does at the moment of AI. Rebollar et al. [57] did not register a rhythm effect on fertility in the second parturition when using controlled lactation as the does' synchronization tool. However, in experiments without any synchronization method [56,58], it was concluded that the reproductive rhythm was related to fertility rate. Anyway, it must be stressed that the current study was not designed specifically to study the effect of breeding system on fertility rate. Besides, an influence of the breeding systems

on the body condition could not be ruled out and further studies would be required to figure out the nature of their relationship.

There were also a positive influence of body protein and energy at the first AI on the percentage of kits born alive at the first parturition. Rommers et al. [15,59] also related a higher protein content at the first insemination, and lower fat content, with a trend towards increasing numbers of kits born alive and percentages of kits born alive at the first parturition. They observed that restricted nulliparous rabbit does inseminated at 17.5 weeks, compared with those fed ad libitum and inseminated at 14.5 weeks, showed higher protein and lower fat content with a similar live weight, and tended to increase the number of kits born alive and the percentage of kits born alive at the first parturition. In primiparous rabbit does, a better body condition (higher body protein, lipid, and energy contents) was also related with changes in metabolic signals (increase in serum protein and leptin concentrations) that might influence ovarian follicle and gamete quality and might be associated with an improved reproductive outcome [60]. Similarly, in sows, ovarian activity and oocyte quality were influenced by the body protein content [61,62]. Another explanation may involve body protein being related to fetal survival, as the increase in litter size has been related to a higher fetal survival, independently of ovulation rate [63].

Live weight and body fat content were also positively correlated with body energy content at the first AI, but negatively with body protein, and had no effect on the percentage of kits born alive. In contrast, Rommers et al. [46] reported that heavier females at the first AI (>4.0 kg) decreased the percentage of kits born alive, although this was combined with an improvement in the litter size at the first parturition. They related it to the development of the reproductive apparatus (larger uterine horns and more *corpora lutea* in the ovaries). Moreover, these results did not agree with those of Quevedo et al. [44], where rabbit does with higher perirenal fat thickness at 3 months of age tended to increase their percentage of kits born alive at the first parturition.

The reduction in the percentage of kits born alive in the second parturition, and the increase in the number of kits born alive, agreed with the results reported by Rommers et al. [46], who also observed, at the second parturition, a higher number of kits born alive and a lower percentage of kits born alive with respect to the first parturition. The different percentage of kits born alive observed in the two farms might be due to the different hybrids used and/or the different environmental management conditions in each farm. No effect of the breeding system or farm on percentage of kits born alive at the second parturition was detected. Weight and body composition at the first AI were not related to this trait. However, higher body protein at the second AI increased the percentage of kits born alive. This result is similar to that recorded at the first AI (on the percentage of kits born alive at the first parturition) and again suggests a positive role of nitrogen content on conception success and fetus viability.

#### **5. Conclusions**

An adequate body chemical composition at the first AI (around 18% protein and 20% fat) allowed a better fertility at the second AI. However, the consecutive reproductive success at the first and second AI did not allow rabbit does to recover body reserves and impaired fertility at the third AI. Body composition also affected the percentage of kits born alive at the first AI and at the second AI that increased with body protein content.

Consequently, rearing management (e.g., time of first AI) is key to avoiding low fertility rates in primiparous rabbit does. Furthermore, when reproductive success is reached, rabbit does may require alternative management strategies to recover their body condition. Finally, determining body composition at moments of AI may be an adequate tool to anticipate the reproductive success possibilities of rabbit does. Further studies considering different ages at the first AI and breeding systems are warranted to confirm these conclusions and to clarify the relationship between breeding systems and body conditions. **Author Contributions:** N.N., M.T. and J.G. contributed to the concept and design of the experiment. M.T., N.N., and J.G. conducted the data processing, with the help of R.E.M. and M.A.I. in the statistical analysis and modelling. All authors contributed to the interpretation of the results. M.T., J.G. and N.N. wrote the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was financed by the project MINECO-FEDER (AGL2008-00627) and the grant from IAMZ-CIHEAM obtained by Myriam Taghouti to develop her MSc. thesis.

**Institutional Review Board Statement:** The study was conducted according to the guidelines of the Declaration of Helsinki, and approved by the Institutional Review Board (or Ethics Committee) of Universidad Politécnica de Madrid (approved retrospectively on 28 July 2021).

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to confidentiality requirements agreed with one of the farms.

**Acknowledgments:** This work was promoted by Jhonny Demey (RIP, formerly of the Universidad Central de Venezuela, Maracay, Venezuela), who we all have in our memories.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **Genotype Imputation to Improve the Cost-Efficiency of Genomic Selection in Rabbits**

**Enrico Mancin 1,†, Bolívar Samuel Sosa-Madrid 2,\*,†, Agustín Blasco <sup>2</sup> and Noelia Ibáñez-Escriche 2,\***


**Simple Summary:** Genotyping costs are still the major limitation for the uptake of genomic selection by the rabbit meat industry, as a large number of genetic markers are needed for improving the prediction of breeding values by genomic data. In this study, several genotyping strategies were examined through simulation scenarios to disentangle the best feasible options of implementing genomic selection in rabbit breeding programs. Most scenarios emphasized the genotyping of candidate animals with a low Single Nucleotide Polymorphism (SNP) density platform. Imputation accuracies were high for the scenarios with ancestors genotyped at high or medium SNP-densities. However, the scenario with male ancestors genotyped at high SNP-density and only dams genotyped at medium SNP-density showed the best economically feasible strategy, taking into account the tradeoff among genotyping costs, the accuracy of breeding values and response to selection. The results confirmed that by combining the imputation technique with a mindful selection of the animals to be genotyped, it is possible to achieve better performance than Best Linear Unbiased Prediction (BLUP), reducing genotyping cost at the same time.

**Abstract:** Genomic selection uses genetic marker information to predict genomic breeding values (gEBVs), and can be a suitable tool for selecting low-hereditability traits such as litter size in rabbits. However, genotyping costs in rabbits are still too high to enable genomic prediction in selective breeding programs. One method for decreasing genotyping costs is the genotype imputation, where parents are genotyped at high SNP-density (HD) and the progeny are genotyped at lower SNP-density, followed by imputation to HD. The aim of this study was to disentangle the best imputation strategies with a trade-off between genotyping costs and the accuracy of breeding values for litter size. A selection process, mimicking a commercial breeding rabbit selection program for litter size, was simulated. Two different Quantitative Trait Nucleotide (QTN) models (QTN\_5 and QTN\_44) were generated 36 times each. From these simulations, seven different scenarios (S1–S7) and a further replicate of the third scenario (S3\_A) were created. Scenarios consist of a different combination of genotyping strategies. In these scenarios, ancestors and progeny were genotyped with a mix of three different platforms, containing 200,000, 60,000, and 600 SNPs under a cost of EUR 100, 50 and 11 per animal, respectively. Imputation accuracy (IA) was measured as a Pearson's correlation between true genotype and imputed genotype, whilst the accuracy of gEBVs was the correlation between true breeding value and the estimated one. The relationships between IA, the accuracy of gEBVs, genotyping costs, and response to selection were examined under each QTN model. QTN\_44 presented better performance, according to the results of genomic prediction, but the same ranks between scenarios remained in both QTN models. The highest IA (0.99) and the accuracy of gEBVs (0.26; QTN\_44, and 0.228; QTN\_5) were observed in S1 where all ancestors were genotyped at HD and progeny at medium SNP-density (MD). Nevertheless, this was the most expensive scenario compared to the others in which the progenies were genotyped at low SNP-density (LD). Scenarios with low average costs presented low IA, particularly when female ancestors were genotyped at LD (S5) or non-genotyped (S7). The S3\_A, imputing whole-genomes, had the lowest accuracy of

**Citation:** Mancin, E.; Sosa-Madrid, B.S.; Blasco, A.; Ibáñez-Escriche, N. Genotype Imputation to Improve the Cost-Efficiency of Genomic Selection in Rabbits. *Animals* **2021**, *11*, 803. https://doi.org/10.3390/ani11030803

Academic Editors: Rosa M. García-Garcia and Maria Arias Alvarez

Received: 2 February 2021 Accepted: 5 March 2021 Published: 13 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

gEBVs (0.09), even worse than Best Linear Unbiased Prediction (BLUP). The best trade-off between genotyping costs and the accuracy of gEBVs (0.234; QTN\_44 and 0.199) was in S6, in which dams were genotyped with MD whilst grand-dams were non-genotyped. However, this relationship would depend mainly on the distribution of QTN and SNP across the genome, suggesting further studies on the characterization of the rabbit genome in the Spanish lines. In summary, genomic selection with genotype imputation is feasible in the rabbit industry, considering only genotyping strategies with suitable IA, accuracy of gEBVs, genotyping costs, and response to selection.

**Keywords:** genomic selection; imputation; litter size; rabbits; genomic simulation

#### **1. Introduction**

The rabbit industry still plays an important role throughout the agricultural sector in some European countries, such as Spain, Italy, France, Hungary, Portugal, Germany, Belgium, Poland and Malta [1,2]. In recent years, the rabbit industry is currently facing a critical period, mainly due to the increase in feeding, management costs and a constant decline in rabbit meat consumption. Hence, farmers and researchers have been looking for promising strategies to improve the current situation, one of them being the optimization of genetic selection by using genomic information. Genetic selection can improve productive and reproductive traits, such as meat characteristics, reducing feeding and management costs, which makes rabbit more appealing for consumers, and hence, will aid the rabbit industry [1]. Reproductive traits, especially litter size, are those with relevant economic weight in the rabbit industry [3]. However, the response to selection for this trait has been relatively low using traditional selection by Best Linear Unbiased Prediction (BLUP), mainly due to its low-heritability [4–6]. Selection using genomic information can be an efficient tool to guarantee a higher genetic gain for traits with low heritability and measured in only one sex such as litter size. [7,8]. Genomic selection has generally showed better accuracy for the predicted breeding value (BV) [9] due to a potentially more accurate kindship estimation between animals. This method has produced positive selection results for traits in dairy cattle [10,11], poultry [12] and pigs [13–15]. However, in rabbits, a highdensity commercial Single Nucleotide Polymorphism (SNP) platform (~200K SNPs) was not available until 2015, which has delayed genomic selection application [16]. Further, additional issues such as the small economic value of paternal rabbits, the costs of the commercial SNP platform, and the short generation interval are still limiting genomic selection as an evaluating method [17]. Strategies allowing us to diminish high genotyping costs are vital in the rabbit industry. The imputation of a low SNP-density platform using a high SNP-density platform has been carried out in other species, keeping or improving the genetic progress. The technique normally consists of genotyping ancestors (grandparents and parents) at high SNP-density in order to assign (impute) the most likely SNP allele to missing genotypes of young candidate animals genotyped at lower SNP-density [10,18,19]. This approach depends on multiple factors concerning the level of SNP density, methods, the structure and size of the reference population, the minor allele frequency of missing or untyped SNPs, the genetic architecture of traits and the particular breeding scheme of a livestock species [20–22].

The aim of this study was to disentangle the most appropriate imputation strategies for implementing genomic selection in maternal rabbit breeding schemes. Under this aim, the imputation accuracy was evaluated from low to moderate SNP-density platforms considering the cost-effectiveness of each strategy. In addition, we investigated how the imputation strategies influence the estimation of breeding values and consequently the response to selection.

#### **2. Materials and Methods**

#### *2.1. Simulation Structure*

We used stochastic simulation to create the populations for developing imputation analyses. The structure of simulations is shown in Figure 1. The simulations were performed by the AlphaSim program [23], available at https://alphagenes.roslin.ed.ac.uk/ wp/software-2/alphasimr/ (accessed on 10 June 2020). The objective was to simulate a rabbit population mimicking a single maternal line under the common rabbit breeding scheme of a small nucleus breeding. The simulated trait was litter size at birth. First, founder genomes were generated. The rabbit genome was simulated by sampling 2000 haplotype sequences for each of 20 chromosomes using the Markovian Coalescent Simulator (*MaCS*) [24]. Each chromosome was 100 cM long, as genetic distance, and included 124.43 × <sup>10</sup><sup>6</sup> base pairs, as physical distance. Chromosomes were simulated according to rabbit population history defined by the *MaCS* program using a per-site recombination rate 8.57 × <sup>10</sup>−9, a per-site mutation rate 1.74 × <sup>10</sup>−9, and an effective population size varying over time (according to "Internal rabbit" as an option of population history in *MaCS*). Later, in the first generation of the foundation of the maternal line (initial population), quantitative trait nucleotides (QTNs) were chosen randomly from the segregating sequence variants and an equal number of QTNs were assigned to each chromosome. The QTNs had additive effects sampled from a Gamma distribution with a shape of 0.60 and a scale of 0.80. These parameters were chosen after exploratory analysis evaluating various values of Gamma distribution parameters against selection response after 20 generations (analyses not shown). The heritability and residual variance were calculated relative to the additive variance in the initial population. The heritability was 0.113 and the trait genetic variance was 0.675 in the base population [5,25].

**Figure 1.** Design of simulations. A first period, generating the linkage disequilibrium (LD) across the rabbit genome using *MaCS* program. The remaining periods are carried out using *AlphaSim*: foundation of maternal line, selection for litter size with Best Linear Unbiased Prediction (BLUP), and with Genomic Best Linear Unbiased Prediction (GBLUP).

After five generations of random mating, for details as in [26], a base population was established using 138 does and 77 sires. The selection was carried out to select 70 does and 35 sires per generation, involving 20 generations with the BLUP method (see the parameters in Table S1). Litter size was assumed as a trait of sex-limited expression, hence, only does presented phenotypic records. In every generation, the 70 does from the previous generation with the highest estimated breeding values (EBVs) were selected to produce the next generation. A total of 35 males were also selected, which stand for 25 principal and 10 surrogate sires used in practice. In addition, two further generations evaluated with the Genomic Best Linear Unbiased Prediction (GBLUP) method were generated to obtain the

last population, which consisted of a progeny with 1500 does and 1500 males. The two selection methods used in these simulations (BLUP and GBLUP) are directly implemented in the *AlphaSim* program.

The high SNP-density platform (HD; 200,000 SNPs) used in genomic selection was generated by AlphaSim. Another two SNP platforms were set up according to SNPdensity size: medium SNP-density (MD; 6000 SNPs) and low SNP-density (LD; 600 SNPs) platforms. All SNPs were assigned proportionally per chromosome. Genetic markers for the MD platform were selected at random from the HD platform and, consequently, genetic markers for the LD platform were selected at random from the MD platform (see the parameters in Table S1). Regarding costs for implementing genomic selection, we only considered the costs of genotyping based on the SNP platforms in the last three generations; therefore, other indirect costs were ignored. Each platform contains 96 "BeadChips" (genotyping proof). We assumed that the cost of the HD platform is EUR 9600. For the other platforms, the costs are around EUR 4800 and EUR 1056 for MD and LD, respectively. This information came from Thermo Fisher Scientific Inc. supplied by genotyping budgets for previous genomic experiments in the Animal Breeding Group at the Institute for Animal Science and Technology in the Universitat Politècnica de València.

An underlying polygenic nature was assumed for litter size, representing the genetic architecture of complex traits [27]. We considered two QTN models, depicting: (1) a polygenic trait controlled by many QTNs with a total of 880 (44 per chromosome; QTN\_44) and (2) a polygenic trait controlled by a small number of 100 QTNs (5 per chromosome; QTN\_5). Simulated data were generated from 36 replicates for every QTN model. The results were summarized over 36 replicates within the QTN model, and presented graphically. The graphics were obtained by the R program [28].

#### *2.2. Imputation Strategies*

All imputation analyses were developed by AlphaImpute program *v1.96*. This program is available at https://alphagenes.roslin.ed.ac.uk/wp/software-2/alphaimpute/ (accessed on 15 August 2020). AlphaImpute uses a hybrid imputation algorithm in which the first step consists of a long-range phasing process, followed up haplotypes library construction, and finally, pedigree-based imputations are performed [29].

Imputation accuracy (IA) was measured by the Pearson's correlation between the imputed allele and the true genotype at untyped SNP markers. The correlation was computed one individual at a time and averaged over individuals [22]. This parameter stands for the genotype probability for reasons sketched out in [30]. Genotype yield was also computed for each imputation strategy below, as the percentage of the SNP allele calls at untyped SNP markers after the imputation process.

To assess the trade-off between IA, genomic prediction accuracy and genotyping cost, a number of hypothetical test scenarios were established, as outlined below and represented in Table 1. We analyzed the simulated data in seven sets of hypothetical scenarios with a replicate of the third scenario (S3\_A). All grand-sires and sires were genotyped at HD platforms, whilst grand-dams and dams were genotyped according to the imputation strategy. The progenies were genotyped at the LD platform except the first scenario (S1), as this scenario used MD platforms to genotype progeny and HD platforms for the genotyping of grand-dams and dams. The second scenario (S2) was like the first scenario, but LD platforms were used for progeny genotyping. The third scenario (S3) used MD platforms for the genotyping of grand-dams, and only half of the progeny was genotyped. S3\_A included the further half of the progeny, having their imputed wholegenomes. The fourth scenario (S4) used MD platforms for the genotyping of grand-dams and dams, whilst the fifth scenario (S5) used LD platforms. The sixth scenario (S6) used MD platforms for genotyping dams, but the grand-dams were non-genotyped. The seventh scenario (S7) had non-genotyped grand-dams and dams. These above-mentioned scenarios summarized all exploratory analyses concerning imputation strategies.


**Table 1.** Structure of imputation strategies and number of genotyped animals for implementing genomic selection in rabbits.

The "S" stands for scenarios. Animals of every scenario were genotyped according to different Single Nucleotide Polymorphism (SNP) density platforms. HD: high SNP-density platform (200,000 SNPs), MD: medium SNP-density platform (6000 SNPs), LD: low SNP-density platform (600 SNPs). *i-*HD: animals with imputed genotypes to high SNP-density. *i-*WG: animals with imputed whole-genome. NG: non-genotyped animals.

#### *2.3. Estimating Breeding Values and Response to Genomic Selection*

Typically, genomic selection entails a training population or reference population (with genotyped and phenotyped individuals) and an evaluated population of young candidates (with only genotyped individuals), using pseudo-phenotypes for sex-limited traits [7,31]. In this study, the rabbit reference population comprised up to 300 females (Table 2). As the reference population was small, genomic prediction was estimated using the Single-Step GBLUP (ssGBLUP) algorithm. It allows us to evaluate jointly nongenotyped and genotyped animals, combining pedigree and markers information into one matrix [32,33]. Otherwise, prediction by GBLUP hinders a greater rate of genetic progress compared to BLUP selection, with numerically small reference populations [11,34].



<sup>1</sup> Selected does are 47% of total females (150). Each female contributes to the progenies with the same proportion of males and females (1:1). <sup>2</sup> Number of males: females with genotypes in each generation. The number varies according to the imputation strategy given standard *BLUPF90* parameters of quality control.

In each simulation, the accuracy of genomic breeding values (gEBVs) was estimated using the imputed genotypes of the eight scenarios described above. In addition, the EBVs of candidate animals were also estimated using only pedigree information (BLUP scenarios). The model was the same for BLUP and ssGBLUP:

$$y = \mathbf{1}\mu + \mathbf{Z}\mathfrak{a} + \mathfrak{e} \tag{1}$$

where *y* is the vector of phenotypes, *μ* is the population's mean, *a* is the vector of additive genetic effects of animals, *e* is the vector of residuals, and *Z* is the incidence matrix for additive genetic effects. Residual effects are sampled from distribution *N*(**0**,*Iσ*<sup>2</sup> *<sup>e</sup>* ). In BLUP, random additive genetic effects are sampled from distribution *N*(**0**,*Aσ*<sup>2</sup> *<sup>a</sup>* ); *σ*<sup>2</sup> *<sup>a</sup>* is the genetic additive variance and *A* is the identity by descent (IBD) relationship matrix constructed from pedigree information. In ssGBLUP, additive genetic effects are sampled from distribution

with *<sup>a</sup>*∼*N*(**0**,*Hσ*<sup>2</sup> *<sup>a</sup>* ). *H* matrix can be interpreted as the (co)variances of multivariate normal distributions of additive effects of both genotyped and non-genotyped animals [32,33]. In ssGBLUP, the inverse of the (co)variance structure of random effect was replaced by *H*<sup>−</sup>1, described as:

$$H^{-1} = A^{-1} + \begin{pmatrix} 0 & 0 \\ 0 & \mathbf{G}^{-1} - A^{-1}\_{22} \end{pmatrix} \tag{2}$$

where *A*−<sup>1</sup> is the inverse of pedigree matrix and *A*−<sup>1</sup> <sup>22</sup> is the inverse of the sub-covariance structure containing only genotyped animals. *G*−<sup>1</sup> is the inverse of the matrix built as described in [35]. According to Garcia-Baccino et al. [36] and Aguilar et al. [37], *A*−<sup>1</sup> was computed accounting for inbreeding in order to avoid inflation (bias), particularly for BLUP scenarios, and to minimize blending problems between genomic and pedigree matrices.

All analyses were performed using the *BLUPF90* suite of programs under their standard parameters of quality control for genomic database [38]. Pedigree information from the 23rd to 28th generations was retained (Table 2). Phenotypic data comprised does from the 23rd to 27th generations. Phenotypic records of progenies at the last generation (28th generation) were not considered because they represented young candidate animals. Although all scenarios (S1–S7) used genotypes of grand-sires and sires (26th and 27th generations) for genomic prediction, a few scenarios (S1 and S2) also used all genotypes of grand-dam and dams. In the other scenarios (S3–S7), genotypes from grand-dams and dams were discarded due to low IA and the large number of missing SNPs [39]. On the other hand, in S3, EBVs of non-genotyped progeny were estimated as means of their parents' EBVs (Table 1).

The accuracy of gEBVs was measured by the Pearson's correlation between predicted and true breeding values (TBVs) on the progeny at the validated population, animals belonging to the 28th generation. To assess the gain of accuracies when genomic information was introduced, the mathematical differences of accuracies between each scenario (S1–S7) and BLUP were calculated within each simulation.

This study emphasizes the genetic improvement via doe selection. Hence, response to selection was calculated by subtracting the TBV mean of young candidates (1500 does), within each imputation strategy, from the TBV mean of the 150 selected individuals. We also computed the percentage of candidate animals correctly selected. Results within each scenario are exposed, comparing IA, the accuracy of gEBVs, the selection response and genotyping cost.

#### **3. Results**

#### *3.1. Simulation Outcomes*

The outcomes from the 36 simulations were similar with regards to the response to the selection of experiments using Spanish rabbit lines and the BLUP method for genetic evaluations: an average of 2.53 kits after 20 generations of selection [4,5]. For genomic selection, the average response to selection across QTN models computed for the first two generations with HD platforms was 0.15 kits (0.13 and 0.17 for QTN\_5 and QTN\_44 models, respectively), which corresponds to a response per generation of 0.08 (ranging between 0.002 and 0.20) across QTN models. These results are in line with outcomes of pig genomic selection experiments [40], as hitherto there have been no empirical genomic selection experiments in rabbits. As expected, the additive genetic variance was smaller in QTN\_5 than in QTN\_44. As expected, scenarios with 44 QTNs per chromosome presented higher additive genetic variance than QTN\_5 ones. This is strictly correlated with the number of QTNs that were fixed during haplotype creation. On the other hand, this discrepancy did not have a clear effect on the response to selection as reported on Figure 2.

**Figure 2.** Relationship between initial additive variance and the response to genomic selection. The initial additive variance is the additive genetic variance at the 25th generation in which the genomic evaluations begin. QTN\_5 stands for the QTN model with 5 QTNs per chromosome and QTN\_44 stands for the QTN model with 44 QTNs per chromosome. The simulated rabbit genome comprised 20 chromosomes.

#### *3.2. Performance of Imputation Strategies*

The average of IA for the animals in the validated population (28th generation) was computed for each scenario. The results are shown in Figure 3. No difference, in terms of IA, was shown between two different QTN models. The greatest IA (0.99) was achieved in the S1. The IA decreased to 0.941 (S2) when the LD platform was used. A great decline in IA was found in S3\_A (0.797) when the half validated population had imputed the wholegenome, unlike S3 (0.935). IA decreased to 0.918 when dams of the training population were genotyped with the MD platform (S4). S6 also presented an intermediate value (0.902) between all scenarios. Using LD platforms on does, IA decreases to 0.858 (S5); whereas IA declines to 0.811 with non-genotyped does (S7). In addition, lower standard deviation was obtained in S1 (0.0003), S2 (0.002) and S3 (0.002). Conversely, S7 and S5 presented the highest values of standard deviation (up to 0.005 and 0.004, respectively).

In the same way as IA, genotype yields across QTN models were higher in S1 (0.999) and S2 (0.984). The values were intermediate for S3 (0.952), S3\_A (0.951) and S4 (0.945). Lower genotype yields were presented in S6 (0.899), S5 (0.893) and S7 (0.866).

#### *3.3. Gemomic Prediction vs. Pedigree-Based Analyses*

The accuracy of gEBVs presented, on average, a higher accuracy of prediction for QTN\_44 than for QTN\_5 (Figure 4). The accuracies of gEBVs of S1 were 0.26 ± 0.015 (QTN\_44) and 0.228 ± 0.014 (QTN\_5), representing mean ± standard error. S2 presented lower accuracy for both QTN models, 0.237 ± 0.014 (QTN\_44) and 0.205 ± 0.013 (QTN\_5). S3 presented similar values compared to S2, 0.237 ± 0.016 (QTN\_44) and 0.193 ± 0.015 (QTN\_5). S4 and S6 were very similar to S3 in QTN\_44, with 0.232 ± 0.016 (S4) and 0.234 ± 0.016 (S6), and had slightly higher values than S3 in QTN\_5, with 0.197 ± 0.015 (S4) and 0.199 ± 0.016 (S6). Conversely, lower accuracy values were found for S5 and S7, with 0.223 ± 0.016 (S5) and 0.22 ± 0.015 (S7) for QTN\_44, and 0.185 ± 0.014 (S5) and 0.18 ± 0.015 (S7) for QTN\_5. The accuracy of gEBVs drastically declined in S3\_A, 0.09 ± 0.016 (QTN\_44) and 0.09 ± 0.013 (QTN\_5), because of low IA presented in progeny

with its imputed whole-genome. The S3\_A values were even worse than BLUP accuracies, which showed 0.202 ± 0.014 (QTN\_44) and 0.166 ± 0.015 (QTN\_5).

**Figure 3.** Imputation accuracy for each imputation strategy within QTN models. "S" stands for scenario. QTN\_5 stands for the QTN model with 5 QTNs per chromosome and QTN\_44 stands for the QTN model with 44 QTNs per chromosome. The simulated rabbit genome comprised 20 chromosomes. Each dot stands for a simulation value within each scenario.

**Figure 4.** Accuracy of the estimated breeding values (EBVs) for each scenario within QTN models. "S" stands for scenario. QTN\_5 stands for the QTN model with 5 QTNs per chromosome and QTN\_44 stands for the QTN model with 44 QTNs per chromosome. The simulated rabbit genome comprised 20 chromosomes. Each dot stands for a simulation value within each scenario. White rhombus stands for the mean by scenario.

The results presented a great variability on simulation-based analysis (on average S.D. of 0.08). Figure 5 shows the mathematical differences of accuracy between all genomic scenarios (S1–S7) and BLUP scenarios for each simulation. S1 presented better results than BLUP over 36 simulations for QTN\_5, whilst it had better results in 95% of simulations for QTN\_44. S3 outperformed in 82% of simulations, whilst S2 outperformed in 75% for QTN\_5. By contrast, S3 was in the 75% of simulations better than BLUP for QTN\_44, whilst S2 was in the 80%. S4 and S6 performed in a rank of 70–80%, whereas S5 and S7 performed in a rank of only 50–60%. S3\_A was better than the BLUP scenario in only 25% of simulations.

**Figure 5.** Histograms of each scenario representing the differences between accuracy from genomic predictions (ssGBLUP) versus genetic prediction (BLUP) per simulation. Right-sides of red dashed lines represent the number of simulations in which genomic selection outperformed BLUP.

The response to genomic selection and the percentages of correctly selected candidates are showed in Table 3. These parameters are strictly correlated with the accuracy of gEBVs in prolific livestock; however, they represent more pragmatic methods of comparison and dissemination for farmers. The correlations were 0.90 and 0.88 between the accuracy of gEBVs and response to genomic selection, and the first one and percentages of animals correctly selected, respectively. A response to selection of 0.105 was found when only pedigree information was used. In S1, the number of kits per generation increased by 22% with respect to BLUP. A slight decline was observed in scenarios S2, S3, S4 and S6 with a value of 0.120, 0.117, 0.114 and 0.116, respectively. S5 and S7 did not show any significant augmentation with respect to the BLUP scenario. As expected, a lower selection response than BLUP was noticed in S3\_A (0.046). The percentages of correctly selected animals were similar to the trend of selection response, thus best performance was obtained in S1 with a value of 30.54%, followed by S2 and S3 with a percentage of 29.45% and 29.36%, respectively. S3 to S7 presented a range between 29.06% and 28.42%. BLUP presents a percentage of correctly selected animals of 27.81%. Even for this parameter, the values in S3\_A were lower than in the BLUP scenario.


**Table 3.** Response to genomic selection and percentages of correctly selected animals across QTN models for each scenario.

<sup>1</sup> Mean of response to genomic selection (kits). <sup>2</sup> Standard error of response to genomic selection. <sup>3</sup> Percentage of animals correctly selected. <sup>4</sup> Standard error of percentage of animals correctly selected.

#### *3.4. Genotyping Costs*

Figure 6 shows the relationships between the price of genotyping and IA, and the accuracy of gEBVs. As expected, the more investments, in terms of genotyping, the more IA and accuracy of gEBVs were found. However, that is not a strictly linear correlation, especially for IA. S1 showed the most expensive investment (EUR 112,000) due to progeny genotyping at the HD platform. The second ranking position was S2 with a genotyping cost of EUR 53,500. S3, S4 and S6 presented a lower cost compared to S2 from EUR 38,500 to EUR 31,000. S5 and S7 are the cheapest scenarios, with costs of EUR 26,800 and EUR 23,500, respectively.

**Figure 6.** Plots of relationship between price of genotyping and imputation accuracy (top figure), and price of genotyping and EBV accuracy (bottom figure). Vertical lines are the mean ± standard error.

#### **4. Discussion**

Our results highlight four main points for discussion: (1) imputation strategies; (2) causes that affect genomic prediction; (3) comparison with other studies; and (4) searching for trade-off: cost and genetic accuracy trends.

#### *4.1. Imputation Strategies*

Genotype imputation was introduced as a tool to increase detection power on association studies for linking results across studies that rely on different genotype platforms [41]. On the other hand, imputation can be a suitable tool to reduce the cost of genotyping in both plant breeding [30,42] and animal breeding programs [18,43,44]. In this case, close ancestors are typically genotyped at HD platforms, whilst progeny are genotyped at lower SNP density. When the marker position of different platforms perfectly overlaps and pedigree information is present, imputation accuracies (IAs) are usually high [29]. In our situation, in which all close ancestors are genotyped with the HD platform, this approach may not produce the same benefit as in the other species due to the prohibitive cost of HD platforms. For this reason, we examined how IAs are affected by the different genotype information of the ancestors. There are several factors influencing IA, highlighting the missing rate of LD platforms as one of the main factors [22,30]. Missing rate is the percentage of SNPs present at an HD platform that are untyped (not covered) at LD platforms. The missing rates were 70% and 99.7% for MD and LD platforms in the current study. Although the first imputation studies suggested values between 50% and 75% [29,45], the values of the current study are in line with other studies that presented very high IA for both plant [42] and animal breeding programs [18,43,46]. The missing rate can be even higher when pedigree information is available for the imputation process. AlphaImpute, which uses both genotype and pedigree information, has been demonstrated to have high imputation performance [18,29,43], even in populations with low levels of linkage disequilibrium [47]. Thus, software and missing rates seem appropriate for imputation in rabbit breeding programs.

The number of animals in the reference population, as a second factor, influences the resolution of haplotypes during the phasing process. A large number of individuals in the reference population ensures high IA in validated populations. However, it can be reduced considerably if the animals from both populations are close relatives, sharing the structure of linkage disequilibrium and haplotypes across wide chromosome segments [22,29]. S7 showed a high IA using only 70 male ancestors, which is explained by the close relationship between them, the female ancestors and progeny under ongoing selection. As expected, S1 presented the highest IA and genotype yields due to the larger number of animals in the reference population (training) and lower missing rate. In general, many studies showed that imputation using HD platforms on ancestors and MD platforms on the progenies produces high levels of IA and concordance rate in dairy cattle [46,48,49], pigs [13,18,19], poultry [12,50], sheep [20] and farmed Atlantic salmon [43]. When the SNP densities of MD (S1) platforms were reduced to LD (S2) in the validated population, IAs decreased only five percentage points. Hence, vast haplotype information keeps retrieving, using LD platforms, because of the high relatedness of rabbits. The IA results of S1 and S2 were similar to those reported in a pig imputation study, especially for Landrace and Yorkshire breeds [19]. In the current study, fewer differences in IA were found when grand-dams were genotyped at MD platforms (less than two percentage points, S3 and S4) compared to S2. This demonstrated that dams at MD platforms were enough to retrieve female haplotypes and to keep high IAs, being noticeable when S6 presented better IA values than S5. The female haplotype resolution is better when dams are genotyped at a higher SNP density than any level of SNP density on grand-dams. On the other hand, the imputation of whole genomes (S3\_A) seems to not be a feasible technique for rabbit breeding programs. IA declined up to 0.797, probably due to the small number of training populations and the higher error rate of imputed SNPs associated with QTNs—very important if they have low minor allele frequency (MAF). Many animals with imputed whole-genomes presented

very low IA—less than 60%. Conversely, some studies showed the benefits of strategies based on imputed whole-genomes in part of training populations [47,51] and evaluated populations using a larger number of genotyped individuals [51].

#### *4.2. Causes That Affect Genomic Prediction*

In this study, ssGBLUP was used as a method for genomic evaluations. This method includes a different sort of information: phenotypes, pedigree, and genotypes. Previous studies demonstrated that ssGBLUP gives potentially more accurate and less biased genomic gEBVs than multistep methods, especially in the presence of small populations and sex-limited traits [31,52].

The accuracy of the EBVs varies greatly within each scenario with respect to IA (Figures 3 and 4). Conversely to the IA, which is mainly influenced by the genotyping strategies, the accuracy of gEBVs is affected by several factors. Some of these are related to the genetic architecture of the traits such as QTNs distribution and their allele frequency [22,53]. The number of SNPs in LD with these QTNs and the allele frequency of those SNPs also played an essential role in the accuracy of genomic prediction [54]. In addition, working with imputed genotypes, it is important to consider the influence of the imputation error of those SNPs and therefore IA [22,53]. In these studies, scenarios with a higher IA value also presented higher accuracy on genomic prediction and the opposite. However, a strong correlation between these two-type accuracies cannot be defined due to the high variability of gEBV accuracy within the scenarios. Additionally, with S3\_A and S7, it was confirmed that low IA values (under 85%) can be deleterious for genomic prediction. For that reason, conservative thresholds for genotype imputation and quality control before genomic selection must be adopted. Furthermore, Cleveland and Hickey [18] showed that differences in gEBVs may be due to both different values of IA but also to the intrinsic genotyping structure of each scenario. There may be animals that have a marginal influence on IA but can significantly affect gEBVs.

As mentioned previously, variability in these scenarios was also caused by different allele frequencies of SNPs and QTNs. High heterogeneity of these factors was observed between simulations due to the random events that occurred during the 28th generation of mating. In some simulations, a larger number of QTNs were fixed, and in some, many SNPs were not associated with the rest of the QTNs. This would also explain the variability between simulations, and why few simulations presented a lower accuracy of gEBVs for genomic scenarios compared to BLUP, even in the scenarios with high IA as in S2 and S3.

Regarding QTN distributions, QTN\_44 presented better performance than QTN\_5. This trend agrees with the study of Zang et al. [53], in which ssGBLUP outperformed when the phenotype was controlled by several genes of equal effect sizes. A method that assumes unequal variances for each marker could suit for the genomic prediction of QTN\_5. Modeling SNP's effect and its variance can potentially give better results in short term selection for this simulation. Iterative ssGBLUP and/or nonlinear weight **A** can easily be implemented and can potentially lead to an increase in prediction accuracy [55,56]. Fernando et al. [57] also propose a single-step Bayesian regression in which it is possible to model the distribution of marker effects in many forms such as *t* distribution, variable selection model, and mixture distributions. Despite this, the method of estimation and of building **G** matrices was kept the same in all simulations for the sake of simplicity and for comparison with other studies.

#### *4.3. Comparison with Other Studies*

S1 represented the typical genotyping strategy used in all livestock species in which parents are genotyped with HD platforms and progeny are genotyped at LD platforms. As expected, S1 exhibits the best accuracies of gEBVs, and the accuracy of genomic prediction presented in S1 is close to that presented when the same candidate animals are genotyped with HD platforms; correlation is almost one in all simulations for both QTN models (data not shown). A similar correlation was found in pigs [18,58] and cattle [21,22,53],

even if the proportion between HD and LD animals was different due to distinct breeding schemes present in these species. Previous imputation studies conducted in commercial pig breeding programs are a good comparison method due to a similar mating system, although the number of animals is lower in rabbits. Scenarios comparable to those present in this study have been reported in Cleveland and Hickey [18]. A similar sharp drop in IA and the accuracy of gEBVs was observed when animals of validated populations were genotyped with MD to LD platforms. In addition, an analogous decline in gEBVs was observed when grand-dams and dams were genotyped with LD or MD, although in both studies IA was quite high for these scenarios. However, Cleveland and Hickey [18] demonstrated that genotyping at HD platforms for animals that are not related with LD animals has little impact concerning IA, but it can significantly affect the accuracy of gEBVs. Nevertheless, this scenario was not included in our study because the limited number of animals would not guarantee the same gap of accuracy presented in pigs. The same consideration made on Grossi et al. [19] can be also made for our study; genotyping reference animals with the LD platform can ensure good accuracy levels of IA and gEBVs when parents are genotyped with HD, if not low levels of prediction have been observed, i.e., S5 or S7.

#### *4.4. Searching for Trade-Off: Cost and Genetic Accuracy Trends*

Cost evaluations are usually performed by comparing imputation strategies against an idealistic genomic selection in which all candidate animals are genotyped at HD platforms. Under this condition, imputation strategies are always a cost-effective technique [18,42,43]. Here, we show only genotyping costs to evaluate investment for moving from traditional genetic evaluation (BLUP) to genomic selection (ssGBLUP). Imputation was further demonstrated as a reliable tool for reducing genotyping investment, with an accuracy higher than BLUP in the majority of cases.

However, some scenarios present a cost-benefit inefficiency when compared with the other scenarios. Clearly, S3\_A cannot be taken as an option due to the low performance of prediction—even worse than BLUP. Results from this scenario suggest that a method based on high thresholds of individual-specific IA must be considered, especially when genomic characterization is available [22], and for strategies with non-genotyped animals [51]. Using animals with imputed whole-genomes for genomic selection is still a challenge due to the imputation error rate [51,59]. In this sense, high error rates of imputed SNPs associated with QTNs are clear in S3\_A, presenting 75% of simulations lower than BLUP scenarios. BLUPF90 reported genotypes duplicated for S3\_A, as *AlphaImpute* copied wholegenomes of full-sibs for some progenies.

The same reason can be applied to S7 and S5 as, even if lower genotype investment was present in these scenarios, these strategies cannot be considered one of the best scenarios because the number of simulations in which genomic selection outperformed BLUP is limited—about half of the cases. Conversely, S3, S6, and S4 presented a considerable increase in accuracy and response to selection with a marginal increment investment.

The opposite situation was presented in S1, and even if a significant reduction in accuracy was observed switching from MD platforms to LD platforms (S2), the big drop in price observed can potentially justify this decline in accuracy. No significant differences in terms of percentage of correctly selected animals and response to selection were observed between S2 with S3 and S6, thus these two scenarios are preferable over S2 due to the lower cost. In addition, S4 can be discarded among the best scenarios, as it presents a lower selection response than S6 and S3 but at a higher price. In conclusion, the best scenarios can be identified as S3 and S6. Genotyping only half of the animals per litter (S3) is an effective strategy to reduce the impact of the cost of genotyping, especially in prolific species such as rabbit and pig. This approach is commonly adopted in the pig industry selection scheme [15,17,18]. On the other hand, S6 demonstrated that genotyping grand-dams (S4) leads to a negligible difference in accuracies when sires and dams are genotyped, with S6 being cheaper than S3 (+ EUR 6750). Thus, we considered S6 as the strategy with the best trade-off for implementing genomic selection in rabbits.

#### **5. Conclusions**

Imputation strategies are feasible in rabbit breeding programs as in other species, as IAs were high, particularly in scenarios with parents' genotype at HD platforms. Furthermore, the slight positive correlation between IA and the accuracy of gEBVs was also demonstrated, and scenarios with IA also have a high gEBV accuracy. The results of the response to genomic selection and the percentage of correctly selected candidates are in the same line as genetic accuracy values. Despite this, gEBV accuracy showed great variance among simulations. Therefore, we must be cautious with the results from simulated data as they are conditioned to several factors of the small reference population (e.g., size, relationships between individuals, inbreeding level, and update). Another comparison between genomic selection and BLUP could be made by enlarging the reference population using dams of the nucleus farms and the multipliers, or even crossbreeding from commercial farms. In conclusion, the adoption of imputation strategies can be an effective strategy for drastically reducing the genotyping cost in rabbits, maintaining an accuracy slightly lower than that with all HD animals. Hence, the best trade-off scenario can be identified in S6; although, as previously stated, these results may change under different conditions.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2076-261 5/11/3/803/s1, Table S1: Parameters and info to run simulations in *AlphaSim* and *MaCS*.

**Author Contributions:** Conceptualization, N.I.-E. and A.B.; methodology, B.S.S.-M. and N.I.-E.; software, B.S.S.-M. and E.M.; validation, B.S.S.-M., E.M. and N.I.-E.; formal analysis, E.M. and B.S.S.- M.; investigation, B.S.S.-M.; resources, A.B. and N.I.-E.; data curation, B.S.S.-M. and E.M.; writing original draft preparation, E.M. and B.S.S.-M.; writing—review and editing, E.M., B.S.S.-M., N.I.-E. and A.B.; visualization, E.M. and B.S.S.-M.; supervision, N.I.-E.; project administration, B.S.S.-M.; funding acquisition, N.I.-E. and A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research study was funded by AGL2017-86083-C2-P1 from "Plan Nacional de Investigación Científica" of Spain-Project I+D. B. Samuel Sosa Madrid was supported by FPI grant, number BES-2015-074194, from "Ministerio de Ciencia e Innovación".

**Data Availability Statement:** The databases used and analyzed in the current study are available from The Figshare Repository. Some scripts are also included in the repository. The link addresses are https://doi.org/10.6084/m9.figshare.13350443 or https://figshare.com/s/2591beabd3ef5bc6d220 (accessed on 9 March 2021).

**Acknowledgments:** We express our profound thanks to Francisco Rosich of "Área de Sistemas de Información y Comunicaciones" (ASIC) for his contributions to pipelines of simulations and imputation analyses.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

