Cheese produced at the individual cow level is the result of a complex process with plenty of elements involved, from milk quality characteristics (e.g. percentage of protein and fat in the milk), milk coagulation properties (MCP), curd firmness (CF) and cheese-yields, such as the quantity of cheese obtained from a given amount of milk (CY) and the loss of nutrients in whey/recovery in curd (REC). Although detailed phenotyping of individual animals is required to understand the biological and genetic background of the traits of interest, inclusion of a large number of traits into a selection index in breeding programs poses difficulties in interpretation and computations.
To overcome the problem of data dimension, and for a better understanding of complex phenomena, multivariate factor analysis (MFA) has been investigated for a variety of traits, such as milk composition, MCP and CY traits in dairy cattle, sheep and goats (Macciotta et al. Reference Macciotta, Cecchinato, Mele and Bittante2012; Manca et al. Reference Manca, Serdino, Gaspa, Urgeghe, Ibba, Contu, Fresi and Macciotta2016; Vacca et al. Reference Vacca, Paschino, Dettori, Bergamaschi, Cipolat-Gotet, Bittante and Pazzola2016), for milk fatty acids in dairy cows (Conte et al. Reference Conte, Serra, Cremonesi, Chessa, Castiglioni, Cappucci, Bulleri and Mele2016; Mele et al. Reference Mele, Macciotta, Cecchinato, Conte, Schiavon and Bittante2016) and for cheese volatile organic compounds (Bergamaschi et al. Reference Bergamaschi, Biasioli, Cappellin, Cecchinato, Cipolat-Gotet, Cornu, Gasperi, Martin and Bittante2015).
Factor analysis belongs to the general framework of multivariate analysis. The main idea of MFA is that n observed variables, x, can be expressed as linear functions of p (p < n) latent variables (factors; Fs): this statistical approach focuses on understanding relationships (the underlying latent concept that the measured variables represent) among a set of observed variables. Thus, Fs can be considered as variables that are not measurable, but they can be extracted from a set of measured-indicator variables by making use of their covariance structure (Bollen, Reference Bollen1989). In such a way, only the underlying concept of interest is kept for further analysis while at the same time data reduction is achieved.
To our knowledge, there are no studies reported in the literature that deal with MFA and milk protein profile. Moreover, milk protein fractions have never been considered in combination with milk technological traits to represent the complexity of cheese-making process. In addition, besides the statistical methodology, no studies have been reported in the literature dealing with the effects of feeding regime and management system on milk proteins and cheese-making related variables. Therefore our objectives were (i) to create a new set of latent phenotypes related to milk quality, detailed protein composition measured by the RP-HPLC method, coagulation properties, and cheese-making traits, by using MFA, and (ii) to assess some individual sources of variation (i.e., stage of lactation and parity) and herd factors as well as the effects of feeding and management systems.
Materials and methods
Animals, herds and dairy farming systems
Milk samples from 1264 Italian Brown Swiss cows reared in 85 herds located in Trento Province (Italy) were collected. A full description of the sampling procedure can be found in (Cecchinato et al. Reference Cecchinato, Cipolat-Gotet, Casellas, Penasa, Rossoni and Bittante2013). In brief, ~15 cows/herd were individually sampled once (evening milking). Samples were processed within 20 h after collection. Information on cows and herds was supplied by the Breeders Association of Trento Province.
The farming systems have been previously analysed and reported in Sturaro et al. (Reference Sturaro, Cocca, Gallo, Mrad and Ramanzin2009, Reference Sturaro, Marchiori, Cocca, Penasa, Ramanzin and Bittante2013). In brief, two main categories of farming systems were distinguished: (i) traditional farms and (ii) modern dairy systems. The traditional systems represent small and old barns where feeding is heavily based on meadow hay and cows are tied. Modern farms, with loose cows and milking parlour, were further distinguished in 3 sub-categories depending on the feeding system: (a) modern dairy system but without use of total mixed ratio (TMR), (b) modern dairy farms with silage-based TMR and (c) modern dairy farms without silage using water to moisturise TMR.
Cheese-making phenotypes
Analysis of milk and milk protein fractions
Individual milk subsamples were analysed for protein, fat, lactose and casein contents using MilkoScan FT6000 (Foss, Hillerød, Denmark). The pH of the subsamples was measured using a Crison Basic 25 electrode (Crison, Barcelona, Spain). Somatic cell count measures were obtained by Fossomatic FC counter (Foss, Hillerød, Denmark) and were logarithmic transformed [SCS = log2(SCC/100 000) + 3] (Ali & Shook, Reference Ali and Shook1980).
Protein fractions [αs1-, αs2-, β- and κ- casein (CN), β-lactoglobulin (β-LG) and α-lactalbumin (α-LA)] were measured by the RP-HPLC method (Bonfatti et al. Reference Bonfatti, Grigoletto, Cecchinato, Gallo and Carnier2008) and were expressed as proportions to the total milk nitrogen (N) content. Further, the phosphorylated form of the αs1-CN was obtained as proposed by (Bonfatti et al. Reference Bonfatti, Cecchinato, Gallo, Blasco and Carnier2011). The remained N milk compounds were also included in the analysis and were calculated by subtracting the sum of the protein fractions from the total N milk content.
Milk coagulation properties and modelling the CF
Measures of milk coagulation properties were obtained using the Formagraph instrument by Foss Electric A/S according to the procedure described in Cecchinato et al. (Reference Cecchinato, Cipolat-Gotet, Casellas, Penasa, Rossoni and Bittante2013). Files containing 360 CF values for each milk sample, recorded every 15 sec for 90 min, were retrieved and used to estimate a set of parameters of CF at time t (CFt) according to equations and methodology developed by Bittante et al. (Reference Bittante, Contiero and Cecchinato2013b). Estimated parameters included: rennet coagulation time (RCTeq, min), potential asymptotical curd firmness (CFP, mm), representing the maximum potential curd firmness of a given sample after infinite time in the absence of syneresis, curd-firming rate constant (kCF, % × min−1) which measures the relative increment of CF toward CFP, that is predominant until reaching the maximum curd firmness (CFmax, mm; at time t max, min) and before syneresis [measured by kSR (% × min−1)] became prevalent. To avoid convergence and estimation problems, CFP was calculated multiplying CFmax by 1·34, that is the coefficient resulting from the linear regression between CFP and CFmax values obtained in a preliminary analysis (Stocco et al. Reference Stocco, Cipolat-Gotet, Bobbo, Cecchinato and Bittante2016). The other three CFt model parameters (RCTeq, kCF, and kSR) were estimated by curvilinear regression using the nonlinear procedure (PROC NLIN) in the SAS software (SAS Institute Inc., Cary, NC).
Individual cheese yield and curd nutrient recoveries
A detailed description of the individual model-cheese processing can be found in Cipolat-Gotet et al. (Reference Cipolat-Gotet, Cecchinato, De Marchi and Bittante2013). The phenotypes were obtained through a model cheese-making procedure on 1500 ml of milk for each cow. In brief, the traits analysed were: (i) three %CY traits, expressing the weight (wt) of fresh curd (%CYCURD), of curd dry matter (%CYSOLIDS), and of water retained in the curd (%CYWATER) as percentage of wt of milk processed, and (ii) four REC traits representing the proportion of nutrients and energy of the milk retained in the curd (RECSOLIDS, RECFAT, RECPROTEIN and RECENERGY calculated as the % ratio between the nutrient in curd and the corresponding nutrient in processed milk). The recovery energy in the curd was calculated as the difference between energy in the milk and in the cheese (NRC, 2001).
Statistical analysis
Multivariate factor analysis
The approach has been described in Dadousis et al. (Reference Dadousis, Cipolat-Gotet, Bittante and Cecchinato2017a). In brief, to avoid severe multicollinearity problems, 3 out of the 29 phenotypes (CFmax, %CYCURD and RECSOLIDS) were excluded. As mentioned above, CFmax is proportional to CFP. Moreover, the %CYCURD is the sum of %CYWATER and %CYSOLIDS, while RECSOLIDS has phenotypic correlation with RECENERGY and %CYSOLIDS greater than 0·9. The remaining twenty-six phenotypes were simultaneously analysed in the following factor model:
$${\bi x} = {\bi \Lambda \xi} + {\bi \delta} {\rm,} $$
where x is a vector containing the measured phenotypes, ξ is a vector of the factors, Λ contains the factor loadings (λ) relating the factors to the original variables and δ is a vector of the residuals.
At a first step, the difference between Pearson and partial correlations of the measured variables was used as a hint to assess the adequacy of the data for MFA. This difference can be viewed as a way to control if the correlation between 2 variables is mediated by other variables, with a high value indicating the existence of a latent structure. The Kaiser-Meyer-Olkin (KMO) measure of sampling adequacy was applied to quantify this difference (Dziuban & Shirkey, Reference Dziuban and Shirkey1974; Kaiser & Rice, Reference Kaiser and Rice1974). Further, MFA was applied. The factors were extracted based on prior knowledge, biological interpretation and the proportion of original variance explained. Moreover, factor rotation (varimax) was used to identify simple structure. Factor loadings >|0·4| were considered as ‘significant’ to explain the factors. The MFA was performed with the psych package (Revelle, Reference Revelle2014) in R (R Core Team, 2013).
Mixed model analysis
To estimate sources of variation related to the factors, the 10 factor scores were analysed with the following model:
$$\eqalign{y_{ijklm} & = \mu + {\rm dairy}\,{\rm system}_i + {\rm herd}_j ({\rm dairy}\,{\rm system})_i \cr & + {\rm parity}_k + DIM_l + e_{ijklm}} $$
where y
ijklm
is the observed phenotype (i.e, the factor scores); μ is the overall mean; dairy system
i
is the fixed effect of the ith dairy system (i = 1 to 4); herd
j
(dairy system)
i
is the random effect of the jth herd (j = 1 to 85)
$\sim N(0,{\bf I}\sigma _h^2 )$
nested within the ith dairy system; parity
k
is the fixed effect of the kth parity (k = 1 to 4 or more lactations); DIM
l
is the lth 30-d class of days in milk (11 classes); and e
ijklm
is the residual random error term
$\sim N(0,{\bf I}\sigma _e^2 );\; $
where I is an identity matrix,
${\rm \sigma} _{\rm h}^2 $
, and
${\rm \sigma} _{\rm e}^2 {\rm \;} $
are herd/date and residual variances, respectively. The significance of the dairy system was tested on the error line of the herd within dairy system, while for parity and DIM class the error line of the residual variance was used. Orthogonal post hoc contrasts (P < 0·05) were built for dairy system: (i) the ‘Traditional’ dairy system was compared with the ‘Modern’ systems; (ii) within the modern systems, the ‘No TMR’ herds were compared with the TMR herds; (iii) within the TMR herds, those that use silage were compared with those that use water.
Results
Factors
Descriptive statistics of the full dataset including all 29 phenotypes are shown in Table 1. The phenotypic Pearson and partial correlations of the 26 traits analysed in MFA are presented in Fig. 1. The average KMO value in our dataset was 0·55. Ten FS, explaining 74% of the original variance, were kept for further analysis. The varimax rotated factor loadings (sorted by maximum variance explained) are shown in Tables 2 and 3.
Fig. 1. Pearson (above the diagonal) and partial (under the diagonal) phenotypic correlations and among the 25 traits used in the factor analysis. On the diagonal the Kaiser-Meyer-Olkin (KMO) measure of sampling adequacy per trait.
Table 1. Summary statistics of milk (yield and quality), protein fractions, coagulation (curd firming) and cheese-making (%CY and REC) traits
† SCS = log2 (SCC/100 000) + 3. Milk protein fractions: CN, casein; LA, lactalbumin and LG, lactoglobulin. Curd firming: RCTeq, estimated RCT; CFP, asymptotical potential value of CF; kCF, curd-firming instant rate constant; kSR, syneresis instant rate constant; CFmax, maximum curd firmness achieved within 90 min; and t max, time at achievement of CFmax. %CY, ratios of the weight (g) of the fresh curd (%CYCURD), curd dry matter (%CYSOLIDS) and curd water (%CYWATER) vs. the weight of the processed milk (g); REC, ratio of the weight (g) of the curd constituent (dry matter, fat, protein or energy, respectively) vs. that of the same constituent in the processed milk (g).
‡ Standard deviation.
§ Coefficient of variation.
Table 2. Rotated factor (F) pattern, factor name and variance explained by the Fs for F1 through F5 †
† Factors have been sorted based on proportion of variance explained. Values above |0·4| in bold. F1: %CY, Factor related to the percentage of individual cheese yield; F2: CFt, Factor related to the curd firmness; F3: Yield, Factor related to the milk yield; F4: Cheese N, Factor related to the milk nitrogen that is present into the cheese curd; F5: as1-β-CN, Factor related to the as1- and β-CN contents in milk, expressed as relative contents to the total milk nitrogen.
‡ SCS = log2 (SCC × 100 000) + 3. Milk protein fractions: CN, casein; LA, lactalbumin and LG, lactoglobulin. Curd firming: RCTeq, estimated RCT; CFP, asymptotical potential value of CF; kCF, curd-firming instant rate constant; kSR, syneresis instant rate constant; and t max, time at achievement of CFmax. %CY, ratios of the weight (g) of the curd dry matter (%CYSOLIDS) and curd water (%CYWATER) vs. the weight of the processed milk (g); REC, ratio of the weight (g) of the curd constituent (fat, protein or energy, respectively) vs. that of the same constituent in the processed milk (g).
Table 3. Rotated factor (F) pattern, factor name, communality (com) † of variables and variance explained by the Fs for F6 through F10 ‡
† Communality, the sum of the squared factor loadings per trait.
‡ Factors have been sorted based on proportion of variance explained. Values above |0·4| in bold. F6: Udder health, Factor related to the udder health of a cow; F7: κ-β-CN, Factor related to the κ- and β-CN contents in milk, expressed as relative contents to the total milk nitrogen; F8: as2-CN, Factor related to the milk as2-CN, expressed as relative content to the total milk nitrogen; F9: as1-CN-Ph, Factor related to the milk as1-CN-Ph expressed as content to the total milk nitrogen; F10: α-LA, Factor related to the milk α-LA.
§ SCS = log2 (SCC / 100 000) + 3. Milk protein fractions: CN, casein; LA, lactalbumin and LG, lactoglobulin. Curd firming: RCTeq, estimated RCT; CFP, asymptotical potential value of CF; kCF, curd-firming instant rate constant; kSR, syneresis instant rate constant; and t max, time at achievement of CFmax. %CY, ratios of the weight (g) of the curd dry matter (%CYSOLIDS) and curd water (%CYWATER) vs. the weight of the processed milk (g); REC, ratio of the weight (g) of the curd constituent (fat, protein or energy, respectively) vs. that of the same constituent in the processed milk (g).
The first factor was heavily and positively loaded on %CYSOLIDS, fat and protein (%) and RECENERGY, thus it was considered as a factor underlying the latent concept of percentage of cheese yield (F1: %CY). The second factor was linked to all CFt traits, but CFP, and the RECFAT, reflecting the curd firmness process (F2: CFt). Due to positive loadings of the instant rate constants (kCF and kSR) and negative to the time traits (RCTeq and t max), the factor was ascribed to a favourable CFt meaning. Moreover, this factor was favourably related to RECFAT. The third factor was associated with the milk, fat and protein daily yields of individual cows, and hence was considered as a descriptor of the milk production (F3: Yield). These 3 factors were almost equally important and together explained 38% of total variability of the 26 milk traits considered. The following 7 factors were less important, each one explaining between 4 to 7% of total variability. The fourth factor was heavily, but negatively, associated with β-LG, the quantitatively most important and variable whey protein fraction, and positively to other N compounds, representing a proxy of casein number, and being positively linked also to RECPROTEIN. Hence this factor was representative of the N found in the cheese (F4: Cheese N). The fifth factor was primary linked to as1-CN (positively) and then to the β-CN (negatively) and was considered as representative of as1- and β-CN importance in milk (F5: as1-β-CN). The sixth factor was positively associated with lactose and negatively with the SCS and the other N compounds. This factor was considered to be an indicator of the udder health status of the cow (F6: Udder health). The seventh factor was strongly and positively associated to the κ-CN and negatively, with a weaker loading, with the β-CN (F7: κ-β-CN). Thus, an increase of this factor was associated to an increased importance of κ-CN. The last 3 Fs were 1 trait-1 factor associations and were named according to the phenotype they were linked to as F8: as2-CN, F9: as1-CN-Ph and F10: α-LA, respectively.
Sources of variation
Table 4 summarises the results of the analysis of variance. The dairy system strongly affected F3: Yield and F1: %CY and had also a negative effect on udder health. The parity of the cow affected four Fs: F3: Yield, F2: CFt, F6: Udder health and also F9: as1-CN-Ph. The days in milk influenced all factors but F4: Cheese N, F9: as1-CN-Ph and F7: κ-β-CN. Factors F3: Yield and F6: Udder health were the only Fs influenced by all effects in the model. On the contrary, none of the effects that the model accounted for affected F4: Cheese N.
Table 4. Analysis of variance (F-values and significance) of the10 factors (F) †
*P < 0·05; **P < 0·01; ***P < 0·001.
† F1: %CY, Factor related to the percentage of individual cheese yield; F2: CFt, Factor related to the curd firmness; F3: Yield, Factor related to the milk yield; F4: Cheese N, Factor related to the milk nitrogen that is present into the cheese curd; F5: as1-β-CN, Factor related to the as1- and β-CN contents in milk, expressed as relative contents to the total milk nitrogen; F6: Udder health, Factor related to the udder health of a cow; F7: κ-β-CN, Factor related to the κ- and β-CN contents in milk, expressed as relative contents to the total milk nitrogen; F8: as2-CN, Factor related to the milk as2-CN, expressed as relative content to the total milk nitrogen; F9: as1-CN-Ph, Factor related to the milk as1-CN-Ph expressed as content to the total milk nitrogen; F10: α-LA, Factor related to the milk α-LA.
Investigating the effect of the dairy system, Table 5 outlines the least squares means (LSM) of the 10 Fs within the 4 different dairy systems and their orthogonal contrasts. Our analysis noted major differences between traditional and modern farms, in favour of the second, for F3: Yield and F1: %CY, and a smaller effect, in favour of the traditional farms, for F6: Udder health. Moreover, for F3: Yield, F1: %CY and F8 as2-CN, the modern farms using TMR showed greater results than those not using TMR. No difference was found in relation to the source of moisture in the TMR.
Table 5. Effects of the dairy system (traditional with tied cows, modern with loose cows), the use of total mixed ration (TMR) within modern farms, and of the moisture source of TMR on the 10 factors (F) †
*P < 0·05; **P < 0·01; ***P < 0·001.
† F1: %CY, Factor related to the percentage of individual cheese yield; F2: CFt, Factor related to the curd firmness; F3: Yield, Factor related to the milk yield; F4: Cheese N, Factor related to the milk nitrogen that is present into the cheese curd; F5: as1-β-CN, Factor related to the as1- and β-CN contents in milk, expressed as relative contents to the total milk nitrogen; F6: Udder health, Factor related to the udder health of a cow; F7: κ-β-CN, Factor related to the κ- and β-CN contents in milk, expressed as relative contents to the total milk nitrogen; F8: as2-CN, Factor related to the milk as2-CN, expressed as relative content to the total milk nitrogen; F9: as1-CN-Ph, Factor related to the milk as1-CN-Ph expressed as content to the total milk nitrogen; F10: α-LA, Factor related to the milk α-LA.
‡ Contrast between the ‘Traditional’ dairy system vs. the three ‘Modern’ ones.
§ Contrast between the ‘Modern No TMR’ dairy system vs. the two ‘Modern TMR’ ones.
¶ Contrast between the ‘Modern TMR Silage’ dairy system vs. the ‘Modern TMR Water’ one.
The LSM of the Fs on which the DIM had an effect are presented in Fig. 2. A ‘lactation like’ pattern (with the peak value in the second month of lactation) was observed for F6: Udder health, and also for F2: α-LA (Fig. 2e, g), revealing, after lactation peak, a decrease in lactose, and an increase of SCS. Lactose was the trait on which the factor F6: Udder health was primary loaded. On the contrary, the latent variable F3: Yield was almost linearly decreased during lactation (Fig. 2c). This pattern is primarily due to daily yield of fat and protein, and not to fresh milk (that peaked in the second month). The strong decrease at the beginning of lactation of milk fat and protein content (and also of their daily yield), paralleled by RECSOLIDS, is reflected by the ‘inverse lactation like’ pattern of the F1: %CY.
Fig. 2. Least squares means of the seven factors across stage of lactation. Description: F1: %CY, Factor related to the percentage of individual cheese yield; F2: CFt, Factor related to the curd firmness; F3: Yield, Factor related to the milk yield; F5: as1-β-CN, Factor related to the as1- and β-CN contents in milk, expressed as relative contents to the total milk nitrogen; F6: Udder health, Factor related to the udder health of a cow; F8: as2-CN, Factor related to the milk as2-CN, expressed as relative content to the total milk nitrogen F10: α-LA, Factor related to the milk α-LA.
A decrease from 1st to 3rd class of DIM was observed for F2: CFt followed by a stabilisation and a slight increase at the end of the lactation (Fig. 2b). An opposite pattern was observed between F5: as1-β-CN and F8: as2-CN (Fig. 2d, f). F8: as2-CN was increasing rapidly at the beginning of lactation when F5: as1-β-CN decreased. A more stable situation was observed thereafter (Fig. 2f).
Figure 3 presents the LSM across the parity levels for the four Fs affected by parity. Factors F3: Yield and F9: as1-CN-Ph had a similar pattern with an increase from the 1st to the 3rd parity, and a slight decline from the 4th parity (Fig. 3b, d). The F2: CFt showed almost no difference between 1st and 3rd parities but an increase after (Fig. 3a). A linear decrease was observed for the factor F6: Udder health related to the mammary gland health status (Fig. 3c).
Fig. 3. Least squares means of the four factors across the levels of parity. Description: F2: CFt, Factor related to the curd firmness; F3: Yield, Factor related to the milk yield; F6: Udder health, Factor related to the udder health of a cow; F9: as1-CN-Ph, Factor related to the milk as1-CN-Ph expressed as content to the total milk nitrogen.
Discussion
Interpretation of the factors
Ten orthogonal latent variables were extracted, out of 26 milk (quantity and quality), milk coagulation and curd firmness indicators, and individual cheese yield traits. The Fs were explaining 74% of the original variability while at the same time drastically reduced the data space. The average KMO was quite close to the value (0·57) reported by Manca et al. (Reference Manca, Serdino, Gaspa, Urgeghe, Ibba, Contu, Fresi and Macciotta2016) using a similar, but smaller, dataset of 12 milk composition, MCP and udder health phenotypes in dairy sheep. Furthermore, the factor model was able to capture basic concepts of the cheese-making process. The same factor scores have been previously used for estimating (co)variance components using standard quantitative genetic model (Dadousis et al. Reference Dadousis, Cipolat-Gotet, Bittante and Cecchinato2017a) and for conducting genome-wide associations and gene-set enrichment analyses (Dadousis et al. Reference Dadousis, Pegolo, Rosa, Bittante and Cecchinato2017b). Results were coherent with the given name of the factors. More precisely, the first four Fs, sorted by variance explained, were able to capture the underlying structure of the cheese yield (%), the curd firmings process, the milk yield and the presence of N into the cheese. Moreover, 4 Fs were associated to the basic milk caseins (as1-β-CN, κ-β-CN, as2-CN and as1-CN-Ph) and 1 factor was related with the, quantitatively, most important whey protein (α-LA). A factor describing the udder health status of a cow, mainly loaded on lactose, other N compounds and SCS, was also obtained.
The meaning of some of the Fs is comparable to previous studies that used MFA, albeit with smaller datasets. For example, in an experimental design study assessing the effect of DIM, body condition and milk protein genotypes on the MCP and the protein composition in Danish Holstein (n = 39) 8 Fs were extracted out of 28 measured variables (Ostersen et al. Reference Ostersen, Foldager and Hermansen1997). Among the measured phenotypes were milk yield traits (milk, fat, protein and lactose), whey and non-casein contents, proportion of caseins, proportion of alleles in genotyped milk proteins, MCP and energy balance and the body condition. The 8 Fs obtained represented MCP properties of milk, milk protein genotypes, the energy and the body condition status of the cows, and effects of stage of lactation on milk yield and protein composition. In another study, eleven milk composition, MCP and udder health phenotypes measured in Brown Swiss cows (n = 1200) were substituted by 4 Fs, describing milk composition, coagulation, acidity and the udder health status (Macciotta et al. Reference Macciotta, Cecchinato, Mele and Bittante2012). In small ruminants, Todaro et al. (Reference Todaro, Scatassa and Giaccone2005) replaced 11 traits related to milk composition, MCP and udder health measured on 117 Girgentana goats with 3 Fs named as ‘slow milks’, ‘milk yield’ and ‘curd firmness’. Recently, Vacca et al. (Reference Vacca, Paschino, Dettori, Bergamaschi, Cipolat-Gotet, Bittante and Pazzola2016), by analysing 9 milk yield, composition and hygiene traits in 1050 Sardinian goats obtained 4 Fs (quality, hygiene, production, and acidity), while Manca et al. (Reference Manca, Serdino, Gaspa, Urgeghe, Ibba, Contu, Fresi and Macciotta2016) in a study on 991 Sarda ewes, extracted 4 Fs (composition and cheese yield, udder health status, coagulation and curd characteristics) from a total of 12 measured phenotypes. All the 5 above mentioned studies used an estimatory MFA with a varimax rotation. The proportion of variance explained by the factor model in those studies ranged between 51·2% (Todaro et al. Reference Todaro, Scatassa and Giaccone2005) to 97% in Ostersen et al. (Reference Ostersen, Foldager and Hermansen1997). Part of this difference, however, could be attributed to the different ways and methods of measuring the phenotypes used in MFA.
Although a direct, factor by factor, comparison between our results and the above mentioned studies is problematic due to differences in the measured phenotypes analysed (number and type), a consistent pattern can be observed: generally, in all the studies MFA clearly distinguishes between milk quality, production, coagulation and health related concepts. These results are encouraging and confirm the ability of MFA to capture basic underlying structure of correlated variables in different species and rearing conditions.
Factor F1: %CY was positively related to fat and protein percentage as well as the %CYSOLIDS and the RECENERGY. This is in agreement with the phenotypic and genetic correlations that have been reported in the literature for these traits (Bittante et al. Reference Bittante, Cipolat-Gotet and Cecchinato2013a) and confirms that cheese produced is not only due to fat and protein content of milk but also to the ability to recover them in curd or lose with whey.
The inverse relationship of RCTeq and t max with kCF, and kSR found on F2: CFt has been previously pointed out in Cecchinato & Bittante (Reference Cecchinato and Bittante2016) and the inclusion of RECFAT on the factor confirm that coagulation properties of milk affect its ability to retain fat in the curd or to lose it with whey. On factor F5: as1-β-CN, the as1- and β- CN were oppositely related. The antagonistic nature of those two caseins has been previously suggested (Bobe et al. Reference Bobe, Beitz, Freeman and Lindberg1999; Bonfatti et al. Reference Bonfatti, Di Martino, Cecchinato, Vicario and Carnier2010). Moreover, the factor F6: Udder health was inversely related to other N compounds that are present in the milk. It is worth noting that within this N fraction, urea was also included. Urea is inversely related to lactose (Miglior et al. Reference Miglior, Sewalem, Jamrozik, Bohmanova, Lefebvre and Moore2007), which was the trait most strongly related to factor F6: Udder health. Herein we considered, for the first time, the milk protein profile along with different milk technological traits. More specifically, a set of 26 relative new phenotypes were considered: 3 production traits (i.e., milk, fat and protein yield), 8 milk protein fractions, 5 traits related to curd firming modelling, 2 cheese yields and 3 nutrients recovery traits. Out of those, only 5 standard milk quality traits were in common with the study of Macciotta et al. (2012) (i.e. fat%, protein %, lactose%, pH and SCC).
To further assess the practical application in dairy farming and in breeding programs of the new set of the latent variables, an investigation on the sources affecting variation on those new traits was followed.
Sources of variation of the factors
Effects of dairy system on the extracted factor scores
Major differences between the different dairy systems were identified for 2 Fs, namely F3: Yield and F1: %CY. Modern farms were associated with a higher milk and solids yield and percentage of cheese obtained for every 100 kg of milk processed. Moreover, within the modern farms, the use of TMR in the feeding system was found beneficial for both latent variables. Only Vacca et al. (Reference Vacca, Paschino, Dettori, Bergamaschi, Cipolat-Gotet, Bittante and Pazzola2016) on lactating goats studied the relationships between Fs and dairy system and farms characteristics: they did not found any associations between Fs and dairy system, but observed an association between altitude of the farm and the production factor and between flock size and the acidity factor.
Nonetheless, the results of our study are in agreement with previous findings on dairy cattle carried out on the measured daily milk yield production traits instead of Fs (Sturaro et al. Reference Sturaro, Cocca, Gallo, Mrad and Ramanzin2009, Reference Sturaro, Marchiori, Cocca, Penasa, Ramanzin and Bittante2013; Bittante et al. Reference Bittante, Cipolat-Gotet, Malchiodi, Sturaro, Tagliapietra, Schiavon and Cecchinato2015). Moreover, although individual cheese traits were not included in the analysis of Bittante et al. (Reference Bittante, Cipolat-Gotet, Malchiodi, Sturaro, Tagliapietra, Schiavon and Cecchinato2015), major differences between the dairy systems were reported for milk fat and protein percentages. These two traits were both loaded on F1: %CY in our study. No significant differences between the dairy systems were observed for the coagulation properties of milk (F2: CFt). The F2: CFt was primary loaded to the kCF, but also to RCTeq, kSR, t max and RECFAT. These findings are consistent with Bittante et al. (Reference Bittante, Cipolat-Gotet, Malchiodi, Sturaro, Tagliapietra, Schiavon and Cecchinato2015), where no significant differences were noted among the dairy systems for kCF and RCTeq. A significant effect of the dairy systems existed in that study for kSR (P < 0·05) and t max (P < 0·01). However, since a factor is a mixture of phenotypes, and in our case F2: CFt was dominated by the kCF, this effect was diluted.
Effect of stage of lactation and parity on the extracted factor scores
The stage of lactation significantly influenced most of the Fs. As expected, the F3: Yield was decreased during lactation, in agreement to Bittante et al. (Reference Bittante, Cipolat-Gotet, Malchiodi, Sturaro, Tagliapietra, Schiavon and Cecchinato2015). The absence of the peak of lactation is probably due to the length of the DIM classes, but it should be noted that, analysing mixed-breed farms in the same area, Stocco et al. (Reference Stocco, Cipolat-Gotet, Bobbo, Cecchinato and Bittante2016) observed the presence of peak of lactation for cows reared in farms characterised by high average milk yield, but not in those with a low milk productivity.
A tendency towards worsening of milk coagulation during lactation is known. This pattern has been reported in Macciotta et al. (Reference Macciotta, Cecchinato, Mele and Bittante2012) using a factor as indicator of milk coagulation (albeit based on the traditional MCP values). Moreover, the same trend has been observed in previous studies using the traditional curd-firmness value at 30 min after the beginning of coagulation, generally known as a30 (Ikonen et al. Reference Ikonen, Morri, Tyrisevä, Ruottinen and Ojala2004). Being consistent, the coagulation ability (F2: CFt) was smoothly decreased (worsening) during early lactation, stabilised between 4th and 8th class of DIM, with an evidence of improvement at the end of the lactation. Moreover, the factor related to the percentage of cheese yield (F1: %CY) showed a decrease from the first to the second month of lactation and then a linear increase. Not surprisingly, this pattern during lactation is consistent with the trend of %CYSOLIDS on which this factor was positively and strongly related (Cipolat-Gotet et al. Reference Cipolat-Gotet, Cecchinato, De Marchi and Bittante2013).
In line with our expectations, the F10: α-LA was decreased during lactation having an opposite trend compared to the F1: %CY. The same pattern has been reported for the measured α-LA in Ostersen et al. (Reference Ostersen, Foldager and Hermansen1997). A smooth decrease of F5: as1-β CN during lactation and the inverse pattern of F8: as2-CN was observed, as expected.
The pattern of the health status of the mammary gland mainly related to the lactose (F6: Udder health) was worsened across the course of lactation. As it is known, lactose in milk is inversely related to mastitis or somatic cell count (Shuster et al. Reference Shuster, Harmon, Jackson and Hemken1991) and SCS has an opposite trend within lactation compared to milk yield (Walsh et al. Reference Walsh, Buckley, Berry, Rath, Pierce, Byrne and Dillon2007). Similar results have been reported using MFA by Macciotta et al. (Reference Macciotta, Cecchinato, Mele and Bittante2012) in dairy cattle and in Manca et al. (Reference Manca, Serdino, Gaspa, Urgeghe, Ibba, Contu, Fresi and Macciotta2016) in dairy sheep. Also in those studies, lactose was the most related trait to the factor underlying the health status of the mammary gland.
The trend of parity effect on F3: Yield was similar to previous results on milk yield traits (Bittante et al. Reference Bittante, Cipolat-Gotet, Malchiodi, Sturaro, Tagliapietra, Schiavon and Cecchinato2015). Moreover, in the same study the kCF values were lower in first lactation cows compared to following lactations, similar to the effect of parity on the F2: CFt found in our analysis. Furthermore, a decline of F6: Udder health with the increase of parity is in agreement with the literature where a decrease in milk lactose percentage and an increase of SCS has been observed in multiparous cows (Yang et al. Reference Yang, Yang, Yi, Pang and Xiong2013). Parity also affected F9: as1-CN-Ph latent variable in our analysis with a pattern among different levels of parities mimicking the trend of F3: Yield.
Conclusions
By applying multivariate factor analysis, 26 original phenotypes were substituted by 10 latent variables. These variables captured important concepts of the cheese-making process as well as basic bovine health indicators related to the mammary gland. Results of ANOVA were in agreement to the given name of the factor and reflected the underlying structure that each factor was representing. This approach, therefore, represents a valuable tool for studying the effects of different production systems, feeding regimes, and health status on the milk protein composition and cheese-making related traits, and for identifying strategies for altering such technological traits in accordance with dairy industry.
The authors wish to thank the Trento Province (Italy), the Italian Brown Swiss Cattle Breeders Association (ANARB, Verona, Italy) and the Superbrown Consortium of Bolzano and Trento for financial and technical support. C. Dadousis benefitted from financial support of the CARIPARO (Cassa di Risparmio di Padova e Rovigo) Foundation (Padua, Italy). The authors also wish to acknowledge Valentina Bonfatti for her cooperation in assessing the milk protein fractions.
