In many developing countries, goats play a vital role in the economy and livelihoods of nomadic farmers (Peacock, Reference Peacock2005). However, their productivity is often low due to several factors, including the genetic potential of the native breeds (Jembere et al., Reference Jembere, Dessie, Rischkowsky, Kebede, Okeyo and Haile2017). Developing breeding programs requires routine genetic evaluation of animals to enable the selection of the best ones to be the parents of the next generation (Scholtens et al., Reference Scholtens, Lopez‐Villalobos, Garrick, Blair, Lehnert and Snell2019). Genetic and environmental factors significantly affect milk yield (MY) and quality in goats (Selvaggi and Dario, Reference Selvaggi and Dario2015), and the knowledge of the genetic and phenotypic aspects of these traits is essential for developing an efficient breeding program. In other words, knowledge of genetic parameters associated with milk traits provides opportunity for more effective selection of breeding animals. This can result in increased milk production and better milk composition.
Estimates of genetic parameters for milk production traits including MY, fat yield (FY) and protein yield (PY) have been reported in several goat breeds such as Spanish Murciano-Granadina goats (Delgado et al., Reference Delgado, Leon, Gama, Lozano, Quiroz and Camacho2006), La Mancha and Oberhasli goats (Garcia-Peniche et al., Reference Garcia-Peniche, Montaldo, Valencia-Posadas, Wiggans, Hubbard, Torres-Vazquez and Shepard2012), Jamunapari goats (Rout et al., Reference Rout, Matika, Kaushik, Dige, Dass and Singh2017), New Zealand dairy goats (Scholtens et al., Reference Scholtens, Lopez‐Villalobos, Garrick, Blair, Lehnert and Snell2019) and Nigerian Dwarf goats (Valencia-Posadas et al., Reference Valencia-Posadas, Lechuga-Arana, Avila-Ramos, Shepard and Montaldo2022). However, there are limited publications on the estimation of genetic and phenotypic parameters of somatic cell score (SCS) in dairy goat breeds (Scholtens et al., Reference Scholtens, Lopez‐Villalobos, Garrick, Blair, Lehnert and Snell2019; Aguirre-Arroyo et al., Reference Aguirre-Arroyo, Castaneda-Bustos, Lechuga-Arana, Gutierrez-Chavez, Angel-Sahagun, Hernandez-Marin and Valencia-Posadas2024).
The Murciano-Granadina goat breed is one of the most important dairy goat breeds in Spain and is distributed extensively in various countries (Miranda-Alejo et al., 2019). In 2015, the private sector imported the Murciano-Granadina goat breed from Spain to the southern region of Iran, which has tropical climatic conditions. This initiative aimed to enhance production efficiency in the native and nomadic goat flocks of the region, primarily maintained under a low-input, low-output production system. The goal was to improve the livelihoods of rural flock holders in the area. To achieve this, purebred Murciano-Granadina does and bucks were either distributed to local flocks or used for crossbreeding with indigenous goat breeds. The estimates of genetic parameters for the reproductive (Mokhtari et al., 2024a) and growth (Mokhtari et al., Reference Mokhtari, Esmailizadeh, Roudbari, Barazandeh, Gutierrez and Mohebbinejad2024b) traits of the Murciano-Granadina goat breed have been recently published. There is no information on the genetic parameters for milk traits in this population of the Murciano-Granadina goat breed. Therefore, this study aimed to estimate the genetic and phenotypic parameters for MY, FY, PY, fat-to-protein ratio (FPR) and somatic cell count (SCC) in the first lactation of the Murciano-Granadina does in Iran.
Materials and methods
The breed and flock management
The Murciano-Granadina goat breed was synthesized in 1975 by mating the Murciana and Granadina goat breeds in the semi-arid areas in southeastern Spain. The main phenotypic characteristics of the Murciano-Granadina goat breed include a straight or sub-concave profile, a medium-sized body with a tendency to lengthen with black or brown uniform coat colour and 77 and 70 cm high at the withers in males and females, respectively (Delgado et al., Reference Delgado, Landi, Barba, Fernández, Gomez, Camacho, Martínez, Navas, Leon, Simoes and Gutierrez2017). The studied flock has been kept under an intensive management system on a private dairy farm. At birth, the newborn kids were weighed and ear-tagged, while information regarding their sex, birth type and the identities of their dams and sires was recorded. Weaning was at about 80 days of age. By considering a mating ratio of 1 buck per 15 does, maiden does were mated to the fertile bucks at about 11 months of age and 25 kg live body weights.
Data, editing protocol and the studied traits
Pedigree information of 33,858 Murciano-Granadina goats originated from 920 sires and 10,769 dams was considered in the present study. The CFC (Coancestry, inbreeding (F) and Contribution) program (Sargolzaei et al., Reference Sargolzaei, Iwaisaki and Colleau2006) was employed to monitor pedigree errors and prepare the pedigree for genetic analyses. The pedigree structure of the studied population of the Murciano-Granadina goats is presented in Table 1. Among the registered animals, those with both parents known, both parents unknown and one parent known comprised 89.23%, 9.91% and 0.86% of all individuals, respectively. Approximately, 65.48% of the animals had no progeny, while 34.52% did have progeny.
Pedigree structure of the population of the Murciano-Granadina goat breed

The original dataset comprised 76,874 test-day lactation records from 7,196 first-parity Murciano-Granadina does, which were collected from 2017 to 2024 at a private dairy farm located in Ghale-Ganj city, in the southern region of Kerman province, Iran. Test-day records of MY, FY, PY and SCC were limited to does for which measurements were taken between 5 and 275 days in milk, with consecutive sampling intervals of approximately 30 days. Does with lower than three test-day records and higher than nine test-day records were excluded from the dataset. The final dataset included 36,958 test-day records of MY, fat percent, protein percent and SCC of 4,632 first-parity Murciana-Granadina does.
The studied traits were total MY, total FY, total PY, average FPR and average SCS. All the studied traits were calculated considering a 275-day lactation period. The centring data method (Ruiz et al., Reference Ruiz, Oregui and Hohenboken2000), also known as Fleischmann’s method, was used for estimating MY, FY and PY according to the following formula:
\begin{equation*}T{Y_i} = \,{P_1}\, \times \,{D_1} + \,\mathop \sum \limits_{i = 2}^k \frac{{{P_i} + {P_{i - 1}}}}{2}\, \times {D_i} + \,{P_{k + 1}}\, \times 15\end{equation*}where TYi refers to the MY, FY and/or PY per lactation, D 1 = interval between kidding and first recording, Pi = yield of the record i, Di = interval between the record i and the record (i + 1) (with i = 1, …, k) and 15 = supposed number of days between the last recording and the dry-off. By aggregating different test-day records into a cumulative value and genetic analysis of the trait by considering a cumulative value per animal, a comprehensive measure of overall milk production performance would be available for an animal.
Does with records on both milk fat and protein percent were considered for calculating FY, PY and FPR. The values of SCC were right-skewed. Therefore, to achieve values which normally distributed, SCC values were transformed into SCS values through the following transformation (Valencia-Posadas et al., Reference Valencia-Posadas, Lechuga-Arana, Avila-Ramos, Shepard and Montaldo2022):
\begin{equation*}SCS = \,lo{g_2}\,\left( {\frac{{SCC}}{{100000}}} \right) + 3\end{equation*}For SCS and FPR, each goat’s lactation record was obtained as the mean of monthly test-day records of up to 275 days. Descriptive statistics for the studied traits are shown in Table 2.
Descriptive statistics for the studied traits in the Murciano-Granadina goat breed

a S.D., standard deviation; C.V., coefficient of variation; Min., minimum; Max., maximum.
b MY, total milk yield; FY, total fat yield; PY, total protein yield; FPR, average fat-to-protein ratio; SCS, average somatic cell score
Statistical model for genetic analysis
The following multivariate animal model was used for estimating genetic and phenotypic parameters for the traits:
where yi represents a vector containing the calculated records for the traits; bi, ai and ei represent vectors of fixed, direct genetic and residual effects for trait i, respectively. Xi and Zi are design matrices that associate the fixed and direct additive genetic effects to the vector of yi, respectively. The fixed effects were kidding season with four levels (spring, summer, autumn and winter), kidding year with eight levels (2017–2024), kidding type with two levels (single and twin) and age at first kidding with two levels (2 and 3 years old). The GLM procedure of SAS 9.4 software (SAS, 2010) was used to determine the fixed effects to be included in the model and to estimate the least squares means for the traits across the levels of fixed effects. The Tukey–Kramer test was used to compare the mean of the traits across the levels of the considered fixed effects.
The direct additive genetic effects were assumed to follow a multivariate normal distribution with a mean vector of null and a (co)variance matrix of G ⊗ A. Here, G and A represent the additive genetic (co)variance matrix and numerator relationship matrices among does, respectively. Additionally, the residual effects were assumed to follow a multivariate normal distribution with a null mean vector and (co)variance matrix of R ⊗ I, where I and R represent the identity and residual (co)variance matrices, respectively. The symbol ⊗ denotes the Kronecker product. Genetic analyses for estimating variance components and genetic parameters were performed by the WOMBAT program (Meyer, Reference Meyer2013).
Results
General considerations
As presented in Table 1, about 6.65% of registered animals were inbred with average inbreeding coefficients of 0.30% and 4.30% in all and inbred animals, respectively. As shown in Table 2, the average values for MY, FY, PY, FPR and SCS were 428.18 kg, 12.55 kg, 8.94 kg, 1.44 and 5.34, respectively. In the present study, FY (68.21%) and FPR (16.25%) exhibited the highest and lowest phenotypic coefficients of variations, respectively.
The influence of non-genetic effects
The least squares means of the investigated traits of the Murciano-Granadina goat across the levels of the considered fixed effects are presented in Table 3. The kidding year and kidding season of does significantly affect all the studied traits (p< 0.01). Kidding type significantly influenced the MY (p < 0.01), but not the other traits (p > 0.05). In the present study, the age of the doe at kidding significantly influenced only MY (p < 0.01), while it had no significant effect on the other traits (p > 0.05).
Least squares means (±S.E.) for the studied traits of the Murciano-Granadina goat breed

x Means with similar letters in each sub-class within a column do not differ from another at p > 0.05, ns: not significant. Means with different letters in each sub-class within a column differ from another at * p < 0.05 and ** p < 0.01.
y MY, total milk yield; FY, total fat yield; PY, total protein yield; FPR, average fat to protein ratio; SCS, average somatic cell score.
The estimates of genetic and phenotypic parameters
The estimates of heritability, genetic and phenotypic correlations for the studied milk traits in the Murciano-Granadina goats are shown in Table 4. All the heritability, genetic and phenotypic correlation estimates were statistically significant (p < 0.05). The heritability estimates for MY, FY, PY, FPR and SCS were relatively low to medium values of 0.15, 0.06, 0.09, 0.10 and 0.26, respectively. The estimates of genetic correlations among the studied traits ranged from –0.38 for FY-FPR to 0.91 for FY-PY, while the phenotypic correlations ranged from –0.05 for MY-FPR to 0.90 for PY-FY. Positive low genetic and phenotypic correlations were found between SCS and the other traits. The genetic correlations of SCS with MY, FY and PY were 0.16, 0.19 and 0.08, respectively. The corresponding phenotypic correlations were 0.03, 0.05 and 0.06, respectively.
Estimates of heritability (on diagonal), genetic (upper diagonal) and phenotypic (lower diagonal) correlations for the studied milk traits in the Murciano-Granadina goat breed

All the heritability, genetic and phenotypic correlation estimates were statistically significant (p < 0.05).
a MY, total milk yield; FY, total fat yield; PY, total protein yield; FPR, average fat to protein ratio; SCS, average somatic cell score.
Discussion
In the present study, a low level of inbreeding was observed in the investigated population of Murciano-Granadina goats. While this level of inbreeding is not currently a cause for concern, it shows some degree of uncontrolled mating within the population. If such uncontrolled mating practices persist, they could adversely influence on the milk performance traits over time because of inbreeding depression. Mean values for the studied traits in the Murciano-Granadina does are within the reported range of previous studies (Scholtens et al., Reference Scholtens, Lopez‐Villalobos, Garrick, Blair, Lehnert and Snell2019; Valencia-Posadas et al., Reference Valencia-Posadas, Lechuga-Arana, Avila-Ramos, Shepard and Montaldo2022). Valencia-Posadas et al. (Reference Valencia-Posadas, Lechuga-Arana, Avila-Ramos, Shepard and Montaldo2022) reported the averages for SCS, 305-day milk, fat and protein in Nigerian Dwarf goats as 2.85, 314.05 kg, 20.94 kg and 14 kg, respectively. In the present study, both FY and PY had higher coefficients of variation than other traits. Within a population, a high coefficient of variation for a trait denotes that there is a large amount of variability in the trait relative to the mean. The higher coefficients of variations for FY and PY in the Murciano-Granadina does indicate substantial variability in these traits. This variation may be attributed to factors such as genetic diversity among animals arising from within- and between-family variations, differences in management practices and varying environmental conditions.
The significant effects of year and season of kidding on MY (Selvaggi and Dario, Reference Selvaggi and Dario2015; Rout et al., Reference Rout, Matika, Kaushik, Dige, Dass and Singh2017; Atoui et al., Reference Atoui, Carabano, Laroussi, Abdennebi, Salem and Najari2024), as well as on FY and PY (Selvaggi and Dario, Reference Selvaggi and Dario2015; Rout et al., Reference Rout, Matika, Kaushik, Dige, Dass and Singh2017), have been well documented. These effects are attributed to variations in climate, food and nutritional quality (Selvaggi and Dario, Reference Selvaggi and Dario2015). Single kidding does have lower least squares mean for MY compared with twin kidding does. The effect of the kidding type can be explained by the fact that litter size changes the endocrines profiles in prepartum does. Specifically, progesterone, estrogen and placental lactogen produced by fetoplacental unit influence mammogenesis and lactogenesis (De Louis et al., Reference De Louis, Djiane, Houdebine and Terqui1980). Litter size effects on MY may be due to the changes in prepartum endocrine profiles caused by the increased foetal and placental mass (Schoknecht et al., Reference Schoknecht, Nobrega, Petterson, Ehrhardt, Slepetis and Bell1991). The significant effects of kidding type on MY in Jonica goat (Selvaggi and Dario, Reference Selvaggi and Dario2015), Jamunapari goat (Rout et al., Reference Rout, Matika, Kaushik, Dige, Dass and Singh2017), Australian commercial dairy goat (Zamuner et al., Reference Zamuner, DiGiacomo, Cameron and Leury2020) and Tunisian local goat (Atoui et al., Reference Atoui, Carabano, Laroussi, Abdennebi, Salem and Najari2024) breeds were reported. They attributed the superiority of multiple-kidding does to single-kidding does in terms of more total MY during a lactation period to the increase in the hormonal activity such as progesterone, placental lactogen and prolactin, which increases the stimulus for the mammary gland to develop during the gestation period. Contrary to us, Almasri et al. (Reference Almasri, Abou-Bakr, MAEM, Kahil, Asaad, Ghoush and Awad2023) reported that the type of birth of kids had no significant effect on the total MY of Syrian Damascus does. Three-year-old Murciana-Granadina does produced significantly higher milk during the first lactation period than the 2-year-old does (Table 3), which may be explained by the maturity stage of the does and the development of the udder’s secretary tissue (Knight and Peaker, Reference Knight and Peaker1982). Scholtens et al. (Reference Scholtens, Lopez‐Villalobos, Garrick, Blair, Lehnert and Snell2019) reported the significant effects of doe age at kidding on the total lactation yields of milk, fat, protein and SCS in New Zealand dairy goats.
The low to medium heritability estimates for all the studied traits implied the higher importance of non-additive genetic and environmental effects on the performance of Murciano-Granadina does in studied population. In other words, additive genetic effects play a modest role in determining these studied traits. Therefore, the potential of direct genetic selection for improving the performance of the does for these traits may be limited.
Delgado et al. (Reference Delgado, Leon, Gama, Lozano, Quiroz and Camacho2006) reported heritability estimates of 0.14, 0.18 and 0.17 for MY, FY and PY in Murciano-Granadina goats, respectively. Valencia-Posadas et al. (Reference Valencia-Posadas, Lechuga-Arana, Avila-Ramos, Shepard and Montaldo2022) estimated heritability values of 0.319, 0.162, 0.155 and 0.102 for average SCS, total MY, FY and PY in Nigerian Dwarf goats, respectively, which were generally higher than the corresponding estimated values in the present study. García-Peniche et al. (Reference Garcia-Peniche, Montaldo, Valencia-Posadas, Wiggans, Hubbard, Torres-Vazquez and Shepard2012) reported high values of 0.48, 0.43 and 0.54 for heritability estimates of MY, FY and PY in La Mancha breed, respectively. They also reported high heritability estimates of 0.61, 0.60 and 0.59 for MY, FY and PY in the Oberhasli breed, respectively.
By applying a meta-analysis study, Jembere et al. (Reference Jembere, Dessie, Rischkowsky, Kebede, Okeyo and Haile2017) reported weighted heritability estimates of 0.33, 0.40, 0.36, 0.52 and 0.54 for lactation MY, FY, PY, fat percent and protein percent in goat, respectively, which were higher than the corresponding estimated values in the present study. Rupp et al. (Reference Rupp, Clement, Piacere, Robert-Granie and Manfredi2011) reported heritability estimates for SCS in Alpine and Saanen goats in France as 0.20 and 0.24, respectively. Mucha et al. (Reference Mucha, Tortereau, Doeschl-Wilson, Rupp and Conington2022) applied the meta-analysis technique and reported values of 0.27, 0.31, 0.30 and 0.21 for pooled heritability estimates of MY, FY, PY and SCS in dairy goats, respectively. Scholtens et al. (Reference Scholtens, Lopez‐Villalobos, Garrick, Blair, Lehnert and Snell2019) reported heritability estimates for MY, FY, PY and SCS as 0.25, 0.24, 0.24 and 0.21 in New Zealand dairy goats, respectively. Aguirre-Arroyo et al. (Reference Aguirre-Arroyo, Castaneda-Bustos, Lechuga-Arana, Gutierrez-Chavez, Angel-Sahagun, Hernandez-Marin and Valencia-Posadas2024) reported heritability estimates of 0.19, 0.24, 0.26, 0.23, 0.11, 0.26 and 0.19 for SCS in U.S. dairy goat breeds of Nubian, Alpine, LaMancha, Saanen, Toggenburg, Nigerian Dwarf and Oberhasli, respectively, which generally were in agreement with the corresponding estimated value in the present study (0.26). The moderate estimated heritability for SCS of the Murciano-Granadina goat in our study suggests potential success in implementing a genetic improvement program. Valencia-Posadas et al. (Reference Valencia-Posadas, Lechuga-Arana, Avila-Ramos, Shepard and Montaldo2022) also reported low and positive genetic correlations for SCS-MY (0.25), SCS-FY (0.18) and SCS-PY (0.23) in Nigerian Dwarf goats. Unlike us, Valencia-Posadas et al. (Reference Valencia-Posadas, Lechuga-Arana, Avila-Ramos, Shepard and Montaldo2022) obtained negative phenotypic correlation estimates of –0.06 for SCS-MY, –0.10 for SCS-FY and –0.08 for SCS-PY in Nigerian Dwarf goats.
Genetic and phenotypic correlation estimates among MY, FY and PY were positive and medium to high, with the corresponding genetic correlations generally higher than the phenotypic ones. Delgado et al. (Reference Delgado, Leon, Gama, Lozano, Quiroz and Camacho2006) reported high and positive genetic and phenotypic correlations (higher than 0.92) among MY, FY and PY in the Murciano-Granadina goat breed. Valencia-Posadas et al. (Reference Valencia-Posadas, Lechuga-Arana, Avila-Ramos, Shepard and Montaldo2022) reported high and positive genetic and phenotypic correlations among MY, FY and PY in Nigerian Dwarf goats. Scholtens et al. (Reference Scholtens, Lopez‐Villalobos, Garrick, Blair, Lehnert and Snell2019) estimated positive and high phenotypic and genetic correlations among MY, FY and PY traits in New Zealand dairy goats, which were 0.81, 0.93 and 0.93 for genetic correlations of MY-FY, MY-PY and FY-PY in New Zealand dairy goats, respectively. The corresponding phenotypic correlations in New Zealand dairy goats were 0.90, 0.96 and 0.92 for MY-FY, MY-PY and FY-PY, respectively.
The high and positive genetic correlation estimates between MY traits, including MY, FY and PY, suggest that selecting for any of these traits may enhance the others. The low positive genetic and phenotypic correlations between SCS and MY traits imply that an increase in SCS could decrease the quality of milk if the selection is based on high-yielding animals; however, these correlations are relatively low which suggests that there would be minimal changes in the SCS by considering these traits simultaneously for selection. In other words, the low genetic correlations of SCS with other traits suggest that SCS may not be a strong predictor of milk production traits, which is consistent with findings in other goat breeds (Rupp et al., Reference Rupp, Clement, Piacere, Robert-Granie and Manfredi2011). The positive but low phenotypic correlations between SCS and MY, FY and PT further emphasize the limited impact of SCS on overall MY and composition. As shown in Table 4, genetic and phenotypic correlations between FPR and MY traits were negative and low to medium estimates with phenotypic correlations lower in absolute values than the genetic correlations. However, genetic and phenotypic correlations between FPR and SCS were positive estimates of 0.20 and 0.04, respectively.
Conclusion
The results of the present study highlight the significant influence of environmental factors on the milk-related traits of Murciano-Granadina does. The low estimates of direct heritability for the examined traits suggest limited genetic variability within the Murciano-Granadina goat population. Consequently, achieving genetic progress for these traits through direct selection may be slow. Therefore, breeding programs should focus on improving environmental conditions to design an effective selection process for this breed while taking controlled mating more strictly into account. The moderate-to-high estimates of genetic and phenotypic correlations among the studied MY traits (MY, FY and PY) indicate that selection for any of these traits could lead to improvements in the other traits. This approach could be beneficial for selecting does based on the total MY of the first-parity Murciano-Granadina does. We found low positive genetic and phenotypic correlations between SCS and MY traits in first-parity Murciano-Granadina does. This suggests that a slight increase in SCS, which serves as an indicator of milk quality, may occur as a result of selecting for higher MY.
Acknowledgements
The authors would like to thank the staff of Ghale-Ganj dairy farm and Fajr Isfahan Agricultural and Livestock Company for their assistance in data collection and flock management.
Funding statement
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Competing interests
The authors certify that there are no conflicts of interest among authors and between authors and other people and organizations.
Author contributions
M.M., A.E. and M.A.F. developed the theoretical framework. M.M., Z.R., A.B. and M.M.S. analysed the data. E.M. was involved in data collection and processing. M.M. wrote the primary version of the manuscript. All authors have contributed to the manuscript revision and have read and approved the submitted version.
Ethical standards
Not applicable.