Genetic analysis of the effects of heat stress before and after lambing on pre-weaning live weight in Spanish Merino lambs.

BACKGROUND: Heat stress (HS) is a major environmental effect on sheep production. Hereby, we estimated the genetic (co)variance component of HS on the pre-weaning performance of 19,022 Merino lambs by analysing the climatological index of temperature and relative humidity (recorded 30 days before lambing and after lambing) using transversal and longitudinal mixed linear models. METHODS AND
RESULTS: The global impact of HS during the last 30 days of pregnancy was -17% for birthweight and ranged between -4% and -8% for live weight at 15, 30 days of age (W30), and average daily gain from birth at 30 days. The results from both statistical approaches showed very similar heritabilities (h2 ), ranging from 0.192 to 0.237 for the direct genetic (D) effects and from 0.072 to 0.082 for the maternal genetic (M) effects, but the antagonism between (D) and (M) was higher when a longitudinal model was used. A significant genotype-environmental effect was also found regardless of whether the climatological covariables were considered in the model. In addition, we employed D and M breeding values for W30 as an example to create a new subjacent index by first using a principal component analysis and employing the leading eigenvalues as a weighted factor that provides the information needed to identify those genotypes that maximise the response for both genetic effects over a wide range of climate-environment levels.
CONCLUSIONS: Our study revealed that the HS indexes of the mother during the gestation period have a significant effect on the growth of the lambs during the early stages of life.

INTRODUCTION

Ten years ago, Segnalini et al. (2013) warned about the likely negative effects on animal production of the marked changes in climate in the Mediterranean region. This negative effect, which was already detected in goats (Finocchiaro et al., 2005; Menéndez‐Buxadera et al., 2012), was also demonstrated in dairy sheep and dairy cows (Bernabucci et al., 2014; Carabaño et al., 2017). But this change in the average weather conditions could be even larger in Spain, where the temperature has increased twice as much as the European average, particularly in the central and southern regions, where most of the flocks of small ruminants are located (Martínez‐Linares, 2009).

The nature of the semi‐automatised recording systems used in dairy animals bred in commercial herds allows its merging with climatic records of the official meteorological stations. In this way, large databases aiming to elucidate the relationship between climate and production over a given period can be built. This strategy was commonly used in numerous studies carried out modelling the genetic variability in heat stress (HS) tolerance in beef and dairy animals (Carabaño et al., 2017), but in most of them, the approach is primarily focused on the direct effects of HS on the animals. However, it should also be taken into consideration that HS may affect the animal earlier if it occurs during the gestation period since it can affect female productivity (Menéndez‐Buxadera et al., 2020; Tao et al., 2012). These latter references showed a negative effect on daily milk production in cattle subjected to HS during pregnancy, particularly during the first part of lactation. But if this negative effect would also be present in meat animals, it will be even greater since the newborn depends basically on its mother's milk during the first weeks of life.

Small ruminants play a significant role in the livelihoods of a considerable part of the human population from socio‐economic aspects (Herrero et al., 2013; Kosgey & Okeyo, 2007). In this sense, meat sheep are considered to be locally adapted, being able to produce even in harsh environmental conditions, and for this reason, they are exploited in marginal areas (Degen, 2007). However, there is also evidence showing that HS has an important negative effect on productive performance and well‐being in this species (Silanikove, 2000), even affecting the viability of the offspring during gestation and the first days after birth (Everett‐Hincks & Dodds, 2008; McCrabb et al., 1993). The previous elements are coherent with the importance of HS in this animal's performance; however, the origin and importance of this effect must be considered during two separate but interacting periods: gestation and the first weeks of life. If these effects are important but not taken into account in the statistical model, the estimates of the animals’ breeding value (BV) could be biased. To our knowledge, previous studies have not been conducted in this aspect.

The picture is even more complex in small ruminants dedicated to meat production. During the pre‐weaning period, the importance of HS depends not only on the direct capacity of the animals to modify their performance by adapting to environmental stress but also on the ability of the mothers to care for and nourish those newborns. In this sense, McCrabb et al. (1993) demonstrated that the HS effect on the mother during pregnancy is another important additional environmental factor to be considered in evaluating the genetic response to harsh environments. As a result, the direct influence of HS on live weight during the first weeks of life needs to be modelled by taking into account all the genetic components, such as the reaction norm (RN) model for direct and maternal effects along the trajectories of the climatological variable (which affects the performance of the animals), the classical permanent environmental maternal components and an additional environmental indirect fixed effect due to the HS during pregnancy of the mother. The fact that all those effects are complex to model together may be the cause of the lack of published results on genetic components for HS in growth or other development traits in this species.

Those traits whose responses are continuous functions of other variables were termed ‘function value traits’ (FVT) by Stinchcombe et al. (2012), who also recommended the use of a random regression model (RRM) to estimate its genetic (co)variance structure due to its intrinsic complexity. In this context, the solution of an RRM of FVT allows us to estimate the genetic parameters and BV for all the animals for genetic direct and maternal effects across all the points of the trajectory of both continuous variables. However, these results are difficult to handle and interpret due to the vast quantity of information produced. In this sense, the use of a principal component analysis (PCA) could be a useful tool since it allows us to significantly reduce the number of variables, with a minimum loss of information (Hair et al., 2013). However, the use of the eigenvectors of the leading eigenvalues derived from PCA are also mutually orthogonal; therefore, they can be used as a weighted factor in a new subjacent index. This methodology, presented by Togashi and Lin (2006), was successfully used in different scenarios of animal breeding, such as egg production (Venturini et al., 2013); beef cattle (Boligon et al., 2016), dairy cattle (Macciotta et al., 2017) and goats more recently (Menéndez‐buxadera et al., 2021).

With all of these in mind, the objectives of this article are as follows: First, to present evidence of the importance of the HS effects before and after lambing on daily growth and live weight in the first weeks of life in meat sheep, using the Merino as a model of the most cosmopolitan breed in the world; second, to estimate the evolution of the genetic (co)variance components of these pre‐weaning growth traits along the trajectory of a temperature and relative humidity (HR) index using a RN model; third, to identify the importance of these effects of HS on the BV of the animals, compared with the classical evaluation model used in the breeding program for this breed. Finally, we demonstrate the feasibility of the use of PCA to reduce the number of variables with a minimum loss of information and thus facilitate the selection process in the breeding program of Spanish Merino sheep.

MATERIAL AND METHODS

Animals

For this study, we accessed the periodic records from the live weight data system officially used by the Merino sheep breeders' society (SCME) in 110 flocks participating in the breeding program of this breed, from which 27 were selected according to being located closer than 15 km of distance to one of the official meteorological stations from the Spanish weather agency. In total, data from 20 different stations were included in this study. Following the recording protocol, three individual estimations of live weight were made for each animal during the first 3 months of age, in addition to the birthweight (W0). These weights allowed us to adjust the live weight at 15 and 30 days (W15 and W30, respectively) as well as the daily gain (DG) between birth and 30 days. The data were edited by deleting records outside the range of ± 3 standard deviations, and finally, the data from 19,022 lambs born over 10 years from 9778 mothers (1698 in the data vector) and 402 parents (95 in the data vector) were used in the study. The pedigree was made with all the available data from the animals and their ancestors in the SCME, making a total of 36,156 animals. The climatic variables were obtained from the official State Agency for Meteorology stations and were the average temperature ambient (T) and HR during the last 30 days of gestation and 30 days after lambing, respectively. These climatic variables were combined to form an index called THI, and the following formula was applied: where T and HR are the average temperature and HR during the indicated period before and after lambing and are labelled THI_b and THI_a, respectively.

Statistical analysis

The data were first studied using a linear fixed‐effects model to identify the most important factors affecting the dependent variables. In this case, we considered the contemporary group as a combination of the flock, year and bi‐monthly lambing (CG, i = 271 classes with more than 10 records in each level); age to lambing of the dams (ED and j = 10 classes of 1 year); sex of the lamb (S with k = 2 levels); litter size (LZ, l = 2 classes). For all dependent variables, the inclusion of THI_b (29 classes, THI = 43 to 72) and THI_a (31 classes, THI = 41 to 72) as fixed continuous variables increased the model's determination coefficients by 3% to 5%.

Based on these preliminary results, the final statistical procedure was performed, which consisted of the comparison of several mixed linear models in which the random effects were modelled differently. The following models were applied using Asreml3 (Gilmour et al., 2009):

In these models, represents the estimates of the W0, W15, W30 and DG for each animal, and fixed contains the same fixed effects (GC, ED, S and LZ) mentioned above. Model 0 is equivalent to the model currently used in the genetic evaluation of this breed, while Models 1 and 1a incorporate the fixed effects of the climate–environment evaluated in two distinct stages over time: during the last 30 days of gestation (b 1 = THI_b) and in the first 30 days of life (b 2 = THI_a), m and p l are random components for the kth animal and lth permanent environmental effect, respectively. The solution of these two models allows us to show the importance of HS before and after lambing on the dependent variables during the first weeks of life adjusted by all the fixed and random effects included in the model compared with the current procedure used in this breed's breeding program, which does not consider the influence of these climatological effects on the dependent variables studied. Finally, Model 2 is a RN model that assumes the existence of the heterogeneity of the (co)variance components of the direct effect of HS on the animal along the scale of THI_a, adjusted by the aforementioned fixed and random effects, as well as the effects of the maternal climate–environment during the last 30 days of pregnancy THI_b. These last two factors were also considered to be continuous variables and were modelled using fixed Legendre polynomials and of order r = 3 and r = 1 for THI_a for the random direct and maternal effects. Models 1a and 2 contain the same fixed effects but are modelled differently: direct (d) and maternal (m) genetic random effects and their covariances (dm) are assumed to vary along the scale of THI_a in Model 2, while they were homogeneous in Model 1a. The effect of the permanent maternal environment (p) was similar in all three models.

The incidence matrices X 1, X 2, Z 1, Z 2 and Z 3 contain elements 0 or 1 to connect each observation with the effects indicated in the subscripts. In Model 2, the Z 1 and Z 2 matrices were replaced by the coefficients of a Legendre polynomial () of order r = 1 for all genetic effects to express the trajectory of THI_a on a scale between −1 and +1. In matrix form, all the models have the same general representation: where A is the numerator of the relationship matrix (n = 36,156 total animals). I and I are the identity matrix for the random effects of the maternal permanent environment (dimension p = 9778 sheep) and residual error (dimension n = 19,022 records), respectively, and ⊗ is the Kronecker product symbol. In Models 0, 1 and 1a, we will have: where G 0 is a 2 × 2 matrix with elements containing variances for direct genetic effects (), maternal effects (), and their covariance (); is the variance of permanent effects of the maternal environment; and is the residual variance common to all observations. Model 2 estimates the same (co)variance components as Models 1 and 1a, but along the THI_a climate scale, using a complex structure that requires an alternative formulation:

In this case, the K random regression coefficient matrices are incorporated into the genetic components, and in this RN model, the K matrix contains the following elements:

Note that K is a symmetric matrix consisting of four sub‐matrices with the same components as in G 0 but expressed by a linear RN containing the (co)variances for direct genetic effects (K), maternal effects (K) and their covariance (K = K) with their respective variances of the intercept ( and ); slope ( and ) and its covariance (; ; and). In this case, subscripts d and m represent the genetic variances of direct and maternal effects, while o and s indicate the intercept and slope of RN, respectively. By manipulating the elements of these matrices, as well as the coefficients of the Legendre polynomial (), the (co)variance genetic components can be estimated along the trajectory of THI_a, following the procedure presented by:

With this formulation, the heritability of the direct () and maternal effects () and genetic correlations between the two effects (r) are estimated by applying classical formulas (Falconer & MCKay, 1996).

BV estimation along the scale of the climatological variable and the use of PCA for each dependent variable

The BVs from Models 0, 1 and 1a of each trait were estimated by the direct solution. In all these cases, the BVs were expressed in a transversal way, assuming no genetic variation in the scale of THI_a.

The solution of RN Model 2 is different, and the results are presented by two genetic functions (gf): for direct effects, gf = [a], and for maternal effects, gf = [a], in each animal. The terms inside the brackets indicate the intercept (o) and slope (s) for direct (d) and maternal (m) genetic effects, respectively. These gf values allow us to estimate the BVs for genetic direct and maternal effects (BVd and BVm, respectively) for each dependent variable along the trajectory of THI_a following the procedure formulated by Jamrozik and Schaeffer (1997): where Φ are the corresponding coefficients of the Legendre polynomial for each ith point in the climatic variable THI_a. As a result, a total of 31 estimations for BVd and a similar number for BVm will be available. This high number of parameters is not easy to handle and, in practice, makes it more difficult to select the most suitable animals. Under these conditions, the use of a PCA could be highly recommended for treating these complex longitudinal genetic (co)variance structures for both genetic effects because it allows us to reduce the number of variables to a greater extent, with a minimum loss of information. PCA was performed using MATLAB software (2019), and the results were used for three different purposes:

First, we identified the existence of genetic (co)variation in the form of a response to the adaptation of animals to the increased level of HS using the elements of the K matrix from Model 2.

Second, to identify the relative variability and directions of BVd and BVm at each point of the wide range of the climate–environment variables measured by THI_a, the result is illustrated by a biplot and Q–Q plot.

Finally, to present an alternative way of facilitating the breeding process, the results BVd and BVm for a wide and representative range of THI_a levels were added (BVt = BVd + BVm), as suggested by, and these new datasets were subjected to a new PCA. The eigenvector (ev of the jth leading eigenvalues () was used as a weighting factor, and a new subjacent index was created for each animal (Ipct), according to Togashi and Lin (2006): where ev are the corresponding eigenvector coefficients and BVt is the total standardised original BV at the selected points on the THI_a scale. This new subjacent index summarised all the total variability of BVd and BVm expressed in the leading. The application of an RN model adds greater complexity, but as we will show later, it can provide important information for the breeding program.

RESULTS

General indicators

The distribution of climate records before and after lambing is presented in Figure 1. In general, this amplitude of THI levels is representative of the climatic conditions of the area where the flocks in this study are located (mainly in the south and southeast of the Iberian Peninsula). As expected, the highest THI values were recorded in July and August, and the lowest values were recorded in December and January of each year. Table 1 shows the basic statistical information for the studied trait together with the HS indicators.

FIGURE 1

Distribution of the number of records and monthly averages of the THI climatic index before and after lambing

TABLE 1

General indicators (and s.e.) from the analysis of heat stress (HS) effects on growth and live weight in the first weeks of life in Spanish Merino sheep

Indicators	Number of records	Average ± s.e.
Birthweight (kgs)	18,897	3.96 ± 0.88
Adjusted weight at 15 days (kgs)	17,836	7.49 ± 1.51
Adjusted weight at 30 days (kgs)	17,854	10.99 ± 2.46
Daily gain (DG) 0 to 30 days (kgs/day)	17,328	0.238 ± 0.07
Temperature before lambing (oC)	19,022	14.74 ± 6.65
Temperature after lambing (oC)	19,022	15.01 ± 4.67
Relative humidity (HR) before lambing (%)	19,022	65.62 ± 14.44
HR after lambing (%)	19,022	66.10 ± 10.87
THI before lambing	19,022	56.9 ± 9.07

Effects of climatic variables before and after lambing on lamb growth

The solutions of Models 1 and 1a allowed us to estimate the deviations of each dependent variable along the trajectory of THI before and after lambing (Figure 2). Both models showed the depressive effects of the climate on the variables studied, but they were evident only when these two factors were considered together on the live weight of the animal (Model 1a). In this sense, it is worth mentioning the lack of changes in the form of response between THI_a = 42 and THI_a = 59, as well as the marked depressive trend in animal performance and from THI_a = 60 up to THI_a = 72. It should be noted that the response curves shown in this figure are adjusted for all the genetic and environmental variations included in Models 1 and 1a.

FIGURE 2

Evolution of adjusted deviations of live weight at 15 and 30 days of age and daily weight gain in Merino lambs during the first weeks of life, adjusted for the direct (lamb's first week of life) and maternal (during the last 30 days of gestation) climate effects

To quantify the global impact of the maternal environmental HS on the live weight of their offspring, we created two climatical zones based on the original dataset: hot (THI_b > 66) and cold (THI_b < 55). This new dataset was analysed using a fixed linear model including the same fixed effect previously described (GC, ED, S and LZ) and the effect of THI_a as a fixed covariable, except for birthweight. The results were highly significant (p < 0.001) for all dependent variables (Table 2). Overall, HS during the last 30 days of pregnancy had a negative effect of −17% for birthweight and between −6% and −8% for W15, W30 and DG.

TABLE 2

The least‐square means of the effect of two contrasted climatological zones during the last 30 days of pregnancy on the live weight (kg) during the first weeks of life of the newborn in Merino sheep

			Trait
Climatological zone	Number of animals	THI	Birthweight	LW at 15 days	LW at 30 days	DG 0 to 30 days
Cold (28.4%)	4839	43 a 50	3.94(0.26)	7.51(0.51)	11.20(0.85)	0.252(0.002)
Hot (29.1%)	5165	65 a 72	3.24(0.31)	6.93(0.61)	10.46(0.10)	0.232(0.003)

Note: Standard error I parentheses

Covariance components and genetic parameters according to Models 0 and 1a

Neither the (co)variance components nor the genetic parameters showed any marked differences in the parameter estimates for ; or for r when compared with the results of Models 0 (which do not consider the climatic effects) and 1a, which include THI before and after lambing (Table 3). The results of the two models are similar, except for DG, which shows differences. At the same time, the correlations between the BVs estimated by both models were close to unity (results not shown), which leaves the previously presented arguments in doubt. However, it should be noted that the (co)variance components and genetic parameters presented in Table 3 refer to two linear models that express the results in a transversal form for each trait.

TABLE 3

Covariance component estimation and genetic parameters* (heritability for direct and maternal and genetic correlation ) effect of live weight traits in the first weeks of life in Spanish Merino sheep considering no climatic effects (Model 0) and adjusting by THI level in the period 30 days before and after lambing (Model 1a)

	Birthweight	Live weight at 15 days		Live weight at 30 days		DG**
(Co)variance component	Model 0	Model 0	Model 1a	Model 0	Model 1a	Model 0	Model 1a
Direct (d)	0.077	0.243	0.242	0.631	0.663	0.074	0.067
Maternal (m)	0.029	0.091	0.093	0.259	0.262	0.023	0.020
Cov d‐m	−0.032	−0.09	−0.093	−0.259	−0.268	−0.031	−0.027
Permanent	0.056	0.174	0.153	0.341	0.337	0.025	0.025
Residual	0.270	0.876	0.873	2.493	2.464	0.203	0.199
Total	0.431	1.271	1.268	3.466	3.458	0.325	0.293
Heritability d	0.178(0.02)	0.191(0.02)	0.192(0.02)	0.182(0.02)	0.192(0.02)	0.227(0.03)	0.228(0.03)
Heritability m	0.067(0.02)	0.072(0.02)	0.073(0.02)	0.075(0.01)	0.076(0.01)	0.071(0.02)	0.0680(0.02)
Correlation d‐m	−0.693(0.08)	−0.611(0.09)	−0.619(0.09)	−0.641(0.08)	−0.643(0.08)	−0.747(0.06)	−0.730(0.07)

the standard error of the genetic parameters between parentheses.

the variance was multiplied by 10–2 to facilitate the presentation.

Covariance components and genetic parameters according to the RN model (Model 2)

The results of the genetic random regression matrix K are shown in Table 3. Unlike the above results, the K matrix of the RN model (Model 2) allows us to estimate the same (co)variance components, but along the scale of the environmental conditions measured by THI_a but also adjusted by the additional maternal environmental effect of THI_b (the HS during the last 30 days of pregnancy), except for W0, in which an adjustment for climatic effects after lambing does not apply. In addition, each of the traits represented in Table 4 and the corresponding variances for direct and maternal genetic effects ( and ), as well as correlations between both effects (Rdo, mo), will have the same meaning as those presented in Table 3.

TABLE 4

Random regression matrix parameters* (Model 2) for growth traits during the first weeks of life in Merino sheep adjusted for HS effects before and after lambing

	Live weight 15 days	Live weight 30 days	DG**
Genetic direct components
σdo2	0.436	0.979	0.077
σds2	0.150	0.894	0.033
σ_mos	0.014	0.125	0.025
Genetic maternal components
σmo2	0.174	0.432	0.021
σms2	0.046	0.073	0.046
σ_mos	0.050	0.094	−0.014
Correlation between direct and maternal components
Rdo, ds	0.054± 0.11	0.133± 0.19	0.487± 0.10
Rmo, ms	0.560± 0.10	0.754± 0.18	−0.447± 0.20
Rdo, mo	−0.549± 0.11	−0.551± 0.10	−0.570± 0.13
Rds, ms	−0.057± 0.16	−0.643± 0.16	0.009± 0.15
Rds, mo	0.146± 0.29	−0.218± 0.26	−0.590± 0.21
Permanent and residual variance
σl2	0.146	0.317	0.024
σe2	0.836	2.340	2.014

meaning of the symbol in the text.

All (co)variance components have been multiplied by 102 except for and (multiplied by 10).

To quantify the magnitude and direction of the genetic (co)variances manifested along the trajectory of the THI_a level, a PCA was applied to this K matrix (Table 5). Interestingly, the PCA results of the three traits were consistent, with most of the genetic variation (52.8% to 62.8%) accounted for by the first principal component (PC1), while PC2 and PC3 explained most of the remaining variation (between 35% and 46%). These three PCs are orthogonal (non‐correlated but act together) and contribute differentially to identifying the existence of a non‐homogeneous pattern of evolution in the genetic (co)variance structure of this population of animals. This is illustrated in Figure 3, which shows the genetic parameters estimated along the trajectory of THI_a, adjusted by the maternal climate HS level during the last 30 days of gestation (Model 2).

TABLE 5

Eigenvalues and eigenvector coefficients for the first three principal components (PC) from the K matrix (Model 2) for growth traits during the first 30 days of life in Spanish Merino lambs

	LW at 15 days			LW at 30 days			DG 0 to 30 days
	PC1	PC2	PC3	PC1	PC2	PC3	PC1	PC2	PC3
Intercept‐direct	0.907	0.212	0.319	0.749	0.503	0.405	0.816	0.370	0.282
Slope‐direct	0.008	0.848	−0.513	0.474	−0.839	0.220	0.412	−0.304	−0.857
Intercept‐maternal	−0.419	0.458	0.632	−0.440	−0.105	0.806	−0.167	−0.513	0.152
Slope‐maternal	−0.043	0.160	0.485	−0.144	0.176	0.371	−0.369	0.712	−0.403
Explained var (%)	62.8	20.5	14.5	52.8	34.7	12.5	57.5	32.2	7.9
Cumulative var (%)		83.3	97.8		87.5	100		89.7	97.6

FIGURE 3

Evolution of heritability estimates for direct and maternal effects and genetic correlation between both for weight at 15 days (a), weight at 30 days (b) and daily gain (c), as well as direct and maternal genetic correlations of the same trait along the scale of THI_a (d, e and f, respectively) using RN Model 2

The heritability estimates for both genetic effects showed a similar pattern for W15, W30 and DG. The main differences were evidenced in r, which increases their antagonism (more negative) as the THI_a level for W15 increases, while the opposite occurs for W30 and DG (Figure 3a–c). On the other hand, the genetic correlations indicate that the expression of the character cannot be assumed to be the same along the THI_a scale (Figure 3d–f), depicting a clear genotype‐environmental interaction (GxE). Since this differential pattern cannot be detected in Models 1 and 1a, it could bias the estimations of BVs and therefore the results of the breeding program.

To further examine the existence of GxE, we estimated the correlations between BVs according to Models 1a and 2 and the possible variations in the form of expressions along the THI_a trajectory. This was performed using only W30 as an example for the analysed trait. Figures 4 and 5 show the distribution of BV frequencies for direct and maternal effects, respectively, as well as the best 600 animals according to their BVs. The results showed a high correlation (0.945 and 0.956 for BVd and BVm effects, respectively). However, this near‐unity relationship is only applicable when the results are expressed in a transversal form (adjusted means of W30), which does not reflect the nature of the problem for traits expressed longitudinally over the trajectory of a climate–environment gradient. These patterns are presented at the bottom of Figures 4 and 5 with the evolution of BVd and BVm for W30 for all and the best 600 animals but estimated by RN Model 2 along the scale of THI_a 30 days after lambing. In this case, the existence of a double environmental genotype interaction (DIGE) in both genetic effects is clear, involving changes in the order of merit, with a depressive tendency for direct genetic effects as the scale of THI_a increases and the opposite for maternal effects.

FIGURE 4

Frequency distribution of the direct breeding value (BV) effect for 30‐day live weight estimated by two models and the pattern formed by the best 600 animals along the scale of THI level after lambing

FIGURE 5

Frequency distribution of the maternal BV effect for 30‐day live weight estimated by two models and the pattern formed by the best 600 animals along the scale of THI level after lambing

Generalising the PCs and construction of the new index

The manifestation of both fg along the trajectory of THI_a are the BVs represented at the bottom of Figures 4 and 5. PCA was used to analyse these results for W30 of direct (BVd) and maternal effects (BVm) estimated, including 12 expressions of the genetic merit of animals (six for each genetic effect) for THI_a covering an ample range of environmental conditions (THI_a = 45–50; THI_a = 60–65 and THI_a = 65–70). The results show that 56.3%, 27.4% and 16.2% of the total (co)variance of BVd and BVm for W30 over the scale of the selected range of THI_a levels are explained by the first three , which is a considerable reduction in the number of variables. Furthermore, the overall results can be easily interpreted by analysing a biplot of the two first components and the corresponding eigenvectors (ev), in which the complex nature of the DIGE is clearly expressed (Figure 6). For instance, the ev for direct and maternal genetic effects are located in contrasting quadrants, while for each stress zone, the ev of both genetic effects went in opposite directions as a consequence of the negative relationships between them (already shown in Table 4 and Figure 3). On the contrary, ev presents a different trend since its value increased from the cold to the hot zone for genetic direct effects, while it decreased for maternal genetic effects. For those reasons, the existence of such DIGE creates a major obstacle for optimising the Spanish Merino breeding program but also for any breed bred under variable environmental conditions.

FIGURE 6

Eigenvector biplot for direct and maternal BV for live weight at 30 days of age in Merino sheep, estimated by a reaction norm model in different heat stress (HS) zones

To design a procedure to reduce these limitations, the dataset for the total BV (BVt) was subjected to a PCA, which showed that 62.5% and 37.5% of the variability of BVt was explained by the two first λ. Their ev (Table 6) was used as the weighting factor of a new total underlying BV (Ipct), which was described in the previous section. The results showed that all the first ev 1 presented a positive value, suggesting a positive response at all the points of THI_a if a selection process based on this PC eigenvector was carried out. However, the pattern for ev 2 is different, with a contribution in the same direction during the less stressful zone of THI_a ≤ 55 but the opposite during the depressive phase (as shown in Figure 2). However, since both eigenvectors are orthogonal (not correlated), they can be added, and therefore, the PCA‐based Ipct subindex will retain the direction and maximum (co)variance structure of the original BVd and BVm used to estimate BVt.

TABLE 6

Eigenvalue (λ) and eigenvector (ev) for the two leading principal component analyses of the total breeding value (BV) for live weight at 30 days of age

	Weighting factor						Eigenvalues
	BVt 45	BVt 50	BVt 55	BVt 60	BVt 65	BVt 70	λ	λ%
evt1	0.227	0.371	0.485	0.486	0.436	0.388	3.75	62.5
evt2	0.636	0.476	0.161	−0.155	−0.345	−0.447	2.25	37.5

BVt 45 to BVt 70 are the total BV as presented in the Materials and Methods section at THI_a = 45 to 70.

According to the results from Table 6, the new subjacent index will be:

The frequency distribution of Ipct and their relationships with the original BVd and BVm are illustrated in Figure 7a,b. The top part depicted the variability of Ipct, which synthesises the joint action of the total genetic merit for BVt (BVd + BVm) for W30 in all the animals, whereas the Q–Q plot of the joint distribution of the synthetic results of Ipct and the original estimates for BVd and BVm in the three stress zones of the THI_a scale is shown at the bottom. Although Q–Q plot is not a statistical test, it can be employed as an auxiliary graphical representation that allows us to evaluate whether the datasets show any changes in the ranking between the underlying Ipct index and the original BVd and BVm. In this case, the results form a straight line, and the proposed alternative can therefore be used as a simple solution to identify with a single value the best animals according to their breeding merit for direct and maternal effects throughout the range of HS levels represented in this database. Similarly, the best 600 animals (top) and their evolution (bottom) are highlighted in Figure 7b, which also shows the variation in the form of the response of their original BVd and BVm obtained from RN Model 2.

FIGURE 7

Frequency histogram of the index Ipct: (a) percentile ranking (left) for live weight at 30 days of age in three HS zones and (b) form of response (right) for live weight at 30 days of age along with THI at 30 days after lambing (based on results of Model 2)

DISCUSSION

The results of this study show that the clear impact effects of thermal stress in Merino lambs during the first 30 days of life can only be correctly quantified when an additional environmental maternal effect for HS during the last 30 days of pregnancy is included in the model. The lambs born from those mothers pregnant in a hot zone (THI_b > = 65) were between 4% and 8% lower for W15, W30 and DG with respect to those pregnant in a cold zone (THI_b < = 55), but the principal negative effect was on W0 (−17%, see Table 2). This behavior is consistent with the results presented by McCrabb et al. (1993) and Everett‐Hincks and Dodds (2008).

The economic importance of this environmental maternal climatological HS can be inferred if it is considered that 29.8% of the sheep included in this dataset were pregnant during the last 30 days of gestation in a hot zone condition (see Figure 1). In general, for an interval of THI_a ≤ 59 or 60, the pre‐weaning performance in this population of Merinos is more or less stable, but after that level, a depressive phase begins, whose maximum expression is presented at the ends of the THI_a scale. However, despite the clear trend detected in our analysis, the specific values of this range should not be generalised, as there may be variations depending on the conditions of each region or production system or on the particularities of genetic merit of the animals (Kerr, 2015).

The importance of this study can be understood if we consider that these environmental effects of HS are not included in the genetic evaluation models of the Spanish Merino breed (nor in any sheep breeding worldwide, to our knowledge). As a result, the results of the breeding evaluation may be biased. This warning has already been given for milk cattle and beef cattle and sheep. Despite that this source of bias is not related to changes in the order of merit between the two types of evaluation models when the results are expressed in a transversal point (as shown in Figures 4 and 5), the effect becomes evident when the analysis is performed along the longitudinal scale of the THI_a. Unfortunately, these variations in the form of response of the ability of animals to adapt to different levels of HS are not used for BV estimation purposes when using models like 0, 1 or 1a.

Currently, the estimated parameters produced by Models 0 and 1a (Table 2) are consistent with the evidence published in the Spanish Merino breed (Molina et al., 2007), as well as for different sheep breeds across the world (review by Safari et al., 2005). However, when Model 2 was applied (Figure 3), the average values of the heritability estimates for direct () and maternal () effects, as well as the correlation between these (r), were similar to those published by and within the same breed but applying different models. According to Gudex (2001), there is genetic variation for adaptation to HS in sheep, and our heritability results lie within the same range cited by this author.

The negative relationships between the direct and maternal effects were evident in all the models and traits analysed (Table 3 and Figure 3), in agreement with the compelling evidence existing on this subject in different species (Meyer, 1992; Wilson & Réale, 2006). In the context of this analysis, which includes maternal and direct climate effects on the animal, this antagonism may be the consequence of a thermoregulation mechanism presented by Laburn et al. (2002), in which the pregnant sheep protects the fetus in stressful situations of extreme heat or cold by allowing a greater (or lesser) flow of nutrients through the placenta. This ability of the pregnant female does not imply total protection, and in any case, such environmental effects can be fully mitigated due to variations in HS levels during the first weeks of lamb life.

Our results also support the approaches outlined regarding the effects of HS during the gestation period of the female, which is a component that can modify the phenotypic expression of DG, W15 and W30 (Figures 2 and 3), as well as the structure of genetic (co)variances in the post‐natal period (Table 3). According to and, these HS effects should be considered in genetic analysis models; otherwise, they could bias the estimated parameters and affect the selection of breeding programs.

In general, the values of r and r for W15, W30 and DG decrease, and r shows its antagonism (Figure 3) along the post‐natal climate–environment scale. In this sense, the double joint interaction of climatic environmental conditions during the gestation period and the first weeks of life suggest that these traits studied cannot be considered as the expression of the same character along the gradient of climatologic ambient levels, since it changes the ranking and stability of the BVs (Table 4; Figures 4 and 5). Therefore, the novelty of this study lies in identifying the existence of this double genotype‐environment interaction along the trajectory of THI_a for the first time in sheep (Figure 6), demonstrating that the best animals in the thermal stress zone will not be the best in the comfort zone for their direct and maternal genetic effects. No previous references have been found in the analysis of this new type of antagonism between direct and maternal genetic effects and tolerance to various levels of HS in sheep, although a similar pattern was published for weaning weight in several cattle breeds.

Finally, most of the previous studies assessing the antagonism between direct and maternal genetic effects on pre‐weaning growth characteristics have been based on traits whose manifestation is transversal. On the contrary, our results are based on variables expressed longitudinally along an ambient level gradient (known as ‘function value traits’, according to) but also agree with this trend. Therefore, the use of Ipct, which was estimated based on this concept, has been proven to be a valid and reliable alternative tool that can improve the outcome of the breeding program of the Spanish Merino sheep. Despite that the results presented in the bottom part of Figure 7a are expressed in terms of percentiles, it should be noted that this subjacent index Ipct absorbs the maximum of the (co)variances and direction of the BV of each animal for direct and maternal effects along the HS trajectory. Therefore, our results suggested the usefulness of using this type of index agrees with results previously presented in the analysis of beef (Boligon et al., 2016) and dairy cattle (Togashi & Lin, 2006), egg production, and more recently, in different autochthonous Spanish goat breeds.

CONCLUSION

This study demonstrated the existence of a significant HS effect affecting the live weight during the pre‐weaning period in Spanish Merino lambs. However, this HS effect acts together with a previous maternal environmental HS during the last 30 days of pregnancy of the mother. Therefore, both factors need to be taken into consideration in the statistical model to obtain an unbiased estimation of BVs. In addition, we demonstrated the existence of a double interaction between direct and maternal genetic effects and along the trajectory of the effects of HS in the first weeks of life in Spanish Merino sheep reared in extensive conditions. Therefore, it is necessary to review the methodology employed for selecting the best individuals from a model whose results are expressed transversally as a deviation from the mean of environmental conditions to a longitudinal approach that uses a function of the scale of environmental conditions based on a RN model. In addition, we demonstrated the validity and usefulness of applying a PCA within the values for BV direct and BV maternal effects along the scale of direct climate–environment levels on the animal producing the records and combining these with the coefficients of the most relevant eigenvectors as a weighting factor. This allows us to estimate a new underlying index that provides the information needed to identify those genotypes that are best able to manifest good performance in a variety of environmental conditions and maximise the response for both genetic effects over a wide range of climate levels.

CONFLICT OF INTEREST

The authors have no conflicts of interest to declare.

ETHICS STATEMENT

All the data and blood samples for DNA extraction were provided by the Asociación Nacional de Criadores de Ganado Merino (ANCGM). No experiments were performed in the individuals.

AUTHOR CONTRIBUTIONS

A Molina: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing ‐ original draft, Writing ‐ review & editing; Sebastian Demyda Peyrás: Data curation, Formal analysis, Investigation, Methodology, Software, Writing ‐ original draft, Writing ‐ review & editing; Manuel Sanchez: Data curation, Investigation, Resources, Validation, Writing ‐ original draft; Juan Manuel Serradilla: Conceptualization, Investigation; Alberto Menendez‐Buxadera: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Resources, Software, Validation, Visualization, Writing ‐ original draft, Writing ‐ review & editing.