Non-parametric analysis of the effects of nongenetic factors on milk yield, fat, protein, lactose, dry matter content and somatic cell count in Murciano-Granadina goats

Abstract

The objective of this work was to evaluate the influence of non-genetic factors on milk yield, its composition (fat, protein, dry matter and lactose) and the count of somatic cells in Murciano-Granadina goat breed, through the application of alternative nonparametric tests to routine parametric tests to explain the variability found in the population with respect the aforementioned traits. 2594 milk yield and composition records belonging to 159 goats from the selection nucleus were analysed. Predictors evaluated were farm, type of birth and parity order, live kids, parturition month, season and year, control number, control type, control month, season and year, days in lactation, days from first control, days from last control to drying, drying month, season and year. All of them presented a significant influence (p < .0001) on all the variables studied with the exception of the number of stillborn kids which did not significantly influence the percentage of each component and the year of drying which seemed not to significantly influence dry matter percentage. Conclusively, nongenetic factors affect milk yield and its compositions in the Murciano-Granadina goat breed. Additionally, the inclusion of the type of milk control and information related to the drying period in models predicting for milk yield and content may provide interesting information, which must be included in genetic evaluations to promote higher and better-quality milk production improving the profitability of autochthonous breeds.

HIGHLIGHTS

• Non-genetic factors may affect milk components more than milk yield.

• Including factors related to lactation cycle can help identifying critical points.

• Drying-off period information promotes milk yield and quality.

• Studying nongenetic factors helps maximising milk predictive model potential.

• Factor combinations studied explain up to 41.8% of variability of milk yield.

Keywords:

• Dairy goats
• environmental predictors
• non-parametric analysis
• local breed
• milk yield and composition

Introduction

Murciano-Granadina goat breed occupies a remarkable position among international most productively interesting dairy goat breeds. The commercial potentialities of this breed give proof of the wide spreading process that the population has experienced in recent decades to almost any corner in the world from the Mediterranean basin, North Africa to Iberoamerican countries (Martínez et al. 2011). Despite average milk yields fall around 600–700 litres in 210 lactation days, it is more than possible for individuals to reach 1000–1300 milk litres in 210 lactation days (Delgado et al. 2017).

During the past decades, worldwide dairy goat milk and cheese yields suffered a dramatic increase deriving from the significant boost experienced in the scope of market needs. In this context, producers’ cooperative movements needed to evolve to fulfil such increased demands. As a result, breeding programmes started focussing on the promotion of milk productivity, quality and in turn its profitability.

These newly aimed breeding programmes not only sought the obtention of more productively efficient individuals, but also of those animals able to provide better quality dairy products (Deroide et al. 2016).

The selection of individuals is simultaneous to the need for the evaluation and quantification of the effects of genetic and environmental or non-genetic factors, and their interactions as such analyses can help detecting the confounding nature of certain predictors or their combinations, which helps isolating and accurately determining the degree at which indeed predictors condition economically important traits (Pizarro et al. 2019c).

Idowu and Adewumi (2017), reported genetic factors account for a high fraction of interindividual, intrabreed and interbreed variability. This heterogeneity sets the foundation for the improvement of milk yield and composition. Contrastingly, non-genetic or environmental factors condition milk yield at quali-quantitative levels, which in turn derives in the alteration of the quality of dairy products. Consequently, controlling such environmental factors plays an essential role, provided they configure the background which either permits or prevents animals from expressing their complete genetic potential (genetic value) for certain economically relevant traits (Pizarro et al. 2019c).

Studies related to milk production and its components routinely use parametric tests for data analysis (Bagnicka et al. 2015). However, they do not take into account that, in a highly selected population, the selection pressure promotes the erosion of the left tail of the distribution, since low-yielding goats are discarded even before being considered at the database. This selection procedure, alters the properties of data distribution from the beginning. Based on this selective context, we could address the imbalance in favour of variability between farms compared to variability within the farm as one of the reasons for the distortion of distribution properties (Pizarro et al. 2019c).

Contextually, these events may lead to situations in which the compulsory assumptions for running parametric tests are violated (nonnormality, heteroscedasticity, sphericity, multicollinearity and the existence of outliers), which makes nonparametric tests the most preferable option for the study of data.

As a first aim, we assessed non-genetic predictors such as farm, birth number, type of parturition, live kids, stillborn kids, month, season and year of parturition, control number, control type, control month, season and year, days in milk, days from first control, days from last control to drying, drying month, season and year on milk yield and composition and somatic cell count (SCC × 10³ sc/mL) in Murciano-Granadina goat breed. Nonparametric tests were applied as a statistical alternative to parametric tests to identify which predictors could be used to issue specific regression equations to linearly predict milk yield, fat, protein, lactose, and dry matter content (in g), and SCC × 10³ sc/mL (quantified in somatic cells per millilitre, sc/mL).

Material and methods

Animal sample

The present research was designed within the scope of a project aiming at assessing the separate effect of nongenetic predictors from genetic biomarkers to develop predictive models that provide deeper insights on environmental and genetic regulation of milk yield, milk composition and SCC × 10³ sc/mL.

The present study is complementary to two studies focussing on the assessment of the genetic fraction (additive and dominance (Pizarro et al. 2019c, 2020b) and epistatic effects (Pizarro et al. 2020a, 2020b)) of 48 SNPs in the goat casein gene complex, which regulates for the expression of milk yield, fat, protein, dry matter, lactose and SCC × 10³ sc/mL.

Contextually, the study of genetic factors may require making large investments provided the costs involved in genotyping and sequencing. These costs make reducing sample sizes to representative dimensions a routine practice. However, as a result of this sample selection process, the statistical techniques which are normally applied, may not fit the distribution properties and other features of the phenotype datasets from selected samples. In these regards, special efforts must be made on determining highly efficient tools, which may save the drawbacks derived from the power or significance issues of statistical nonparametric alternatives.

To this aim, the sample selection process followed in this study was applied on the animals (n = 2359479 milk yield and components records from 151997 goats) present in the herdbook of the National Association of Breeders of Murciano-Granadina breed goats (CAPRIGRAN). All these registered animals were ranked considering the most recent and updated official breeding values for milk yield and composition reported for all the animals published in 2015.

The overall performance for milk yield and content was evaluated to rank the animals in the productive records to implement sample selection. To this aim, as more than one trait was included, combined selection index (ICO) procedures were applied to summarise the value of each individual comprising its partial values for milk yield and composition following the premises in Van Vleck (1993). We decided not to include dry matter in the ICO, as redundancies may occur deriving from the relationship of this trait and fat or protein content. To determine the weights to apply to each trait, we considered the phenotypic relationship across milk yield and composition traits (except for dry matter) scoring their relevance as selection criteria when the breeding goal was milk yield and quality. In matrix notation, the weights to be applied on the selection index combining the partial scores of each modality were obtained as, $\mathrm{b}={\mathrm{P}}^{-1}\mathrm{g},$ where b is the vector of weights to be applied to each production or content trait, P is the phenotypic (co) variance matrix, and g is the vector of genetic (co)variances of every trait with each other. MatLab r2015a (The MathWorks Inc., 2015) was used to compute all selection index combinations. With the aim to prevent a potential distortion of the effects of the factors to be evaluated in this study and considering the market demands for the simultaneous goals of quantity efficiency and quality efficiency (Mu et al. 2019), the weights for milk yield, fat, protein and lactose were considered equal and followed the proportion of 1: 1: 1: 1, respectively. The combined index used (ICO), could be represented following the equation: $$\mathit{ICO}=\frac{{\mathit{PBV}}_{\mathit{milkyield}}{W}_1}{{\mu }_{\mathit{milkyield}}}+\frac{{\mathit{PBV}}_{fa\mathrm{t}}{W}_2}{{\mu }_{\mathit{fat}}}+\frac{{\mathit{PBV}}_{\mathit{protein}}{W}_3}{{\mu }_{\mathit{protein}}}+\frac{{\mathit{PBV}}_{\mathit{lactose}}{W}_4}{{\mu }_{\mathit{lactose}}}$$ where PBV is the predicted breeding value for each of the traits and animals included in the matrix; W₁ is the weight for milk yield, W₂ for fat, W₃ for protein and W₄ for lactose in g and normalised at 210 days; and $\mu$ the mean for each of the traits included in the ICO computed in g and at 210 days. After ICO was computed for each of the animals included in the matrix, we sorted a total of 200 animals from the whole routine milk recording of Murciano-Granadina goat breed in a ranking considering their ICO value obtained at the previous genetic evaluation. Animals with extreme PBVs may be less efficient and less balanced than we could expect at first. Furthermore, not all traits are affected to the same degree by selection for these extremes. For these reasons, 200 animals were ranked as follows: we chose 67 females presenting the lowest ICO values in the rank, 66 females with values around percentile 50, and the 67 females presenting the highest ICO values in the rank, so as to perform an adjusted representative sampling of the aggregate genotype distribution in the population. The samples belonging to animals with missing or incomplete phenotype registries were discarded, hence the final set for genotyping consisted of blood samples from 159 stud book registered goats out of the 200 animals that were initially considered.

Sample information and management

The animals belonged to twenty-eight different farms in southern Spain. Phenotypical records for milk yield and components were registered randomly during different periods, from 2005 to 2018. Age of individuals ranged from 1.24 years to 9.15 years.

Milk yield and composition standardisation

Murciano-Granadina production system characterises by two annual kidding seasons, with milking periods ranging from 210 to 240 days (Delgado et al. 2017). Total milk yield (in g), composition (in %) and SCC × 10³ (sc/mL) were estimated until 210 days of lactation adapting the procedures described in Pizarro et al. (2019a). Then, percentages for milk protein, fat, dry matter and lactose were transformed into g.

Real yield (RY_j) or ${RY}_j=d_{1j}{Y}_{1j}+30\mathrm{ }{\sum }_{i=2}{Y}_{ij}+\left[d_{2j}-30(n_j-2)\right]{Yn}_j$ was used to calculate individual milk yield and components. RY_j is the real yield of individual (goat) j; Y_1j is milk yield at first control; n is total control number; Y_ij is milk yield at control ij, Yn_j is milk yield at the last control.

Approximately five milk controls were performed per individual on average. Farm internal policies for the implementation of Official quality controls depended on the different milking systems described by the Royal Decree Law 368/2005, on 8 April 2005, which regulates the official control of the milk yield for the genetic evaluation in the bovine, ovine and caprine species of the Spanish Ministry of Agriculture (2005), except for the first control and the last which were assessed per goat individually. Then, days between kidding date (KD) and the date when the first control took place after that kidding (FC) were calculated using, d₁ = FC-KD and the days between the penultimate (PC) and the last control (LC) as follows, d₂ = LC-PC.

Each goat parturition date and consecutive control date records until 210 days of lactation had been completed were considered to normalise milk yield per goat to prevent interindividual differences that may be ascribed to differences in lactation period among other predictors. Normalised milk yield per goat at 210 days was computed as NYj = d₁P₁+A + B, with NYj being normalised yield for goat j with P₁ being milk production in the first control, $A=30\mathrm{ }{\sum }_{i=1}^{n_j-2}\frac{1}{2}\left(P_i+P_{i+1}\right),$ and $B=\left[LC-30(n_j-1)\right]\frac{1}{2}\left(P_{n_j-1}+P_{n_j}\right),$ with $P_i$ being milk production in the i^th control and $n_j$ total number of controls of goat j. Then normalised productions at 210 days were computed as $MP210={\sum }_{i=1}^{n-1}\left[\left(\frac{{\mathit{pldc}}_i+{\mathit{pldc}}_{i+1}}{2}\right)·\mathrm{ }{I}_{i,\mathrm{ }i+1}\right]$ with MP210 being the accumulated milk yield until 210 days of lactation; pldc_i is milk yield under milk control i; pldc_i+1 is milk yield in the next control and I_i,i+1 is the days interval between two consecutive controls.

Fat, protein, dry matter, lactose and SCC × 10³ determination

Monthly milk sampling was performed at the Official Milk Quality Laboratory, Cordoba, Spain for the quantification of goat milk total protein, fat, dry matter and lactose using a MilkoScan analyser ™ FT1, while Fossomatic™ FC Somatic cell counting was used to test for SCC × 10³ sc/mL. After purging data, 2594 records from 159 goats were used for statistical analysis. Fat, protein, dry matter, lactose g were quantified after the percentages obtained for normalised milk yield at 210 days.

Preliminary statistical testing

The variables in our study (milk yield, fat, protein, dry matter, lactose contents and SCC × 10³ sc/mL) were numerical and continuous, while the factors considered in our model (farm, parturition season, month and year, birth number, control season, month, and year, number of controls, milking system, alive kid and dead kid number, birth type, days in milk, days to 1st control, days from last control to drying, drying month, season and year) were considered as categorical or ordinal, respectively (Table 1).

Table 1. Category description for non-genetic environmental factors considered in Murciano-Granadina goats.

Factor	Type	Levels
Farm	Nominal	25 farms
Parturition month	Ordinal	January, February, March, April, May, June, July, August, September, October, November, and December
Parturition season	Ordinal	Spring, Summer, Autumn and Winter
Parturition year	Ordinal	2006–2018
Birth type	Nominal	Abortion in lactation, Simple, Double, Triple, Quadruple
Parity order	Numeric	1–9 births
Alive kid number	Numeric	0–4 kids
Stillborn/Dead kid number	Numeric	0–2 kids (whichever the cause of their death is)
Number of control	Numeric	1–33 controls
OMET (Type of control)	Nominal	A4 (milking happens each 4 weeks), AT4M (milking happens each 4 weeks always in the morning), AT4T (milking happens each 4 weeks always in the afternoon), A6 (milking happens each 6 weeks), AT6M (milking happens each 6 weeks always in the morning) and AT6T (milking happens each 6 weeks always in the afternoon)
Days in milk	Numeric	0–182 Days
Days to first control	Numeric	0–225 Days
Days from last control to drying	Numeric	0–60 Days
Control month	Ordinal	January, February, March, April, May, June, July, August, September, October, November, and December
Control season	Ordinal	Spring, Summer, Autumn and Winter
Control year	Ordinal	2006–2018
Drying month	Ordinal	January, February, March, April, May, June, July, August, September, October, November, and December
Drying season	Ordinal	Spring, Summer, Autumn and Winter
Drying year	Ordinal	2007–2018

Milk yield and components are commonly preconceived to follow a normal distribution according to published research. However, properties of the data obtained from the field often challenge researchers, as the methods used or the conditions under which such records are collected, promote data to violate parametric assumptions. Hence new alternatives for the assessment of such data offer new opportunities to adequate the tools used to the nature of the data to be processed (Young et al. 2018). Data was filtered and animals whose records fell outside the ranges reported for the breed in bibliography were discarded. Parametric assumptions of normality, homoscedasticity, sphericity and multicollinearity were tested on complete historical records in the official dairy control of the Murciano-Granadina breed until 2018 (n = 2359479 records of 151997 goats) for milk yield, content (fat, protein, dry matter and lactose) and SCC × 10³ sc/mL variables. The Shapiro–Francia’s normality test routine of the ‘test and distribution graphics’ package of the Stata Statistical Software version 14.2 (StataCorp (2016), College Station, TX) was used to test the normality (Supplementary Table S1). The rest of parametric assumptions (Levene’s and Mauchly’s W tests and Tolerance and Variance Inflation Factor, respectively) were performed using SPSS Statistics for Windows, Version 25.0, IBM Corp (2017).

Afterwards, assumptions of normality were tested again clustering goats at the same lactational stage as a normal distribution of data could be expected as suggested by literature (Pizarro et al. 2019a).

Statistical analyses

Following these preliminary analyses, parametric assumptions were tested on our field data. As our field data violated parametric assumptions, a nonparametric approach was suggested. Hence, Kruskal–Wallis H test (nonparametric alternative to the One-Way ANOVA or One-Way ANOVA on ranks) was performed to study the potentially existing differences between levels of the same factor as suggested by Kay et al. (2019). A summary of the levels for all the categorical factors included in the model is shown in Table 1.

Then after conducting a Kruskal–Wallis H with three or more groups (k), we computed the strength effect of the factors on the variables tested. Kruskal–Wallis H produces Chi-square values with k-1 degrees of freedom and Chi-square values were transformed into F values using the expression modified from Murphy et al. (2014). Then, we calculated partial eta squared as Lakens (2013).

Partial eta squared (ηp²) was calculated considering eta squared measures the proportion of the total variance in a dependent variable that is associated with the membership of different groups defined by an independent variable. In this context, partial eta squared is a similar measure to eta squared, but which differs in the effects of other independent variables and interactions being partialed out. Additionally, although SPSS cannot calculate Cohen's f directly, it may be obtained from partial eta-squared. Cohen’s f is an extension of Cohen's d, which is the appropriate measure of effect size to use for a t test. Effect size is a quantitative measure of the ability to explain the variability that can be found in a variable cause by the different possibilities that a certain factor can take (Salkind 2010).

Once effects of significant factors had been identified, Dunn test was performed to analyse pairwise comparisons to detect differences between levels or categories within the same factor. In statistical hypothesis testing, incorrectly rejecting the null hypothesis based on the observed data is called making a Type I error (Salkind 2010). In multiple comparison, there is a presumable increase in Type I error deriving from the inclusion of multiple variables in models. The Bonferroni correction was applied to compensate for such an increase in Type I error. As all the variables had been previously reported to be non-normally distributed Shapiro–Francia’s tests (p < .001), within the same factor. Shapiro–Francia’s tests were carried out with the .sfrancia routine of Stata Statistical Software version 14.2 (2016) (StataCorp (2016), College Station, TX). All non-parametric tests were carried out using the independent samples package from the non-parametrical task of SPSS Statistics for Windows, Version 25.0, IBM Corp (2017) (IBM SPSS Statistics, Armonk, NY).

Then, CATREG with the Optimal Scaling procedure from the Regression task from SPSS Statistics for Windows, Version 25.0, IBM Corp (2017) was applied to issue specific regression equations to predict how milk yield (in g), fat, protein, lactose, and dry matter content (in %), and SCC × 10³ (sc/mL) linearly depended on the predictors which nonparametric tests had determined to be significant (p < .05). That is, CATREG standardised coefficients were evaluated to determine the number of standard deviation units that a certain independent variable linearly increased or decreased depending on the values of certain combinations of statistically significant predictors. As standardised coefficients were considered, factors measured using different scales or units can be used. As reported in IBM Knowledge Center (2019), the idea behind optimal scaling is to assign numerical quantifications to the categories of each variable, thus allowing standard procedures to be used to obtain a solution on the quantified variables.

Results

Supplementary Table 2 reports descriptive statistics with the aim to show the properties of the data used. Tables 2 and 3 show a summary of the results for Kruskal Wallis H Ranks test, partial eta squared and Cohen's f for all the levels of the factors of farm, parity order, birth number, number of alive and dead kids, month, season and year of birth, control number, type of control, month, season and year of control, days in milk, days from last control to drying, days to first control, month, season and year of drying affecting milk yield (g), fat (%), protein (%), dry matter (%), lactose (%) and somatic cell count (sc/mL) and their biological translation. Extended information in regards to the numeric results for all the test performed is provided in Supplementary Tables S3 and S4. Supplementary Table S3 shows complete numeric results for Kruskal Wallis H Ranks test, which was perform to detect differences in the distribution of milk yield and components across levels of the factors tested. The values of partial eta squared and Cohen’s f are also provided in Supplementary Table S3. Supplementary table S4, reports a summary of the numeric results for medians, Dunn’'s test and Bonferroni's significance correction of the aforementioned 19 effects (to determine which were the levels that accounted for the differences detected by Kruskal Wallis H, medians are provided to show the magnitude of differences and Bonferroni correction is applied to correct for the effect of variance inflation derived from including a lot of factors. All the factors integrated into the model significantly influenced all the variables studied (milk production (g), protein, fat, protein, dry matter, lactose (%) and somatic cells (sc/mL), with the exception of the number of kids born dead, which did not significantly affect fat, protein, dry matter and lactose content, and the year of drying, which did not significantly affect dry matter content. Then Supplementary Table S5 reports the results of CATREG. Supplementary Table S5 shows a model summary for categorical regression models and regression equations for the milk production and content quality related variables tested using standardised regression coefficients (β) for all factors. These models were designed at a previous stage of the study in which genetic effects through the analysis of SNPs were included (Pizarro et al. 2019b). Standardised coefficients are evaluated to determine how each unit of a certain independent variable linearly increases or decreases depending on the values of certain combinations of predictors. Then, R² and Adjusted R² of the model tell how these significant combinations of factors relate to describe and predict what happens and will happen in the future of dependent variables. R² and Adjusted R² are provided below the model for each variable in Supplementary Table S5.

Table 2. Summary of the results for Kruskal Wallis H Ranks test, partial eta squared and Cohen’s f for the levels of the factors of farm, parity order, birth number, number of alive and dead kids, month, season and year of birth affecting milk yield (g), fat (%), protein (%), dry matter (%), lactose (%) and somatic cell count (SCC × 10³ sc/mL) and their biological translation.

Variables Factors	Milk yield	Protein	Fat	Dry matter	Lactose	Somatic cells count
Farm	Significant differences across the 28 farms (600–5,347 g for milk yield, 3.005–4.660 % for protein, 4.290–7.210 % for fat, 12.530–17.740 % for dry matter, 3.790–4.980 % for lactose and 112.000–1970.000 sc/mL for somatic cells
Parity order	First parturition goats produced the lowest amount (1,832 g). Production then increased from the second birth to the fourth (2,696 g), to fall beyond this point (up to 1,400 g, ninth lactation)	Small increase as parity order increased (SD 0.092%)	Increased beyond fourth parturition (5.500 %)	No significant differences	Primipara goats produced the highest content (4.810%), followed by 2nd parturition does (4.780%). From 5th parturition on, lactose declined (until 4.590%)	Closely linked. Increases with parity order (from 378.000 to 1969.000 sc/mL from 1st to 9th parturition)
Birth number and number of alive kids	Higher as litter size (more than two alive kids) increased (from 2,052 to 3,240 g). In goats that aborted, but begin lactation, production was significantly higher (2,418 g) than in single-birth goats (2,052 g) and similar to double-birth goats (2,434 g)	Goats that aborted and continued in lactation, produced slightly higher percentage of protein (3.670%). No difference for goats giving birth to 1– 3 alive kids (SD 0.066 %)	Decreased linearly as the litter size increased (more than one alive kid) from 5.440 to 4.880 % for fat, from 14.600 to 13.750 for dry matter and from 4.760 to 4.610 % for lactose			Lower somatic cells counts were present in simple birth goats (452.500 sc/ml) increasing with alive kid number up to 1012.000 sc/mL (4 alive kids). Goats that aborted but lactated, somatic cells increased compared to single births (782.000 sc/mL)
Born dead kids	When some of the kids were born dead, milk production increased (2,490 g). However, production decreased significantly in births where both twins were born dead (900 g)	No effect was found				The birth of one dead kid implied an increase (1285.000 sc/mL), unlike what happened when two kids were born dead (384.000 sc/mL)
Month of birth	Goats kidding in January and March produced more milk (2,815–2,820 g). However, yield decreased in goats born during June to October (from 2.717 to 1,699 g)	The opposite trends described for milk yield			Goats kidding in the month of October and November produced more lactose (4.840–4.870 %)	Births in February related to lower somatic cells counts (295.000 sc/ml), with a tendency to rise from March, with a peak in July (1053 sc/mL), which decreased towards December up to 493.000 sc/ml
Season of birth	Goats kidding in autumn and winter reported the highest milk yields (2,575–2,686 g), decreasing during spring and summer (1,960 to 2,226 g)	Decreased from spring-summer to autumn-winter from 3.450 to 3.620 % for protein, from 5.530 to 5.000 % for fat, from 13.925 to 14.790 % for dry matter and from 4.670 to 4.810 % for lactose				Lowest somatic cells counts recorded in winter (444.000 sc/mL), increasing towards spring (747.500 sc/mL) to decrease in summer (595.000 sc/mL)
Year of birth	An increase was observed over the years up to reach 2,981 g, decreasing progressively towards 2018 (1,344 g)	Constant over the years with a peak in 2018 of 5.450 % in protein, 8.140 % in fat and 5.450 % in 18.920% in dry matter			Progressively decreased until 2018 (from 4.870 % in 2009 to 4.410 % in 2018)	Fluctuated over the years (from 305.000 to 4632.000 sc/mL), significantly increasing from 2017 on

The results presented are significant at p < .05 at minimum [for specific statistical information see Supplementary Tables S2, S3 and S4.]

Table 3. Summary of the results for Kruskal Wallis H Ranks test, partial eta squared and Cohen’s f for the effects of all the levels of the factors of control number, type of control, month, season and year of control, days in milk, days from last control to drying, days to first control, month, season and year of drying affecting milk yield (g), fat (%), protein (%), dry matter (%), lactose (%) and somatic cell count (SCC × 10³ sc/mL) and their biological translation.

Variables Factors	Milk yield	Protein	Fat	Dry matter	Lactose	Somatic cells count
Number of Control	Relatively constants (Variation coefficient of 34.552 % for milk yield, 10.597 %, for protein 7.206 % for fat, 5.408 % for dry matter and 2.417 % for lactose. Milk yield and lactose described the opposite trend to protein, dry matter and fat with median g and percentual contents progressively and slightly decreasing from first to last control					Progressively increased from third control on (from 377.500 sc/mL to 2.868 sc/mL)
Type of Milk Control	A4, AT4T and AT4M controls collected a higher amount of milk (μ: 2,346 g) compared to A6, AT6M and AT6T controls (μ: 1,440 g)	Highest quantities recorded for A6 control (3.950 %)	The highest fat (μ: 5.998 %) and dry matter (μ: 15.040 %) contents collected in AT4T, A6 and AT6T controls		A6 control reported the lowest contents (4.255 %), while AT4M and AT6M controls reported the highest values (μ: 4.820 %)	The AT4M control reported the lowest values (508.000 sc/mL). No differences reported for the rest of control types except for AT6T (1021.000 sc/mL)
Control month and season	Remained stable throughout the year (Variation coefficient of 7.455%), with the exception of November to January, when it slightly decreased (from 2,700 to 2,112 g). Higher milk yields collected from Spring controls (2,500 g). No differences between autumn and summer. Decreased in winter up to 2,200 g	Decreased during the months from April to September, with July and August being the months with the lowest contents (3.220 % for protein, 4.475 % for fat, 13.230 % for dry matter and 4.600 % for dry matter). Highest values for these variables were reached during winter months (3.700 % for protein, 5.760 % for fat, 15.100 % for dry matter and 4.840 % for lactose while the lowest ones were reported during spring (3.330 % for protein, 4.620 % for fat, 13.430 % for dry matter and 4.610 % for lactose)				Remained stable from January to May (Variation Coefficient of 4.686%), slightly increasing from June-December (Variation coefficient of 13.511%), except for August (509.500 sc/mL). Overall, slightly higher in spring-summer (μ: 685.500 sc/mL) than in autumn-winter (μ: 541.000 sc/mL)
Control year	Lower amounts were reported prior to 2013 (μ: 2,038 g). Peak reached in 2014 (3,020 g), then levels decreased in 2017 (1,888 g) slightly recovering in 2018 (2,352 g)	Relatively stable over the years (SD 0.147 %), increasing from 3.440 % in 2008 to 3.760 % in 2018	Remained stable over the years (Variation coefficient of 5.064% for fat and 1.541 % for dry matter)		Decreased over the years from 4.855% in 2006 to 4.535 % in 2018	Increased over the years from 317.000 sc/mL in 2008 to 1480.000 sc/mL in 2018
Days in Milk (DIM)	Slight decrease in production from day 30 (before 30 DIM, μ: 2,380 g, after 30 DIM, μ: 2,049 g)	Maintains stable values throughout lactation (Variation coefficient 7.841%)	Significantly increased (before 30 DIM, μ: 5.164 and 14.282 %, after 30 DIM, μ: 5.465 and 14.676 % for fat and dry matter, respectively)		Maintains stable values throughout lactation (Variation coefficient 3.376 %)	Describes oscillating values throughout the lactation period (Coefficient of variation of 79.658 %)
Days from last control to drying	Remained relatively stable with slight significant increases as days from last control progress to drying (Coefficient of variation of 27.022% for milk yield, 7.856 % for protein, 12.645% for fat, 5.857 % for dry matter and 3.407 % for lactose contents. Lactose content decreases as period between last control to drying shortens from 4.965 % (2 months prior to drying) to 4.330 % (3 days prior to drying)					Decreased as number of days from last control to drying reduce from 112.500 sc/mL (55 days prior to drying) and 2579.000 sc/mL (7 days prior to drying)
Days to First Control	Levels become irregular when there is a month or more to first control (Variation coefficient: 30.709%)	Stable levels with slight increases as days to first control increase with variation coefficients of 7.637 % for protein, 11.069% for fat, 5.711 % for dry matter and 3.722 % for lactose				Highest counts (2320.000 sc/mL) when it took over two months to first control
Month and season of drying	Goats drying in October produced the highest yields (2,989 g), followed by goats drying in February, March, and September to December (μ: 2,723 g). Goats drying in April produced the lowest milk yield (1,287 g). Goats drying in autumn produced the highest milk yields (2,870 g)	April (spring) was the month reporting the highest contents (3.845 % for protein, 5.930 % for fat and 15.260 % for dry matter). Contents slightly decreased from September to December (autumn to winter) around 0.150 %			Stable levels across months, only slightly decreasing from August to November, to slightly increase again until July-August after which the same cycle repeats (SD 0.059 %)	Increased from January to April (with a peak in spring of 1182.000 sc/mL), decreasing afterwards in autumn, when the lowest values are reached (488.000 sc/mL)
Drying year	Progressive increase until 2014 (3,100 g), when levels decreased until 2,055 g in 2017	Progressive increase over the years from 3.150 to 3.860 %	No significant effect (p > .05)		Decreased over the years from 5.060 to 4.545 %	Progressive increase over the years from 299.000 sc/mL to 1221.000 sc/mL

The results presented are significant at p < .05 at minimum [for specific statistical information see Supplementary Tables S2, S3 and S4.]

Discussion

Selection efficiency for milk quality related traits can only be achieved through the identification and inclusion in models of extrinsic and intrinsic factors which may condition the non-genetic performance of individuals. Isolating these effects, we can determine their exact impact and provide information that can be very useful to achieve desirable levels in profitable traits, regulate such effects, or counteract situations that may be economically detrimental. The main aims of the models used in milk industry should focus on how to control the factors influencing milk yield and composition and which the repercussions from these productive requirements will be at long-term improving the accuracy of the predicted breeding values (Pizarro et al. 2020b).

All the traits considered were within the standards found for other breeds of goats (Brito et al. 2011). Our results may suggest a high variability exists between farms. Such variable patterns could be attributed to the conglomerate of factors that the farm involves as a factor itself. In these regards, some authors have suggested climatic variations at a given place, farm management policies (traditional, intermediate or modern), altitude (plains, hills, mountains), farm size (small, medium, large), the nutritional quality of the feed provided to the animals or the composition of the herd (Jaafar et al. 2018; Vacca et al. 2018) may presumably be some of the factors on which the differences reported for our results may base. For instance, in the case of somatic cells, the size of the farm and the composition of the herd may have been decisive as, under confined animal conditions, typical of intensive systems, animals are both more exposed to intramammary infections and to a greater likelihood of large scale transmission by direct contact between animals (Lôbo et al. 2017).

Birth number has also a significant impact on milk production, as milk yield increases with birth number, as a result of an increased hormonal activity, metabolic activity, secretory cell activity, udder volume and intake of nutrients used in milk synthesis (Knight and Peaker 1982; Carnicella et al. 2008). The effect of birth number on milk yield might be attributed to the proportion of mammary alveoli from the previous lactation. Those alveoli would then be added to alveoli that developed in subsequent lactations. This continuity was interrupted as the age of the goats increased (Merkhan and Alkass 2011). Fat content is significantly and directly related to parity order as a subsequent increase in fat content progressively occurs as goat’s parturition number increases. This finding was supported by Strzałkowska et al. (2010), who reported goat milk characterises by a low level of fatty acids (<1.0 mEq/L) in relation to the SCC × 10³, animals’ age and stage of lactation. These authors suggest that highest contents of fatty acids and fat were recorded in the last stage of lactation, while the highest milk yield was observed in primiparous does. The milk obtained from this group of goats was, however, the most susceptible to lipolysis process. Contrastingly, Vacca et al. (2018) did not observe a significant effect of the parity order on the content of fat and protein, which, agrees with the findings found for protein in our study. According to Addass et al. (2013) parity order influences the amount of lactose produced, with goat in their third parturition producing more quantity. However, in our study such increased lactose content was reported by goats in their first two births, linearly descending as parity order increases, which was also described by Vacca et al. (2018). Oppositely, the higher the values for parity order, the higher the SCC × 10³ sc/mL are as well (Vacca et al. 2018). SCC × 10³ sc/mL is an indicator of mammary gland health state and lactose can be reduced with mastitis while SCC × 10³ sc/mL increases (Hanuš et al. 2010). According to the research reported by Paape et al. (2007), from the fifth birth onwards, the SCC × 10³ sc/mL increases significantly, which may be a consequence of anatomical-physiological changes in the mammary gland promoting transfer of leukocytes from the blood to the milk (Mehdid et al. 2019).

Regarding the type of birth, an increase in milk production was observed as the number of kids born increased as supported by other authors (Carnicella et al. 2008; Ibnelbachyr et al. 2015). The increase in the amount of milk produced as litter size increases may be ascribed to hormonal stimulation derived from multiple foetuses, which could be explained by levels of placental lactogen, progesterone and prolactin during gestation, given their stimulating action on the mammary glands (Lôbo et al. 2017). However, Carnicella et al. (2008) and Brito et al. (2011) did not report an influence of the type of birth on the content of fat and lactose, contrary to what is described in the present work. Perhaps the fact that these authors only recorded animals with single and twin-parities may be the basis for such findings. The counts of somatic cells are higher as goat accumulates more births, a fact which was corroborated by Pleguezuelos et al. (2015) for the same breed.

According to Galina et al. (1995), multiple births can extend the interval between subsequent births due to the stress suffered by the goat, since to maintain the gestation of one or more kids, a a greater mobilisation of nutrients during gestation and lactation may be required. Hence, a longer period of recovery of the goats until the next conception may be necessary as well. A failure to complete this recovery process (an incorrect or incomplete recovery) can damage and result detrimental for the next lactation, which may influence milk quality and composition.

However, in the context of our research, the kids in this study were artificially nursed from birth. In cases when kids are not left to nurse with the dam, according to Filho et al. (2004), even if in multiple births occur, goats receive a constant stimulus during milking, so the effect of the type of birth could be rather ascribed to differences during gestation than to subsequent lactations. This may also explain the differences observed in the present work for the production of milk regarding the number of kids born alive and dead.

The effect of the month of birth has a significant effect on the production of milk and its components, with goats born in the months of January, March, April, May and June reaching the highest production. However, Crepaldi et al. (1999) observed significant differences in favour of the months of January and February with respect to March and April. The increase in production in the winter months may be due to lower temperatures, greater availability of food, lower incidence of disease and lower prevalence of external and internal parasites (Assan 2015).

The season in which the births take place also has a significant impact on the production of milk and its components in the Murciano-Granadina goat, with the goats that give birth in autumn and winter being the ones that produce the greatest quantity of milk. However, Gómez et al. (2002) found that the highest yields were observed in spring births in the Murciano-Granadina breed. According to Carrizosa et al. (1993) the productive differences between lactations are mainly due to differences in the duration of lactations, since autumn births correspond to prolonged lactations and therefore higher production. However, the reduction in production in the spring and summer seasons may be due to the high temperatures recorded in southern Andalusia during the second half of spring and summer, which may generate thermal stress in the animals with the consequent reduction in feed intake and the increase in additional activities such as maintenance and homeothermy, while promoting chemical reactions in the body and the production of thermal shock proteins, which consume large amounts of adenosine triphosphate (ATP) (Salama et al. 2014).

Goats display the maximum appetite rates between 0° C and 10° C, gradually decreasing as the temperature increases from 40 °C (Appleman and Delouche 1958). These results do not agree with those described by Ibnelbachyr et al. (2015), who observed higher milk production during spring and summer. However, as described by Ibnelbachyr et al. (2015) fat, dry matter and lactose content are higher in goats born in spring and summer, contrary to which was described by Zahraddeen et al. (2007) who suggested that during dry seasons, fat and lactose content decreases as a result of high temperatures and reduced feed intake. Some authors attribute this difference in seasons to climatic conditions or the availability of pasture (Ishag et al. 2012).

As can be seen in the results, the effect of the year of parturition on milk production is significant, a fact that was previously described by Pleguezuelos et al. (2015) in the Murciano-Granadina goat. Such variations may be the consequence of climatic changes, fluctuations in nutrient availability and herd composition (age and status) (Ishag et al. 2012; Ibnelbachyr et al. 2015). Some authors also attribute these annual variations in milk production and its components to the characteristics of pastures, which vary based on environmental conditions (Piergiovanni and Casassa 1982). The reasons described above would also explain the variations observed for milk production and its components derived from the control year factor.

According to Lôbo et al. (2017) the number of controls performed on goats significantly affects the production of milk, fat, lactose and somatic cells with the exception of protein, describing how the lactose content decreases as the number of controls increases, while fluctuations are observed for the rest of traits.

The type of milk control used has not been considered by any other previous research. Our results are indicative of the fact that the application of each of them may significantly and differently condition the variables studied. Therefore, selecting the appropriate milk control to be used in dairy farms is of utmost importance as it may enhance the information considered by the models implemented within the breeding programme of the breed.

Herds of the selection nucleus are selected not only for their production and quality data but for the genetic connections among them. Our results suggest the lower fat, protein, dry matter and lactose contents are recorded during the dairy controls carried out during the months of July and August, with the months between November and February being the months with the highest milk production. These results are in line with those described by Kljajevic et al. (2018) who observed that the lowest values of fat, protein, dry matter and lactose content were reported during the month of August, while the highest ones were reported in the month of December in Saneen goats. Other authors suggest that the lowest values for fat and protein in milk are provided during the months of June and July (Goetsch et al. 2011) or June to August (Mayer and Fiechter 2012), while the highest ones are obtained in November and December or October, respectively.

Mayer and Fiechter (2012) reported seasonal variations are a consequence of the physiological changes that the goat experiences during lactation. In the case of the Murciano-Granadina goat, a significantly higher quantity of milk was registered during spring. However, Lôbo et al. (2017) reported the highest yields during winter months in studies conducted in Brazil (southern hemisphere). Similarly for fat, protein, dry matter and lactose content, winter is the season in which the highest values are significantly recorded, as described by Lôbo et al. (2017) (in summer, southern hemisphere) and contrary to de Souza et al. (2009) (in winter, southern hemisphere). For the protein content, previous studies carried out in Murciano-Granadina goats corroborate a greater quantity during winter months (Dı́az et al. 1999).

The average values observed for somatic cell content in the present work are similar to those commonly described in the dairy goat sector (270·10³ and 2,000·10³ cells/mL) (Granado et al. 2014). The controls carried out in spring and summer reported significantly higher levels, as supported by de Souza et al. (2009) in studies conducted in Brazil (Southern Hemisphere). However, Lôbo et al. (2017) reported SCC × 10³ sc/mL may fluctuate depending on the season.

Milk production has also been reported to be influenced by the days in milk. Our results suggest a decrease in milk production occurs from the day 30 of lactation onwards. According to León et al. (2012) in the Murciano-Granadina goat, maximum yield occurs on average 45 days after lactation, which is similar for other populations of dairy goats, where the peak usually takes place between the fifth and eighth week of lactation (Morand-Fehr and Sauvant 1980; Kala and Prakash 1990). Compared to other Spanish breeds these peaks occur in the second and third week in Majorera (Fresno et al. 1994), in the fourth week in Malagueña (Herrera et al. 1984) and Tinerfeña goats (Fresno et al. 1994).

On the one hand, fat and dry matter content significantly increase with days in milk, as it has been reported for other breeds, such as Polish White Improved goats and Draa (Olechnowicz and Sobek 2008; Ibnelbachyr et al. 2015). However, Mestawet et al. (2012) reported fat content to be higher at the beginning and end of lactation, while other authors observed a decrease in production as lactation progresses (Mahmoud et al. 2014).

On the other hand, protein maintains stable values through lactation, as described by El-Tarabany et al. (2018) suggesting such stability can be attributed to the ability of native breeds to preserve body condition at different stages of lactation. However, Mestawet et al. (2012) and Ibnelbachyr et al. (2015) observed a higher protein fraction at the beginning and end of lactation, which contrasts the findings by Strzałkowska et al. (2009) who reported an increase as lactation progressed. Lactose also describes a similar behaviour to that of milk production, which was also corroborated by other authors (Olechnowicz and Sobek 2008; Mahmoud et al. 2014; El-Tarabany et al. 2018).

Somatic cell content presents oscillating values throughout the lactation period. However, de Souza et al. (2009) described a lower incidence in early lactations which increased through lactation. Similarly, El-Tarabany et al. (2018) described how total leukocytes decreased in the initial stage of lactation compared to mid and late stages. This reduction in early stage can be attributed to the increase in cortisol (as a result of the increase in milk yield towards maximum level), which is probably responsible for the deterioration of the immune response (Caroprese et al. 2012).

The number of days in milk to the first milk control have also reported a significant impact on the production of milk and milk components. Such effect has also been reported for the days from the last control to the drying of the goats. However, increasing stage of lactation, as measured by increased days in milk and months of the year with the highest mean days in milk, was associated with increased SCC × 10³ sc/mL in goats. The drying period before parturition is necessary to regenerate damaged and senescent epithelial cells to maximise the production of the next lactation, since this depends on the number of epithelial cells and their activity (Knight 2000). The effect of the month of drying, the drying season and the year of drying have not been further studied in literature, despite there is a clear effect on the production of milk and its components as suggested by our results, which should be considered to perform appropriate farm management policies. For instance, systematic treatment of goats at drying-off is an efficient method for the cure of subclinical mastitis and control of SCC × 10³ sc/mL at the beginning of the following lactation (Poutrel et al. 1997).

Conclusions

The thorough analysis of the influence of non-genetic factors on the production of milk and its components can be crucial for the maximisation of the efficiency of the evaluation models used in breeding programme in the Murciano-Granadina goat breed. In addition, as it is observed in the present work, the inclusion as fixed effects of factors such as the type of dairy control and parameters related to the drying-off event, which are not likely to be studied nor considered, provide interesting information, which should be included in genetic evaluations, given their role as promoters of a greater milk yield and a better compositional quality. The improvement of the profitability of highly selected breed, which becomes more challenging as selection degree increases, can be maximised through the evaluation of certain combinations of factors which may detect or identify critical points in the productive cycle to control, which in turn can help us making the most of the productive potential of breeds.

Ethical approval

The study followed the premises described in the Declaration of Helsinki. The Spanish Ministry of Economy and Competitivity through the Royal Decree-Law 53/2013 and its credited entity the Ethics Committee of Animal Experimentation from the University of Córdoba permitted the application of the protocols present in this study as cited in the fifth section of its second article, as the animals assessed were used for credited zootechnical use. This national Decree follows the European Union Directive 2010/63/UE, from the 22nd of September of 2010. Furthermore, the present study works with records rather than live animals directly, hence no special permission was compulsory.

Acknowledgements

The authors thank the National Association of Breeders of Murciano-Granadina Goat Breed, Fuente Vaqueros (Spain) and the PAIDI AGR 218 research group for their support and assistance.

Disclosure statement

No potential conflict of interest was reported by the author(s).