The use of milk Fourier-transform infrared spectra for predicting cheesemaking traits in Grana Padano Protected Designation of Origin cheese

The prediction of the cheese yield (%CY) traits for curd, solids, and retained water and the amount of fat, protein, solids, and energy recovered from the milk into the curd (%REC) by Bayesian models, using Fourier-transform infrared spectroscopy (FTIR), can be of significant economic interest to the dairy industry and can contribute to the improvement of the cheese process efficiency. The yields give a quantitative measure of the ratio between weights of the input and output of the process, whereas the nutrient recovery allows to assess the quantitative transfer of a component from milk to cheese (expressed in % of the initial weight). The aims of this study were: (1) to investigate the feasibility of using bulk milk spectra to predict %CY and %REC traits, and (2) to quantify the effect of the dairy industry and the contribution of single-spectrum wavelengths on the prediction accuracy of these traits using vat milk samples destined to the production of Grana Padano Protected Designation of Origin cheese. Information from 72 cheesemaking days (in total, 216 vats) from 3 dairy industries were collected. For each vat, the milk was weighed and analyzed for composition (total solids [TS], lactose, protein, and fat). After 48 h from cheesemaking, each cheese was weighed, and the resulting whey was sampled for composition as well (TS, lactose, protein, and fat). Two spectra from each milk sample were collected in the range between 5,011 and 925 cm −1 and averaged before the data analysis. The calibration models were developed via a Bayesian approach by using the BGLR (Bayesian Generalized Linear Regression) package of R software. The performance of the models was assessed by the coefficient of determination (R 2VAL ) and the root mean squared error (RMSE VAL ) of validation. Random cross-validation (CVL) was applied [80% calibration and 20% validation set] with 10 replicates. Then, a stratified cross-validation (SCV) was performed to assess the effect of the dairy industry on prediction accuracy. The study was repeated using a selection of informative wavelengths to assess the necessity of using whole spectra to optimize prediction accuracy. Results showed the feasibility of using FTIR spectra and Bayesian models to predict cheesemaking traits. The R 2VAL values obtained with the CVL procedure were promising in particular for the %CY and %REC for protein, ranging from 0.44 to 0.66 with very low RMSE VAL (from 0.16 to 0.53). Prediction accuracy obtained with the SCV was strongly influenced by the dairy factory industry. The general low values gained with the SCV do not permit a practical application of this approach, but they highlight the importance of building calibration models with a dataset covering the largest possible sample variability. This study also demonstrated that the use of the full FTIR spectra may be redundant for the prediction of the cheesemaking traits and that a specific selection of the most informative wavelengths led to improved prediction accuracy. This could lead to the development of dedicated spectrometers using selected wavelengths with built-in calibrations for the online prediction of these innovative traits.


INTRODUCTION
The manufacture of cheese is one of the main profitable ways for adding value to milk.In Italy, more than 80% of the milk is processed into cheese, of which ~50% are Protected Designation of Origin (PDO) cheeses, with the remaining portion mostly used to produce local and typical cheeses.Because of the tight regulations imposed on the origin of the milk and production methods, the cheeses with labels have distinct qualities.Indeed, their quality highly depends on the farming system, geographical area, and season (Bertoni et al., 2001;Bergamaschi et al., 2020).The PDO restrictions increase cheese manufacture costs along the dairy chain.However, due to the interaction of human and environmental elements that are distinctive to a particular location, these productions present the potential to obtain a higher price (Bouamra-Mechemache and Chaaban, 2010).A representative example of an Italian PDO product is Grana Padano, an artisanal, cooked, extra-hard, long-ripened cheese from partially skimmed raw milk produced in several provinces of the Po Valley (Northern Italy).The several steps involved in producing this cheese have been refined over the centuries and, since 1996, are regulated by the European Union.Nowadays, Grana Padano has the highest production rate of any PDO cheese in the European Union (CLAL, 2020) and is therefore a suitable case study to represent the labeled products with similar characteristics.
The analysis of milk composition is necessary for PDO cheese production to monitor the production chain and also to use that information to establish milk price for the farmers (Galgano et al., 2001).The monitoring of cheesemaking efficiency is traditionally assessed by using milk-based predictive formulas (Emmons and Modler, 2010) to estimate cheese yield (%CY), the most important trait in the economics and profitability of dairy industries.Unfortunately, the use of milk quality traits alone cannot offer a realistic prediction of the overall cheesemaking process, as the efficiency of the cheesemaking is defined not only by %CY, but also by the recovery of individual milk constituents in the curd (%REC), and their loss in the whey (Banks, 2007;Cipolat-Gotet et al., 2013).However, information on the rate of milk nutrient recovery in the curd and, in general, on the cheesemaking efficiency of commercial cheeses is very scarce.This is essentially due to cost limitations for sample collection and analysis under operational conditions.Among the available opportunities to fill this gap, infrared (IR) spectroscopy has become broadly used along the dairy chain for the determination of the most varied composition of milk (Wojciechowski and Barbano, 2016;Bresolin and Dórea, 2020;Tiplady et al., 2020) including the prediction of cheesemaking traits (Ferragina et al., 2013) and the composition of different dairy products (Woodcock et al., 2008;Stocco et al., 2019;Bittante et al., 2022).The ever-growing development of new user-friendly spectrometers with built-in calibration (De Marchi et al., 2014) and the novel approaches for the treatment of the milk spectra (Tiplady et al., 2019) provide a new opportunity to get quick and accurate results to monitor the cheesemaking process.Fourier-Transform IR (FTIR) spectroscopy has often been used to predict the cheesemaking traits from individual milk samples (Ferragina et al., 2013;El Jabri et al., 2019;Stocco et al., 2023).Fagan et al. (2008) conducted a study at the production level to predict cheesemaking traits using near-infrared sensors installed in the vats, along with other relevant traits such as cheese curd syneresis.Nevertheless, no study has been published yet on the prediction of different measures of %CY and %REC traits from milk FTIR spectra collected in the vat and used to produce Grana Padano PDO.Therefore, the aims of this study were to (1) characterize %CY and REC in the Grana Padano production; (2) investigate the feasibility of using FTIR bulk milk spectra to predict %CY and %REC traits, (3) quantify the effect of the dairy industry on prediction accuracy, iv) explore the contribution of the milk spectrum for the prediction of these cheesemaking traits.

MATERIALS AND METHODS
No human or animal subjects were used, so this analysis did not require approval by an Institutional Animal Care and Use Committee or Institutional Review Board.

Origin and Collection of Samples and Spectra
A total of 216 milk and whey samples (aliquot of 50 mL per each sample) were collected from the cheesemaking vats during the mixing of the milk, before the start of the cheesemaking process while whey samples were collected after the separation of the curd.Sampling was performed between February 2021 and January 2022 from 3 cheese industries in the north of Italy, comprised within the Grana Padano PDO cheese production area.The cheesemaking process complied with the production regulation stated by the Consortium (Reg.EU 584/2011).Samples were analyzed for chemical composition at the Breeders Association of Lombardy Region (Italy) by using a FTIR spectrometer (MilkoScan FT 6000, Foss Electric, Hillerød, Denmark), over the spectral range from 5,011 to 925 cm −1 .Spectra were stored as absorbance (A) using the transformation A = log(1/T), where T is the transmission (Figure 1).Two spectral acquisitions were carried out for each sample from the same vial, and they were averaged after visual editing for outlier detection before data analysis.Somatic cell score [calculated as log 2 (SCC/100,000) + 3], differential SCC (DSCC), and logarithmic bacterial count [LBC; log 10 (total bacterial count/1,000)] were also measured by the breeder association as innovative traits to monitor the health status of the cows in PDO cheese (Jayarao et al., 2004).Averaged infrared spectra of vat milk over the spectral range from 5,011.5 to 925.9 cm −1 , with intervals from minimum to maximum (A, without any pretreatment; B, after standardization).Spectral regions are identified as short-wavelength infrared (SWIR), short-and mid-wavelength infrared (SWIR-MWIR), first and second mid-wavelength infrared (MWIR-1, MWIR-2, respectively), and mid-and longwavelength infrared (MWIR-LWIR).

Grana Padano PDO Protocol of Production
The milk to produce Grana Padano PDO is collected from herds located in several provinces of the Po Valley.The raw milk is stored at a temperature of not less than 8°C and then is skimmed via natural creaming before being homogeneously mixed in 1,000-L vats.Vat milk is added of an aliquot of the whey from the day before is added as a starter to decrease its pH and then heated between 31 and 33°C.When this temperature is reached, the calf rennet is added.After the coagulation of the milk, the curd is cut into small granules to improve the expulsion and separation of the whey.The curd is then cooked until 56°C to promote further whey expulsion.The curd is left to rest in the vat, immersed in the whey, for a maximum of 70 min from the end of the cooking phase, to form a compact mass.Then, the curd is divided in 2 portions, extracted from the vat and placed into cylindrical molds.During this phase, the cheese is pressed to remove additional whey and compact the mass.After 2 d in the molds, the cheese is dry-salted or immersed in brine.This step can take from 14 to 30 d, depending on the saline solution and the size of the cheese wheel.After brining, the cheese wheels are transferred to the ripening rooms.Grana Padano must be ripened for at least 9 mo, and some varieties can be aged for up to 20 mo or more.

Definition of Cheesemaking Traits
The weights of the milk and the curd after brining were used to calculate the weight of the whey, whereas the weight of each curd component was calculated by subtracting the weight of that component in the whey from the corresponding one in the processed milk.This information allowed us to calculate 7 cheesemaking traits, as proposed by Cipolat-Gotet et al. (2013).The formulas are reported as follows: 1. %CY CURD = [weight of fresh curd (g)/weight of milk processed (g)] × 100; the percentage of the weight of fresh curd compared with the weight of the processed milk.2. %CY SOLIDS : [weight of curd TS (g)/weight of milk processed (g)] × 100; the percentage of the weight of curd dry matter compared with the weight of the processed milk.3. %CY WATER : [weight of curd moisture (g)/weight of milk processed (g)] × 100; the percentage of the weight of curd water compared with the weight of the processed milk.4. %REC PROTEIN : [weight of curd protein (g)/protein content of milk (g)] × 100; the percentage of the weight of the curd's protein compared with the protein content of the milk.Any potential outlier of each trait exceeding 3 standard deviations (SD) from the mean were excluded from the analysis of the measured traits to ensure accuracy.

Editing of the Spectra and Chemometric Models
Before spectra analysis, the raw absorbance values of each wavelength in the FTIR spectra (Figure 1a) were centered and scaled to a null mean and a unit variance (Figure 1b).Then, samples having a large spectral distance (i.e., Mahalanobis distance > 3) were considered as outliers and removed from the calibration dataset.Four samples were identified as outliers and removed from the data.The spectra were not subjected to any other mathematical pretreatment.The calibration models were developed using a Bayesian approach (Ferragina et al., 2015), in particular using the Bayesian Generalized Linear Regression (BGLR; de los Campos and Perez Rodriguez, 2015) package available in R software (R Core Team, 2021).Each trait was regressed to 1,060 wavelengths using y = μ + ∑ 1,060 j x ij β j + e i , where μ is the overall mean, x ij are the FTIR wavelengths, i is the sample (from 1 to 216) and j the wavelengths (1,060 from 5,011 to 925 cm −1 ), β j are the regression coefficients, and e i is the residual assumed to be independently and identically distributed with a normal distribution with mean equal to 0 and variance equal to s 2 e .The Bayes B model implemented in the package we used incorporates prior information about the model parameters and updates this information based on the observed data to estimate the effects of wavelengths on the phenotypes (Ferragina et al., 2015).

CVL Models
To investigate the feasibility of using FTIR bulk milk spectra to predict %CY and %REC traits, a

Molle et al.: PREDICTION OF CHEESEMAKING TRAITS
random cross-validation (CVL) was applied, in which 80% of the total records were randomly selected and used to build the equation (calibration set; CAL), and the remaining 20% of records were used to test the model (validation set; VAL).To account for sample variability, the procedure was repeated 10 times for each trait and the results were averaged over the 10 replicates.The SD across the 10 replicates was also calculated.To quantify the dairy plant effect on the model's predictive ability, a stratified cross-validation (SCV) procedure was applied, developing the model using the samples of 2 plants and validating it on the samples of the third plant (procedure repeated 3 times, each dairy plant representing one-third of the dataset).The coefficient of determination (R 2 VAL ) and the root mean squared error of validation (RMSE VAL ) were used to assess the models performances.For each trait, (1) the equation coefficients corresponding to each wavelength (n = 1,060), and (2) the Pearson correlations between the trait and the absorbances at each wavelength were assessed.The results were plotted dividing the overall spectral range into 5 spectral regions, as previously presented by Ferragina et al. (2017): short-wave infrared (SWIR, often reported in the literature as near-infrared), the transition between short-and mid-infrared (SWIR-MWIR), first and second mid-infrared (MWIR1 and MWIR2, respectively), and transition between mid-and long-infrared spectral regions(MWIR-LWIR).
To assess the importance of single wavelengths for the prediction of the cheesemaking traits and the effect of the wavelength selection on the prediction accuracy, the wavelengths have been ranked by their Bayes β coefficient for each trait.A CVL (80% CAL-20% VAL sets randomly assigned; 10 replicates) was applied for each trait, using 10 different subsets containing the first 1,2,5,10,25,50,100,200,400, and 800 wavelengths, respectively.As for the previous CVL procedure, the averages of R 2 VAL and RMSE VAL values derived from the 10 replicates were used as indicators of prediction accuracy to identify an optimal selection of the number of wavelengths compared with the prediction statistics obtained using the entire spectrum.

RESULTS
Table 1 reports the descriptive statistics of vat milk composition analyzed for this study.Milk fat in the vat was ~2.64%, a percentage equal to that of casein (2.63%), leading to a fat: casein ratio equal to 1.00 (the regulation for Grana Padano PDO establishes this ratio between 0.80 and 1.05).The SCS was on average 0.53 (corresponding to ~18,000 cells/mL of milk) with a DSCC of 62% and an average 1.85 for LBC (corresponding to ~71,000 colony-forming bacterial units/ mL of milk).
Table 2 reports the descriptive and prediction statistics of the cheesemaking traits achieved via CVL by using vat milk spectra.In particular, %CY SOLIDS and %CY WATER were on average 4.72% and 3.86%, respectively, giving a %CY CURD of ~8.54%.Regarding %REC traits, the average values were ~75%, 85%, 42%, and 57% respectively for %REC PROTEIN , %REC FAT , %REC SOLIDS , and %REC ENERGY .The prediction statistics showed R 2 CAL values ranging from 0.59 to 0.75 for the %CY traits and from 0.39 to 0.63 for the %REC traits.The %CY CURD had a R 2 VAL value of 0.54 and a RMSE VAL of 0.30% while the SD of the trait was 0.44%.In the present study, %REC PROTEIN appeared to be the most predictable among %REC traits, with R 2 VAL of 0.48 and RMSE VAL of 0.53% while %REC FAT exhibited lower prediction accuracy.Table 3 shows the descriptive and prediction statistics for the SCV across the 3 dairy plants.Dairy A had the highest percentage of all 3%CY traits, followed by dairy C, and then dairy B. Dairy C had the highest %REC PROTEIN (75.34% ± 1.13), followed by dairy A (74.87% ± 0.28), and then dairy B (74.56% ± 0.48).Dairy B had the highest percentage of milk %REC FAT (88.55% ± 1.44), followed by dairy A (83.89% ± 1.90), and then dairy C (83.35% ± 2.15), whereas dairy A showed the highest percentage of %REC SOLIDS (42.37% ± 0.88), followed by dairy C (41.79% ± 2.10), and dairy B (41.30% ± 1.16).Dairy A had also the highest %RE-C ENERGY (56.67% ± 0.99), followed by dairy B (56.76% ± 0.82), and dairy C (55.89% ± 1.90).In the SCV, the R 2 CAL were relatively high for %CY CURD , %CY SOLIDS , %REC PROTEIN , and %REC ENERGY .The R 2 VAL in the SCV ranged from 0.05 to 0.34 and the RMSE VAL was generally higher than RMSE CAL , ranging from 0.20 to 3.99.

Optimum Number of Wavelengths for the Prediction of Cheesemaking Traits
By reducing the number of the spectrum wavelengths to only those with the highest absolute β values, it is possible to reduce the quantity of data required for the predictive model with a low or null effect on the accuracy and repeatability of the prediction.Figure 2 illustrates the relationship existing between the number of spectrum wavelengths used for the prediction of cheesemaking traits and the performance of the equation models in terms of R 2 VAL and RMSE VAL .To achieve this objective, the individual wavelengths were selected by ranking the most important Bayes β coefficients in absolute value for each cheesemaking trait.As expected, the R 2 VAL increased and the RMSE VAL decreased with an increasing number of informative wavelengths used for the development of the predictive model.However, the maximum performance of the model was achieved with a different and relatively small number of wavelengths according to the examined trait.Moreover, the patterns of R 2 VAL and RMSE VAL among cheesemaking traits were different (Figure 2).For example, for %CY SOLIDS 10 wavelengths were sufficient to achieve the R 2 VAL and RMSE VAL comparable to that obtained using the entire spectrum.In opposite, 200 and 100 wavelengths were needed to reach the maximum prediction accuracy for %CY CURD and %CY WATER , respectively.In the case of %CY CURD and %CY WATER, the prediction accuracy gradually improved to 200 and 100 wavelengths, respectively, but then additional wavelengths reduced the prediction accuracy (Figure 2).Similar scenarios were observed also for the %REC traits, even though the reduction of the model  performance was at a different number of wavelengths, according to the trait considered.
Another interesting finding was that among the first 10 most informative wavelengths, different for each trait, the largest part of them belonged to the water regions of the spectrum (Supplemental Table S1, https: / / doi .org/ 10 .6084/m9 .figshare.25075805,Stocco, 2024c), further confirming the importance of maintaining the entire spectrum water range for the prediction of these traits, before selecting the most informative wavelengths.For example, it was found that the wavelengths ~3,470 cm −1 were highly informative for the prediction of %CY WATER and %REC FAT .The selection of the most important wavelengths showed that most of the peaks with the highest Bayes B are located in the water region, suggesting their usefulness for an accurate prediction of cheesemaking traits.These findings clearly demonstrate that the usual removal of water peaks from the calibration may cause the loss of significant chemical information.The wavelengths usually attributed to fat had also large importance, as 3 of them were present in the first 10 wavelengths for the prediction of %CY CURD , and 2 for %CY SOLIDS and %CY WATER .Moreover, among the bands usually attributed to protein, ~1,100 and between 3,100 and 3,300 cm −1 , 3 wavelengths for %CY CURD , 4 for %REC FAT, and one for %REC PROTEIN were in the top 10 wavelengths, respectively.However, it is important to consider that the first 10 wavelengths in absolute value for each trait were not always sufficient to ensure the best prediction accuracy, as shown in Figure 2.This was especially clear for the cheesemaking traits for which 100 to 200 wavelengths were fundamental to reach high R 2 VAL and low RMSE VAL .According to these findings, using only specific wavelengths of the FITR milk spectrum could guarantee sufficient predictive ability for the indirect measure of the cheesemaking traits in field conditions for Grana Padano PDO production.This can significantly reduce the time required for data analysis and processing and would also open the possibility of developing specialized portable and low-cost mid-infrared spectrophotometers characterized by a low number of measurable wavelength devices to monitor informative traits during the cheesemaking process.

Milk Composition and Measured Cheesemaking Traits
Despite the abundance of cheese varieties, only some of them are recognized and largely exported worldwide.Among these, Grana Padano represents, alone, more than 20% of the milk use and one-third of the total  Italian PDO production (CLAL, 2020).Several studies have been carried out with the aim to offer tools for the monitoring of the process, from the mimicking of the coagulation pattern to the description of milk microbiota (Rossetti et al., 2008;Pretto et al., 2013;Stocco et al., 2015), but no study has been done on the prediction via FTIR spectroscopy of the cheesemaking traits in this type of cheese.One of the most peculiar aspects of the production of Grana Padano cheese is represented by the natural creaming of the milk before cheesemaking, as seen by the relatively low-fat content and its low variability (Table 1).This practice derives from the old tradition of villages where it was important to recover milk fat for nutritional purposes (Gobbetti et al., 2018).Today the skimming practice is used to standardize and reduce the fat content in milk and, depending on the temperature and duration, it helps to increase the autochthonous mesophilic microbiota.The values for SCS and LBC in this study are generally low compared with raw full-fat milk (Jayarao et al., 2004: Abeni et al., 2005).Indeed, natural creaming has a cleaning effect, consisting in the ascend of bacteria, spores, and somatic cells to the surface of the milk together with fat globules (Caplan et al., 2013).Unlike SCS and LBC, DSCC seemed not to be influenced by this effect, as its value was in line with previous studies reporting DSCC content in raw full-fat milk from individual cows (Bobbo et al., 2019;Stocco et al., 2020).

Prediction of Cheesemaking Traits via CVL Procedure
The average for cheesemaking traits measured in this study was different from those usually observed in previous studies based on laboratory cheesemaking procedures (Cipolat-Gotet et al., 2013;Ferragina et al., 2013;Bittante and Cipolat-Gotet, 2018).The measurements used in the current investigation were obtained from a process that was carried out in accordance with the PDO manufacturing specifications for Grana Padano.Specifically, for %CY CURD and %CY WATER , the proportion of water to the total fresh curd was lower compared with that of dry matter, and this is in contrast with most of the studies in the literature (Bittante et al., 2014;Cipolat-Gotet et al., 2016;Stocco et al., 2017).Grana Padano PDO can be classified as low-moisture cheese because of the dimension of curd particles (size of rice grains) after cutting, the intensity of the heat treatment during curd cooking, and curd-handling techniques (i.e., curd extraction).These procedures that characterize the cheesemaking of Grana Padano induce a higher loss of water in the whey (Gobbetti et al., 2018).Similar to previous studies (Cipolat-Gotet et al., 2013), the variability of %CY WATER (coefficient of variability = 8%) was the highest among the %CY and %REC traits, because the moisture retention in the curd is the most difficult trait to control at both the laboratory and dairy industry levels.Indeed, water retention in cheese varies and it is influenced mainly by processing conditions, such as the type and concentration of rennet, the cutting time and intensity, the draining and pressing of wheels, the salting technique, the length, and the climatic conditions of ripening (Remeuf et al., 1991;Janhøj and Qvist, 2010;Everard et al., 2011).Proteins play a fundamental role in the coagulation and syneresis processes that characterize cheesemaking (Emmons et al., 2003), and the loss of proteins in whey reduces the %CY (Hallén et al., 2010).Almost all the predictive formulas for estimating %CY have been based on knowledge of the protein (or casein or para-caseinate) and fat contents of milk (van Slyke and Price, 1952;Banks et al., 1981;Emmons et al., 1990) or the sum of them (Verdier-Metz et al., 2001).All these formulas assume that the recovery of milk protein and fat in the curd is constant.However, previous findings reported that these recoveries are not linearly dependent on the chemical composition, but rather on a composite effect of the chemical composition, technological properties of the milk, and the different cheesemaking steps involved in the process, such as the time and size of curd cutting (Janhøj and Qvist, 2010;Cipolat-Gotet et al., 2013).The lower average of both %REC PROTEIN and %REC FAT (Table 2) found here were related to the milk fat: casein in the vat, the curd cutting (i.e., higher losses of fat in the whey), and the curd cooking temperature (i.e., higher protein denaturation) in the Grana Padano PDO production.
Moving to the prediction statistics, for all %CY and %REC, the R 2 CAL was overall lower than in other studies related to the prediction of the same traits collected on individual milk samples (Ferragina et al., 2015(Ferragina et al., , 2017;;Stocco et al., 2023).Due to the small number of examined cheese factories and individual vats compared with the trials created at the laboratory level, the collection of phenotypes under field conditions implies low variability in the dataset.Furthermore, the partial skimming of the milk in the Grana Padano PDO leads to a standardization of the chemical composition among the vats, thus reducing the overall original variability of milk composition and spectra.With a residual predictive deviation (RPD; as the ratio between SD and RMSE VAL ) higher than 1, the calibration for %CY CURD allows to achieve a prediction with a smaller error than the own variability of the trait.This outcome is encouraging and could be further improved with a larger and more variable dataset, including data from multiple cheese industries, to enhance the overall variability of the training set and enable the development of a more robust model.The inclusion of 3 cheese industries aimed at achieving an adequate variability of the dataset.In contrast, the 3 different cheese industries differed in the chemical composition and the operative procedures, both of which affected the prediction accuracy of the SCV (Table 3).The combination of intragroup low variability and large intergroup differences can have an effect on the robustness of the calibration.In particular, the variability of the calibration set is one of the most important factors affecting the prediction accuracy (Grelet et al., 2021).Moreover, it is commonly accepted that complex phenotypes such as cheesemaking traits are predicted by being associated with signals in the spectrum associated with components having a relationship with the targeted trait (i.e., milk fat and protein).
The R 2 CAL values indicate that the calibration is sufficiently robust.Our results for the prediction (in term of R 2 VAL and RMSE VAL ) are similar to those reported in previous findings in which cheesemaking traits in individual bovine milk samples were predicted via CVL.For instance, Ferragina et al. (2015) found R 2 VAL of 0.71, 0.65, and 0.28 and RMSE VAL values of 1.03, 1.44, and 3.13 for %CY CURD , %REC PROTEIN , and %REC FAT , respectively, using a similar statistical model and CVL procedure in Brown Swiss milk samples.Similarly, Bittante and Cipolat-Gotet (2018) obtained R 2 VAL and RMSE VAL values of 0.79 and 0.43, respectively, for predicting %CY SOLIDS using a similar dataset.El Jabri et al. ( 2019) used milk FTIR spectra to compare the use of Bayesian linear regression and partial least square (PLS) regression for predicting cheesemaking traits from Montbéliarde cows, with the best prediction models achieving R 2 VAL values of 0.85 and 0.91 for %CY CURD and %CY SOLIDS , respectively.El Jabri et al. (2019) In that specific study, the water regions of the milk spectrum were removed before applying the PLS models.This step may help to improve the R 2 VAL by removing the noisy spectral range which contains most of the variability between the samples, allowing the prediction model to fit better the data.However, the removal of those regions may lead to the loss of some of the useful information, as stated by Wang and Bovenhuis (2018), and as also shown in the Supplemental Figure S1 (https: / / doi .org/ 10 .6084/m9 .figshare.25075904.v1,Stocco, 2024a) by the high Bayes B peaks falling in these water regions.In our study, we calculated also the RPD as a metric of model validity to compare models across different studies.In our study, the RPD for the CVL ranged from 1.11 for the trait predicted with the lowest accuracy to 1.72.The RPD presented in Table 2 demonstrates the significant role of water retention in %CY traits and the influence of its variability on prediction accuracy.Indeed, the %CY SOLIDS values were predicted with the highest R 2 VAL (0.66) and the lowest RMSE VAL (0.16%), whereas the %CY WATER values displayed the lowest predictive performances (R 2 VAL = 0.44 and RM-SE VAL = 0.23%).This lower prediction accuracy for %CY WATER agrees with other studies investigating the prediction of cheesemaking traits in milk from different species.For example, Stocco et al. (2023) showed R 2 VAL of 0.64 and RMSE VAL of 1.43 for %CY WATER in caprine species, having the RMSE VAL value 3 times higher more than that of %CY SOLIDS .
All the %REC traits were characterized by lower prediction accuracies compared with %CY in this study and also if compared with other studies (Ferragina et al., 2013;El Jabri et al., 2019).In general, %REC traits are highly influenced by the relationship among milk components, coagulation properties, and operating conditions, rather than the concentrations of individual components in the milk.As an example, the ability of the curd to retain fat is mostly influenced by the strength of the protein network rather than the fat composition of the milk itself (Bittante et al., 2013).This can also be seen in the spectra, where the wavelengths usually associated with protein (between 3,100 and 3,300 cm −1 ), as described by Stuart, (2004), showed a relatively high peak of Bayes B for the prediction of %REC FAT .Moreover, %REC FAT is largely influenced by the time of curd cutting, the dimension of the particle curd size, and the cooking temperature (Janhøj and Qvist, 2010).These factors cannot be reflected by the milk spectrum, but they largely affect the variability of %REC FAT , and therefore its prediction accuracy.For these reasons, %REC FAT had the highest RMSE VAL (2.46%, Table 2).The variability of %REC FAT in turn contributes to the variability of %REC SOLIDS and %REC ENERGY , which were characterized by unsatisfactory prediction accuracies as well (Table 2).Among all traits, the ones that depend on others (e.g., %CY CURD affected by %CY WATER and %CY SOLIDS ; %REC SOLIDS and %REC ENERGY affected by %REC PROTEIN and %REC FAT ) showed a higher variability of the slope of the equation, probably due to the combinations of more factors affecting the variability of the prediction accuracy.

Prediction of Cheesemaking Traits from Milk Spectra via SCV Procedure
To evaluate the effect of the dairy plant on the prediction accuracy of cheesemaking traits, an SCV procedure was tested (Table 3).The milk processed by various dairy facilities adhering to the PDO cheese manufacturing requirements, as in the case of Grana Padano, is expected to adhere to a defined processing technique.However, Grana Padano is still an artisanal cheese, and many factors, from milk quality to operating conditions, could inevitably differ between plants, thus influencing cheesemaking efficiency.The variability among the dairy industries could be observed from the descriptive statistics in Table 3, where %REC FAT exhibited the highest SD values for 2 industries (A and B), whereas for dairy C the traits related to the milk TS were the most variable.These differences affected the prediction accuracy of the cheesemaking traits, leading to a relatively low R 2 VAL and high RMSE VAL (Table 3).To illustrate this situation, the variations in cheesemaking traits among different dairy industries were evaluated using a Tukey test (data not displayed).For instance, the B and C dairy factories exhibited significant differences in all the analyzed traits, which were indicative of diverse prediction scenarios, as depicted in Supplemental Figure S2 (https: / / doi .org/ 10 .6084/m9 .figshare.25075910.v1,Stocco, 2024b).In addition to the variability of the studied traits caused by the operative conditions, the poor performance of the validation among dairy plants may be also due to the low intrinsic variability of the dataset caused by the natural creaming undergone by the bulk milk before its use in the vats.
As expected, the SCV highlights higher overfitting (difference among R 2 CAL and R 2 VAL ) for all the traits when compared with the CVL results.The poorest predictions were obtained for %CY WATER and %REC FAT (Table 3), for the same reasons as for the CVL, underling different operative procedures (e.g., cutting and temperature), directly influencing the curd firmness, and then the %REC FAT and %CY WATER.These differences can be seen in Supplemental Figure S1 where the cheese industries A and B showed different patterns of prediction bias.The dairy industry C showed an average higher variability for all the cheesemaking traits, which improved the overall prediction accuracy for the 2 other dairy plants when C was included in a calibration set.Those considerations insist on the fact that the broad variability found in the samples is needed to build robust calibrations.
The results of SCV revealed that it is challenging to accurately predict these traits, as demonstrated by the low RPD.To enhance the accuracy of the prediction model, it may be beneficial to increase the number of dairy industries included in the training set and the number of samples for each dairy plant, covering as much variability as possible.This could lead to the development of a more robust model to handle the variability of the data in the VAL set.

CONCLUSIONS
The promising results of this study demonstrated that the %CY and %REC traits in the Grana Padano PDO production can be predicted using the FTIR spectra of vat milk before cheesemaking.This finding presents the opportunity for a new range of quality monitoring tools to be used at the field conditions to monitor the efficiency of cheese manufacture instead of predictive formulas based on milk composition.The SCV allowed us to test the importance of the dairy plant effect factors on the prediction accuracy of the cheesemaking traits.This study also shows that the need to manage vast amounts of data might be reduced by choosing the most significant wavelengths in a FTIR spectrum to guarantee correct predictions for several cheesemaking features.

2
Interval = minimum and maximum of observed values.3 N CAL = number of items in the calibration set; R 2 CAL and R 2 VAL = coefficient of determination for the calibration and the validation ± SD; RMSE CAL and RMSE VAL = root mean square error for the calibration and validation ± SD; slope of the calibration equation ± SD; RPD = residual predictive deviation, ratio between SD and RMSE VAL ± SD; prediction bias = the average difference between the predictions and the labels in the dataset in absolute value ± SD.

Figure 2 .
Figure 2. Evolution of the accuracy of prediction at the increasing of the number of wavelengths used to build the calibration for each cheesemaking trait.Cheese yield (%CY; the weight of fresh curd, curd solids, and curd water as a percentage of the weight of milk processed); milk nutrients recovery in the curd (%REC; the weight of the curd component [protein, fat, TS, energy] to the same component in milk, multiplied by 100.R 2 VAL = coefficient of determination for the validation set; RMSE VAL = root mean square error for the validation set.

Table 1 .
Molle et al.: PREDICTION OF CHEESEMAKING TRAITS Descriptive statistics of vat milk quality 1 2DSCC = differential SCC.

Table 2 .
Molle et al.:PREDICTION OF CHEESEMAKING TRAITS Descriptive statistics of cheesemaking traits 1 obtained from the dairy cheese factories and prediction statistics deriving from the cross-validation procedure using FTIR %CY; the weight of fresh curd, curd solids, and curd water as a percentage of the weight of milk processed); milk nutrients recovery in the curd (%REC; the weight of the curd component [protein, fat, TS, energy] to the same component in milk, multiplied by 100).

Table 3 .
Molle et al.:PREDICTION OF CHEESEMAKING TRAITS Descriptive statistics (mean ± SD) of cheesemaking traits 1 obtained from the dairy cheese factories and prediction statistics (mean ± SD) deriving from the stratified cross-validation procedure using Fourier-transform infrared spectra from vat milk samples Molle et al.: PREDICTION OF CHEESEMAKING TRAITS