The Impact of Age at First Lambing on Milk Yield and Lactation Length in a Population of Istrian Sheep under Semi-Intensive Management

Simple Summary The impact of age at first lambing on dairy traits has been poorly investigated in sheep, especially in semi-intensive husbandry systems. Insufficient information on this issue, especially scientifically proven, leads to many speculations among breeders. The aim of this study was to examine the impact of age at first lambing on days in milk, daily milk yield, and total milk yield. The study was conducted on field data routinely collected in a population of Istrian sheep bred under selection. It was determined that the prolongation of the first mating of lamb ewes to the second year of life is not beneficial to milk production. In conclusion, ewes reared under semi-intensive dairy orientated systems could be successfully bred in their first year of life in order to reduce their unproductive non-milking phase. Abstract This study aimed to examine the impact of ewe’s age at first lambing (AFL) on days in milk (DIM), average daily milk yield (DMY), and total milk yield (TMY). Symmetrical bimodal distribution of AFL enabled classification of maidens in those mated in the first (47%) or second year of life (53%). After accounting for all available sources of phenotypic variability with the linear mixed model for repeated records, it was estimated that AFL had a statistically significant effect only on DIM (p < 0.001). The litter size had a significant effect only on TMY (p < 0.001), while the effect of the parity was significant for all the examined traits (p < 0.001). The results of the study suggest that prolongation of age at first mating to the second year of life is not justified in dairy-orientated sheep farms. However, more evidence on this issue is needed for generalization, especially considering some other traits that can impact profitability of dual-purpose sheep farms (reproduction traits, growth rate of lambs, etc.).


Introduction
Dairy-orientated sheep husbandry systems, dominantly present in the Mediterranean basin, mainly rely on local breeds well-adapted to the specific environment and seasonal availability of pastures [1]. The majority of systems "tagged" as a dairy are actually a dual purpose, with a substantial portion of incomes provided from lamb's meat [2]. In France, Greece, Italy, and Spain (representatives of the dairy sheep sector), the average milk yield varies from low to medium (85-216 L/ewe/lactation) [3]. A similar range of milk yields is present in Croatia where dairy sheep farming relies on three foreign (Lacaune, East-Friesian, and Travnik pramenka) and two autochthonous breeds (Pag sheep and Istrian sheep).
The Istrian sheep is a dual-purpose breed originating from the Istrian peninsula whose formation officially started in 1771. Many foreign breeds (rams) were involved in creation 2 of 10 of the breed: Gentile di Puglia, Bergamasca, Southdown, Merinolandschaf, Merino, Awassi, and East-Friesian [4,5]. In total, 1632 individuals are currently included in the national selection program [6]. BLUP has been used as selection tool for about a decade now [7] with a test-day repeatability animal model [8]. The transition to a single-step genomic BLUP is underway.
Planning of the breeding in flocks, especially in seasonally sensitive systems, presents a challenge for farmers since both economical and biological considerations need to be carefully balanced [9,10]. The optimum age at first lambing in dairy sheep heavily relies on husbandry system, local dairy facilities policies, inherent characteristics of the breed, etc. However, with all the information available, both from the practice and science, it is still somewhat difficult to determine the optimal age for lamb ewes to be bred. The breeding of ewe lambs as soon as possible without repercussions on their productivity is an efficient way of making dairy sheep operations more productive by reducing the unproductive phase of the females. On the other hand, breeding lamb ewes at a very young age poses risk for incomplete ewe's body development and soundness.
Conversely to the abundant information of optimal age at first rearing in beef [11][12][13] and dairy [14][15][16][17][18][19] cattle, information on this issue for dairy sheep is limited. Very few systematically conducted studies on this have been published so far for dairy sheep, and none of them was related to semi-intensive (pasture-based) rearing. Gootwine and Pollott [20,21] found that breeding Awassi and Assaf ewe lambs later in life led to increased milk production in their first lactation, but they concluded that it was not economically justified. Hernandez et al. [22] in a study conducted on Lacaune breed reported the optimal age at first lambing is 13-15 months. However, both of the above-mentioned studies were conducted in intensive management systems. Absence of scientific evidence leads to many speculations on this issue, thus the aim of this study was to examine the effect of age at first lambing on some dairy characteristics of Istrian sheep reared in a semi-intensive (pasture-based) breeding environment.

Animals
The study was conducted using field phenotypic records collected on 1891 Istrian ewes born between 2009 and 2018 in 101 flocks. The ewes under study had been performance tested for the purpose of the national selection program on dairy traits. All the data used in the study were provided by the Croatian Ministry of Agriculture. The ewes had been raised in a semi-intensive husbandry system. Natural pastures had been the main source of food all year round, except in the couple of cold winter months (depending on the year) when hay represented the major component of the rations (1.5-2 kg/day). In addition to forages (pasture and hay), ewes were fed with concentrates (corn or corn/barley) in the last two months of gestation and during lactation (mainly while being milked or after milking, usually 200-500 g/day).

Data
In total, 4725 records for each trait (days in milk (DIM), total milk yield in the milking period (TMY), and average daily milk yield (DMY)) were used in the analysis. Testday records were collected monthly in accordance with the regular alternate scheme (morning/evening system) based on ICAR rules [23]. The TMY per lactation was estimated using the Fleischmann's method [24] from test-day records for all ewes that had at least three milk recordings within lactation. The first recording was conducted between the 5th and 30th day after weaning and thereafter for approximately 30 days (28-34) until the end of lactation (i.e., when daily milk secretion fell below 0.2 kg). The age at first lambing (AFL) was available for all ewes under study. The distribution of AFL had two peaks centered at 14 and 24 months of age (  Prior to the estimation of the effect of AFL on the examined traits, a thorough analysis of distribution of phenotypic records by parity, litter size, and lambing season was conducted. Several adjustments of the original data were performed in order to retain all available data in the analysis as follows: multiple births were set to Class 2+, parities after 7th were set to Cass 7+, October lambings were set to November lambings, while April lambings were set to March lambings. The frequency distribution of phenotypes by parity, litter size, and lambing season after these adjustments is shown in Figure 2.  Prior to the estimation of the effect of AFL on the examined traits, a thorough analysis of distribution of phenotypic records by parity, litter size, and lambing season was conducted. Several adjustments of the original data were performed in order to retain all available data in the analysis as follows: multiple births were set to Class 2+, parities after 7th were set to Cass 7+, October lambings were set to November lambings, while April lambings were set to March lambings. The frequency distribution of phenotypes by parity, litter size, and lambing season after these adjustments is shown in Figure 2.  Prior to the estimation of the effect of AFL on the examined traits, a thorough analysis of distribution of phenotypic records by parity, litter size, and lambing season was conducted. Several adjustments of the original data were performed in order to retain all available data in the analysis as follows: multiple births were set to Class 2+, parities after 7th were set to Cass 7+, October lambings were set to November lambings, while April lambings were set to March lambings. The frequency distribution of phenotypes by parity, litter size, and lambing season after these adjustments is shown in Figure 2.

Statistical Analysis
The inferential statistical analysis is suited to the unbalanced measurement experimental design with repeated measurements. Linear mixed models (LMM) represent an extension of simple linear models to allow accounting for both fixed and random effects. These models are particularly useful when there is non-independence in the data [25]. Accounting for random variability in the model was performed to correct the standard errors for non-independence of the data and not to explore variability within the classes of random variables. Prior to the analysis, the significance level α was set to 0.05. After testing for the significance of the available fixed effects and comparing performances of several different statistical models, the following mixed linear model (Equation (1)) was used in the analysis of the total milk yield (TMY) and daily milk yield (DMY): where Y ijklmn is the nth observation and µ is the overall mean in groups ijklm. Parity (P i , I = 1, . . . , 7+), litter size (LS j , j = 1, 2+), and age at first lambing (AFL k , k = 1, 2) were fitted as class, while length of suckling (s) and milking (m) period as continuous fixed predictors. The herd-year-month of lambing (hym) and animal (a) were fitted as random effects. The effects length of suckling (s) and length of milking (m) period were fitted as fixed continuous predictors (linear regressions). It might seem at first that one of these covariates is redundant in the model, but it turned out that both covariates were informative and, moreover, practically uncorrelated (see the paper by Kasap et al. [26] for more details). The reduced statistical model (Equation (2)), without covariates s and m, was used in the analysis of the days in milk (DIM).
Separate random effects terms were considered independent. In both models, the covariance structure was unstructured. We did not impose any constraints to the values, so each variance and each covariance were estimated uniquely from the data which provided the best possible fit (very close to what the data reflected).
The statistical analysis was conducted in R programing environment [27]. Figures were obtained with package "ggplot2" [28]. The package "lme4" [29] was used in estimation of the AFL on the examined traits. The R function "lmer()" is described in detail in [30]. Estimation of the least square means (LSM) and their standard errors was obtained with the package "emmeans" [31].

Descriptive Statistics
Distributions of the analyzed traits classified by litter size, parity, and AFL are shown in Figure 3. On average, multiple-bearing ewes had longer lactation, produced more milk in the milking period, and had higher DMY. DIM consistently increased with parity while TMY and DMY only up to the third parity and then decreased. Ewes mated for the first time in the second year of life (AFL2) had higher DIM but lower DMY and TMY. AFL2 showed greater homogeneity in phenotypic expression for all the examined traits.

Inferential Statistics
The results of the inferential statistical analysis (Table 1) are generally in line with the raw means discussed above (descriptive statistics). The estimated marginal means (LSM) presented in Table 1 were averaged across the levels of the other fixed class effects for the sake of simplicity and better insight into the magnitude of the particular effect. The estimates (LSM) obtained from the same statistical model, but simultaneously taking into account all the fixed class effects, without averaging on the other effects, are presented in Table 2. These results enable comparisons of the estimates across all possible levels of class fixed effects, and they correspond to the summarized estimates in Table 1. Description of the results in this form gives the opportunity to easily get information on

Inferential Statistics
The results of the inferential statistical analysis (Table 1) are generally in line with the raw means discussed above (descriptive statistics). The estimated marginal means (LSM) presented in Table 1 were averaged across the levels of the other fixed class effects for the sake of simplicity and better insight into the magnitude of the particular effect. The estimates (LSM) obtained from the same statistical model, but simultaneously taking into account all the fixed class effects, without averaging on the other effects, are presented in Table 2. These results enable comparisons of the estimates across all possible levels of class fixed effects, and they correspond to the summarized estimates in Table 1. Description of the results in this form gives the opportunity to easily get information on expected performances of Istrian sheep by simultaneously considering these three sources of phenotypic variability.

Age at First Lambing
The effect of the main interest in this study, AFL, significantly affected only DIM (p < 0.001). It was estimated that lactation lasted 14 days longer in AFL1, but without positive impact of AFL on milk yield (both TMY and DMY). In fact, it was estimated that early mated lamb ewes (AFL1) tend to produce negligibly more milk than AFL2. Close inspection of the results (Table 2) revealed that AFL had a negative effect on the milk yield traits (DMY and TMY) and a positive effect on the DIM irrespective of the parity and litter size.

Age at First Lambing
The effect of the main interest in this study, AFL, significantly affected only DIM (p < 0.001). It was estimated that lactation lasted 14 days longer in AFL1, but without positive impact of AFL on milk yield (both TMY and DMY). In fact, it was estimated that early mated lamb ewes (AFL1) tend to produce negligibly more milk than AFL2. Close inspection of the results (Table 2) revealed that AFL had a negative effect on the milk yield traits (DMY and TMY) and a positive effect on the DIM irrespective of the parity and litter size.

Litter Size
The litter size had a positive effect on DMY and TMY and a negative effect on DIM. However, the effect was statistically significant only for TMY (p < 0.001). Multiple bearing ewes produced~10 L milk more in lactation. As determined above for AFL, the estimated effect of parity was consistent across other class fixed effects in the model ( Table 2).

Parity
Parity had a significant effect on all the traits in the analysis (p < 0.001). Estimates for milk yield (DMY and TMY) reached a maximum in third parity and thereafter decreased, while estimated for DIM increased consistently up to the last class of parity (7 + ). The maximal difference in adjacent estimates was determined between fourth and fifth parity for DMY and TMY and between first and second parity for DIM.

Discussion
Semi-intensive dairy orientated sheep breeding systems dominate over intensive dairy systems in many Mediterranean countries, where the majority of European sheep milk is produced. In such conditions, particularly welcomed efforts are those that tend to increase productive potential of the animals with little or no direct investments. The scientific knowledge on the impact of AFL on ewe's dairy performance is scarce, and there are still many unanswered questions and doubts on this issue, not only among breeders but also among scientists. The question of optimum AFL has been comprehensively addressed by some studies (e.g., [21,22]), but investigated only under intensive management systems where undesired environmental effects can easily be alleviated. However, ewe's performance in conditions where ewes dominantly graze and receive limited amounts of concentrate (strong impact of natural environmental conditions) might have different outcome. Estimation of this effect is pretty challenging in many non-intensively managed sheep populations, mainly due to a lack of available phenotypic data and all other necessary information on systematic effects. Fortunately, this breed has been in the performance recording system for dairy traits for more than a decade. The availability of a large number of systematically collected phenotypes and other important information gave us the opportunity to get some answers with the appropriate statistical model. The population of Istrian sheep under study exhibited a wide range of AFL (from 10 to 29 months), although with a bimodal distribution with peaks at 14 (47%) and 24 (53%) months. The average DIM was 183 days (52 days of suckling and 131 days of milking period). The average DMY and TMY yield in milking period was 1.08 L (sd = 0.50) and 142.5 L (sd = 80.44), respectively. Classification of the observations by parity and litter size was in accordance with our expectations (based on our own experience and many reported results on this issue from the literature). On the other hand, the results pertaining to the AFL, which was the main goal of the study, only partially match with the previous findings on this issue. The lamb ewes mated in their first year of life (AFL1) had a two weeks shorter lactation period, but DMY and TMY were only marginally affected by the AFL. Actually, it turned out that ewes mated earlier in life (AFL1) had negligibly higher DMY and TMY compared to ewes mated in their second year of life (irrespective of the parity). However, the difference between the groups of AFL was not statistically significant for the both of the analyzed traits. Contrary to these results, Pollott and Gootwine [21] in their study on Assaf breed with the AFL range between 10 Animals 2021, 11, 1604 8 of 10 and 28 months found that AFL did not affect DIM, but it considerably affected TMY. They reported that older ewe lambs at the first parity produced more milk in the first lactation, but less milk in second lactation (the estimated regression coefficients of TMY on AFP were 3.7 L/month and −1.82 L/month, respectively). Hernandez et al. [22] reported similar results for Lacaune breed. They also found a positive relationship between AFL and the first lactation TMY and negative relationship between AFL and the second lactation TMY. In addition to the plenty of important findings from their study, it is also worth mentioning here that ewes mated for the first time later in life reached their maximum TMY in their first lactation. This is in accordance with the previous reports for Lacaune [32], as well as for some other sheep breeds (e.g., Latxa, [33]). The authors argued that lambing at younger age probably had lower first lactation yields because of directing greater amounts of energy towards growth during their first lactation compared to ewes lambing at older age when they reached a mature size. Our results are not completely consistent with those reported above because our lamb ewes mated later in life (AFL2) did not exhibit better performance in the first lactation (TMY and DMY). In addition, regardless of the AFL, the ewes at the first parity did not outperform those in later parities. However, the results pertaining to the impact of the AFL on the later lactations correspond with their reports to substantial extent. Our experience and the results of this study suggest that breeding lamb ewes early in life (AFL1) has no detrimental impact on the milk yield, neither in the first nor in the later parities. Based on these results, as well as those reported earlier, we hereby argue that prolongation of the mating (until the second year of life) is not justified by itself if there are no particular detectable obstacles (e.g., poor body development, poor condition, or some kind of unsolved health problems). The authors above also concluded their work on this issue in the same manner because it was found that the milk gains in the first lactation accompanied by increased AFL could not compensate losses in the future lactations and shorter lifetime production. We hereby argue that the results of this study nicely reflect the impact of the AFL on DIM, TMY, and DMY in the population of Istrian sheep reared in a semi-intensive husbandry system. The results should serve as scientific evidence for future breeding activities in this population and some other populations reared in the Mediterranean semi-intensive dairy orientated system.

Conclusions
The results of this study suggest that breeding maiden ewes in the first year of life does not diminish their lifetime dairy performances. To reduce the unproductive phase of the animals, breeding should not be prolonged to another year if soundness, body development, and condition are not compromised.
Although the applied statistical approach partially accounted for the additive genetic effects in revealing the direction and magnitude of the target effect (AFL), more evidence is needed for generalization.
Analysis under the framework of the genetic animal model could provide more evidence on this issue, particularly in the context of the single step genomic BLUP. Extension of this research is advised to take into account some other important traits that can affect the overall success dairy-orientated sheep facilities, such as fertility, growth rate of lambs, etc.