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.


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 3 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 3 sc/mL (quantified in somatic cells per millilitre, sc/mL).

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 3 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(Pizarro et al. , 2020b and epistatic effects (Pizarro et al. 2020a(Pizarro et al. , 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 3 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, b ¼ P À1 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: where PBV is the predicted breeding value for each of the traits and animals included in the matrix; W 1 is the weight for milk yield, W 2 for fat, W 3 for protein and W 4 for lactose in g and normalised at 210 days; and l 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 3 (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.
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 1 ¼ FC-KD and the days between the penultimate (PC) and the last control (LC) as follows, d 2 ¼ 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 1 P 1 þA þ B, with NYj being normalised yield for goat j with P 1 being milk production in the first control, A ¼ 30 , 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 h i 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 3 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 TM FT1, while Fossomatic TM FC Somatic cell counting was used to test for SCC Â 10 3 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 3 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).
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 3 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  (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 (gp 2 ) 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 3 (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 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 singlebirth 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%  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 (b) 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 2 and Adjusted R 2 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 2 and Adjusted R 2 are provided below the model for each variable in Supplementary  Table S5.

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 3 , 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 3 sc/mL are as well (Vacca et al. 2018). SCC Â 10 3 sc/mL is an indicator of mammary gland health state and lactose can be reduced with mastitis while SCC Â 10 3 sc/mL increases (Hanu s et al. 2010). According to the research reported by Paape et al. (2007), from the fifth birth onwards, the SCC Â 10 3 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 omez 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)  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 3 and 2,000Á10 3 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 3 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 on 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 3 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 3 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.