UT of Jammu and Kashmir have a considerable livestock population (1.79% of the country’s total livestock and 5.15% of its sheep) but, productivity remains lower than global averages. The current mutton shortfall is estimated at around 611.16 million kg, equating to a deficit of approximately 1400 crore rupees. Corriedale sheep were introduced from New Zealand in the 1970s, initially managed in state-owned breeding farms (Khan et al., 2020) and later distributed to neighbouring areas through collaborative breeding initiatives. Ewes exhibit strong maternal characteristics, with a twinning rate from 5% to 25%. Lambs weigh between 3 to 4.5 kg at birth, while adult rams and ewes have a body weight of 79 to 125 kg and 59 to 82 kg respectively (Janbaz et al., 2026). They are valued for early maturity, robust physique, meat production and adaptability to the region’s conditions (Khan et al., 2022).
       
Pre-weaning and post-weaning weight are economically important in sheep production (Janbaz et al., 2026) and exhibit moderate to highly heritability, offering wide opportunities for genetic improvement. Precision and unbiased estimates of genetic parameters are essential for effective selection and breeding programs. Traditionally, growth traits measured at different ages have been analyzed as separate traits and consequently dealt through univariate and multivariate analyses (Fogarty et al., 1994; Safari et al., 2005; Bhatia and Arora, 2005). However, growth being a longitudinal data needs special statistical treatment since their covariance pattern varies among repeated measurements over time.
       
Random regression model (RRM) can estimate variance components and genetic parameters for any point in the trajectory of the growth curve. These models allow estimation of genetic parameters and implicitly models the correlation among successive measures for traits over time, without any constraints imposed by multi-trait discrete models. RRM easily handles irregular measurements, accounts time specific environmental effects and genetic differences in each animal’s growth trajectory.
       
The variation in growth traits over time is less comprehensively measured in sheep industry. The random regression model is extremely valuable to the sheep industry as it allows for the quantification of genetic merit at an infinite number of age points. In sheep industry, utilizing growth curve variation for optimal selection decisions is of immense value so, it becomes relevant to assess the growth trajectory of Corriedale sheep population by Random Regression Model.

Description of data, location and farm management
 
The research was conducted on 5492 lambs born to 101 sires (rams) and 1615 dams (ewes) in different generations and spanned over 37 years on Corriedale sheep maintained at Mountain Sheep and Goat Research Station (MRCS&G), SKUAST-Kashmir, under organized conditions. The farm is situated in the Ganderbal District of J&K, India, with longitude and latitude coordinates of 74.47^oE and 34.14^oN, respectively and an elevation of 5300 ft above the mean sea level. The climate at the farm is temperate, characterized by cool winters (average daytime temperature of 2.5^oC with night temperatures below freezing point) and warm summers (average temperature of 24^oC). The spring season is marked with regular rainfall, while the autumn season tends to be predominantly dry. The animals are reared under a semi-migratory and semi-intensive production system. During the day, they are allowed to graze and are kept in enclosed paddocks at night. The animals from mid-June to mid-September, migrate to highland pastures (Laderwas, Sonamarg, at an altitude of 11,800-14,000 ft above mean sea level) for summer grazing. The breeding season begins in September and lasts until November, with flushing conducted prior to mating season. Ewes are included in breeding plan after reaching 18 months of age. Rams are kept separate from ewes during grazing and tupping is allowed inside paddocks during night only. Following mating, rams are marked with different wool colors on their briskets and ewes receive rump stamps every morning. Close inbreeding is avoided by planned and controlled mating, via housing ewes with specific only and separately rearing male and female flocks during the day.
       
Information regarding lamb Ids and their inheritance (Ram and Ewe records) were collected from history-cum pedigree sheets maintained at the farm. Other information collected included the date of birth, sex of lamb and dam parity. Animals were recorded for birth weight (BW), weaning weight at 3 months (WW), six-month weight (6M), nine-month weight (9M) and twelve-month weight (12M). Digital balance was used for recording the body weight of the lambs. For the final analysis, the growth data were categorized into 9 periods of 4 years with the last period spanning 5 years intervals. The data recorded were divided into two seasons of lambing (Jan-Feb and Mar-Apr), sex of lamb (male and female) and dam parity (1-9). The characteristics for this dataset are illustrated in Table 1.


 
Statistical analysis of data
 
Least squares analysis of variance was used to assess the significance of non-genetic factors (Period and Season of birth, Sex of lamb, Parity of dam) by emmeans package in R. The genetic analysis was done by including the random effects and significant non genetic factors in the Random Regression Model. Single trait RRM was used for analysis, as body weights were recorded at fixed and discrete age points. According to Schaeffer (2004), researcher can decide which classification variables or mathematical functions to use. All RRMs fitted Legendre polynomials of age in months as independent variables. This method has been previously applied in genetic evaluations (Schaeffer, 2004; Skorput et al., 2014; Venkataramanan, 2016; Mahala et al., 2020; Dige et al., 2021). The data was fitted with four sets (Direct additive, Individual permanent environment, Maternal genetic and Maternal permanent environment) of random regression coefficients. Maternal genetic and maternal permanent environmental effects were modeled allowing their contribution to vary across the growth trajectory through age-dependent Legendre polynomials. Direct and maternal genetic effects were assumed to be proportionate to the numerator relationship matrix (A). Further analysis was conducted, considering polynomial (k) fits from 2^nd to 4^th orders of, including a constant term and powers of age up to (k-1) for the four random effects, to identify the best model that describes the data. The general model is represented as follows:

 
Where,
Y_ij: Body weight of i^th animal at j^th month of age.
T: Age in the original scale for which (t * _ij) is calculated as
(t * _ij): Ages standardized between -1 and +1, derived as

 

 

T_min : Earliest date or youngest age;
T_max : Latest date or oldest age.
ϕ_m (t * _ij): m^th Legendre polynomials of age.
X_ij : Set of fixed effects.
β_m : Fixed regression coefficients for modeling the population mean.
α_im, γ_im, δ_im and ρ_im: Random regression coefficients for direct genetic, maternal genetic, direct and maternal permanent environmental effects, respectively.
k_A-1, k_M-1, k_p-1 and k_C-1:  Corresponding order of polynomial fit for each effect.
e_ij: Residual effect.
       
The data were analyzed using the Average Information Restricted Maximum Likelihood (AIREML) algorithm by the WOMBAT software of Meyer (2007). The analysis took into account polynomial (k) fits from 2^nd to 4^th order, which included a constant term and age powers (up to k-1) for four random effects, all aimed at identifying the best model to describe the data. To rank the models, Logarithm of the REML function (LogL), Akaike’s information criterion (AIC) (Akaike, 1998) and Bayesian Information criterion (BIC) (Schwarz, 1978) were calculated.
 

AIC = -2 log L + 2p
BIC = -2 log L + plog (N)

 
Where,
p= No. of parameters.
N= No. of observations.
logL= Maximized log likelihood.
       
AIC and BIC take into account number of parameters and sample size (Foulley and Robert-Granie, 2002). However, Likelihood ratio test (LRT) was used to find the best model (Wolfinger, 1993). The residual variances independently distributed with heterogeneous values varying with age were considered. Changes in error variances were modeled as a step function, dividing age into five distinct classes each spaced out by 3 months. The variations in growth patterns over time were explained using the eigenvalues and their corresponding eigenvectors from the genetic regression coefficients matrix (Kirkpatrick et al., 1990).
       
Eigen function for j^th age and i^th eigenvalue of random regression coefficient matrix was calculated as:
 

ϕ_je_i

 
Where,
e_i= Eigenvector for the i^th eigenvalue.
       
The genetic covariance between ages was calculated as (Jamrozik et al., 1997).
 

G_m = ϕ_mj ϕ′_mj 

 
Where,
G_m=  Covariance matrix for m = additive, maternal genetic, maternal permanent environment and individual permanent environmental effects
ϕ_mj= Legendre polynomials for the random effect of m and j^th age groups.
       
Estimated breeding values EBV = ϕ_j a′_i for j^th point in the growth curve between birth to 12 months of age was calculated as k_i random solutions obtained for i^th animal and Legendre polynomial corresponding to the j^th age group.
where,
ϕ_j = Legendre polynomial for the j^th age group.
a′_i = Transpose of the vector of random solutions for the i^th animal.

The predicted growth curve in Corriedale sheep for 360 days against the observed 5 data points is depicted in Fig 1. Table 1 summarizes the descriptive statistics for the body weight traits of Corriedale sheep. The least-squares mean (LSM) for BW, WW, 6M, 9M and 12M were 3.55±0.04, 11.31±0.27, 17.76±0.30, 18.88±0.17 and 22.97±0.34 Kg, respectively. The non-genetic factors viz. period and season of birth, sex of lambs and parity of ewes were significant on all the traits under study. The body weight measurements at specific time points were utilized as repeated measures for each animal along the growth curve. As age progressed, the number of records declined due to factors such as culling and mortality. The number of animals with at least 3, 4 and 5 records were 3780 (68.1%), 3081 (55.5%) and 2494 (44.9%), respectively. The model comparison criteria are detailed in Table 2 and 3. The R4334 model, with a 4^th-order fit for direct genetic and individual permanent environmental effects, a 3^rd-order fit for maternal genetic and maternal permanent environmental effects and five heterogeneous residual variances, was selected as the best-fit model.






       
The eigenvalues corresponding to the intercept, linear, quadratic and cubic coefficients from the best model are presented in Table 4. The eigenfunctions, derived from the eigenvectors of the genetic variance-covariance matrix (Kirkpatrick et al., 1990) provide valuable insights into the selection effect along the growth trajectory (Lewis and Brotherstone, 2002). The first eigenfunction for the direct additive genetic effect accounted for 91% of the total genetic variation suggesting selection based on this primary growth pattern would lead to improvement across all ages (Fig 2).




       
Variance ratios (heritabilities) for body weights at different ages are shown in Table 5. The heritability estimates derived from RRM for growth at BW, WW, 6M, 9M and 12M were 0.22±0.03, 0.21±0.03, 0.18±0.02, 0.32±0.03, 0.28±0.03, respectively. The body weight correlations were calculated using the respective genetic covariance matrices at different ages (Table 6, 7 and 8). A positive relationship was observed in the estimates of body weights across distinct age classes, except for the birth weight estimates when compared to live weights at 6, 9 and 12 months of age.








       
This study successfully modeled the growth curve of the Corriedale sheep using a random regression model and demonstrated its usefulness for estimating genetic parameters at any point in the growth curve. The model selection process, based on LRT, AIC and BIC identified R4334 as the optimal model.
       
The analysis started with a polynomial order of 2 for all random effects, sequentially increasing up to order 4. Different random regression models had varying polynomial orders, but consistent fixed effects to model random effects. The R4334 model emerged as the best-fit model based on LRT, AIC and BIC (Table 2 and 3). This model featured a polynomial order of 4 for random effects (Animal direct and environment effect) and order 3 for (Maternal genetic and permanent environmental effect) with 5 distinct levels of heterogeneous residual variance components. The likelihood ratio test (LRT) revealed that adding complexity to the models can significantly enhances the model’s ability to explain the data. Further evaluations with more complex models R4444 and R4434 either showed no significant improvement or resulted in a worse fit, confirming that additional complexity beyond RR4334 isn’t justified. Therefore, RR4334 emerges as the most effective model, providing the best balance between fit and complexity without the risk of overfitting.
       
The direct additive genetic effect shows a high and positive intercept of 91%, explaining most of the variance (Fig 2). Our findings are consistent with (Fischer et al., 2004; Ghafouri-Kesbi et al., 2008; Kariuki et al., 2010, Barazandeh et al., 2012, Mohammadi and Farhadian, 2017, Ghiasi and Mokhtari, 2018; Mahala et al., 2020). The linear coefficient linked to additive genetic, maternal genetic, maternal permanent environmental and individual permanent environmental accounted for 5.65, 1.42, 0.51 and 11.07 per cent of the total variation, respectively.
       
The first eigen function of eigenvalue accounted for 91% of the additive genetic variation across different ages. The trajectory for this function indicated positive and upward trend as age increased. This suggests that focusing on this aspect of variation provides scope for overall improvement at all ages. Lewis and Brotherstone (2002) found similar findings regarding the first eigen function in Suffolk sheep. The 2^nd and 3^rd eigenvalues for the direct additive genetic effect were significantly lower, indicating a limited potential for genetic improvement in animal growth through selection based on these eigen functions. The second Eigen function was positive until about three months and negative thereafter. However, Venkataramanan (2016) reported a consistent trend with positive values and a slight increase in the trajectories for both the first and second Eigenvalues. The 3^rd eigenfunction showed a decreasing trend up to 9 months and thereafter linearly positive up to twelve months. The 4^th eigenfunction started with a negative value, became positive at weaning, thereafter showed decreasing trend until nine months and finally increased linearly up to twelve months of age. Both the third and fourth eigenfunctions didn’t exhibit much variability, meaning they contributed very little to the overall variation in the trait. The polynomial orders for RRM were established based on recommendations from earlier studies on growth (Meyer, 2004; Fischer et al., 2004; Molina et al., 2007; Ghafouri-Kesbi et al., 2008; Mahala et al., 2020). Beyond this point, using higher-order polynomials led to computational issues with convergence, likely due to the limited number of records available at older ages, as noted by Arango et al. (2004).
       
In this study, the highest estimates of direct heritability were found at 09 months of age. In contrast, the heritability trends were not consistent with previous findings (Fischer et al., 2004; Venkataramanan, 2016), which reported a steady increase in direct heritability up to 12 months of age. The high heritability at 9M indicates that this age is under strong genetic control and is an ideal time to make a selection decision in this population. The variance attributed to maternal genetic and maternal permanent environmental factors throughout the growth trajectory were quite low, suggesting that while maternal influence on birth weight was present, it diminished after weaning in RRM. Several researchers (Molina et al., 2007; Ghafouri-Kesbi et al., 2008; Mahala et al., 2020) have also noted a decline in maternal genetic and permanent environmental effects as the animals aged. It’s worth mentioning that some studies have reported higher estimates than ours, likely due to differences in breed (Barazandeh et al., 2012; Kheirabadi and Rashidi, 2016; Venkataramanan, 2016). Nonetheless, our findings indicate that after weaning, an animal’s genotype plays a more significant role in determining body weight compared to other random effects.
       
The individual permanent environmental variance contributed the most to the variance in growth traits with highest estimate at nine months of age. This could be attributed to the potential overestimation of additive effects in later ages, possibly due to end-effect of polynomials in Random Regression Model (Ghafouri-Kesbi et al., 2008).
       
Consistent with our findings (Fischer et al., 2004; Kariuki et al., 2010; Barazandeh et al., 2012; Mahala et al., 2020) showed similar trends in phenotypic correlation, The genetic correlations between the youngest and oldest were generally low to moderate, while it was higher between subsequent ages. This pattern aligns with previous studies on genetic correlation for growth data using RRM (Lewis and Brotherstone, 2002; Fischer et al., 2004; Aziz et al., 2005; Ghafouri-Kesbi et al., 2008; Barazandeh et al., 2012; Venkataramanan, 2016; Mahala et al., 2020). The findings of our study support earlier findings (Ghafouri-Kesbi et al., 2008; Venkataramanan, 2016), suggesting that maternal effects at a younger age are influenced by different genes than those expressed at a later age.

The growth trajectory was effectively modeled using the Random Regression Model, which enables estimation of variance components across the entire growth curve. Likelihood ratio test, Bayesian information and Akaike information criteria suggested R4334 as the optimal Random Regression Model. Random Regression Models are highly suitable for genetic evaluation of growth traits under field recording systems with unevenly distributed data across all age groups.

The authors are indebted to Hon’ble Vice Chancellor SKUAST-Kashmir for providing the necessary infrastructure facilities and to SKUAST-K start up (Kashvet Innovations Private Limited) and ANRF-ARG (ANRF/ARG/2025/009380/LS) for the financial support that enabled the successful completion of the project.

The authors declare that no actual or potential conflict of interest could inappropriately influence the work.