The milk fingerprint of Sardinian dairy sheep: quality and yield of milk used for Pecorino Romano P.D.O. cheese production on population-based 5-year survey

The Pecorino Romano P.D.O. is the main sheep cheese produced in Italy and the first one among the sheep cheeses, in terms of quality and value, exported from the European Union. About half of the sheep milk produced in Italy is processed into this type of cheese by 36 dairies belonging to the Pecorino Romano Consortium. Eight million records of biweekly analyses of milk collected within a 5-year period from farms delivering their milk to the aforementioned consortium were analysed in this work. Monthly evolution curves were plotted for fat, protein, lactose, pH, NaCl, SCC, bacterial load and principal fatty acids (FAs). Due to the seasonal production systems of Sardinian sheep, monthly evolution of milk fat and protein contents and cheese yield are directly linked to the lactation curve pattern and the pastures quantity and quality. Also, the FA profile of milk is affected by grass availability and quality in both early and mid-lactation, whereas it is influenced by the energy balance of ewes in late lactation. Cheese yield equation was computed based on fat and protein contents and considering the variability among dairies in technological processes used in transforming Sarda sheep milk to Pecorino Romano P.D.O. These data could be a relevant basis to set-up future grids of milk payments based on quality standards. Moreover, they could be useful to formulate administrative policies on the dairy sector with the prospective to improve milk quality of Sardinian sheep destined to the Pecorino Romano production.


Introduction
Based on the official website of its Consortium of Protection (Consorzio per la Tutela del Formaggio Pecorino Romano DOP 2020), the Pecorino Romano P.D.O. (PR) is described as 'a hard cheese, made with fresh whole sheep milk, derived exclusively from farms located in the area of production (Sardinia, Lazio and the Grosseto province in Tuscany)'. It has a thin crust with the colour of pale ivory or straw, and sometimes it is covered by a special protective material. It is an aromatic cheese with a lightly spicy and tangy taste in geographical indications (PGIs) and protected designations of origin (PDO).
With 18,170 tons and 112 millions of Euros, PR represents the first sheep cheese exported in volume and value from the European Union (CLAL 2020). PR is the main cheese produced from the Italian dairy sheep milk, with half of the national sheep milk production processed into this type of cheese (Pulina et al. 2018). The 95% of Italian PR is produced in Sardinia from the Sarda local breed of sheep, which counts 2,317,000 lactating ewes (31st December 2019, Anagrafe Nazionale Zootecnica 2020).
In Sardinia, most of the dairy sheep farms are located in the plains or in the hills, at altitudes below 500 m. The breeding system is semi-extensive, where ewes graze on natural and/or artificial pastures, which represent the main source of nutrients, constituting about 80% of the dry matter (DM) annually ingested by the flock. The remaining 20% of DM consumption is represented by roughages and concentrates. Silagebased feeding systems with limited grazing are not commonly observed except for intensive farms with feedstuff self-production.
Based on the Anagrafe Nazionale Zootecnica (2020), the average size of Sardinian dairy sheep farms is 208.3 ewes and 5.6 rams. Milk production per ewe per lactation ranges from 120 to 150 litres in 120 days in primiparous, from 150 to 230 L in secondiparous and from 180 to 260 L in multiparous, both in 180 days of milking (Pulina et al. 2018).
Cheese yield depends on milk fat and protein contents. Pirisi et al. (1994) have empirically formulated a PR yield equation based on these compounds. Monthly evolutions of milk fat and protein contents are deeply associated to the availability and quality of herbage during the milking season and to the physiological evolution of the lactation curve of sheep. Fatty acid (FA) profile of milk fat, which represents an important feature of the nutritional characteristics of sheep cheese, is also related to pasture quality and energy balance of lactating ewes (Nudda et al. 2014). Increase of milk somatic cell counts (SCCs) enhance the PR yield at 24 h due to higher humidity of the curd (Pirisi et al. 1996(Pirisi et al. , 2000, thus misrepresenting the real cheese yield, that will be lower after long ripening time. Higher bacterial counts could interfere on cheese ripening because sporophytes have higher changes to survive the milk-heating process (Fuka et al. 2013).
PR Consortium of Protection has recorded, for 5 years, the monthly amount of produced cheese, the quantity of milk processed into cheese by the 36 dairies belonging to the P.D.O. area; for 28 of these dairies, bulk milk composition, including FA profile whose nutritional value is transferred to cheese, were bi-weekly registered at farm level. Recorded data is considered to be representative of the sheep milk produced in Sardinia and constitutes the seasonal fingerprint of its quality and aptitude to transformation in PR cheese.
This article, which analyses the largest data set ever examined on sheep milk quality, was edited to answer to some of the issues of concern for PR producers. In particular, monthly variation of milk composition is of interest because it helps to program the output of milk processing plans; moreover, the content and evolution of FA profile are important to determine the nutritional and technological prosperities of PR. Thus, aims of this article were: (a) to describe and explain the monthly evolution of Sardinian sheep milk quality and the magnitude of its variability over 5 years; (b) to evaluate the influence of month and year on PR cheese yield; (c) and to calculate the PR yield prediction equations based on fat and protein concentration of milk, considering the variability in cheese making processes adopted by the different dairies into an unique a dimensional variable.

Materials and methods
A total of 1096 records was collected from 36 dairies in a 5-year period from 2015 to 2019, for 7 months, from January to July. For each month and each dairy, total milk collected, the quota processed into PR and the cheese obtained at 24 h after curding, were considered.
The data on sheep milk composition for 28 out of 36 dairies producing PR were also elaborated. The milk composition data [fat, protein, casein, and lactose percentage, fatty acids (individual FA: C18:1t11, C18:3n3 and CLAc9,t11; groups of FA: SFA, UFA, MUFA and PUFA, standing for: saturated, unsaturated, monounsaturated and polyunsaturated FA, respectively) content, SCC, total bacterial counts (TBC), urea, pH and NaCl] were provided by the laboratory of ARAS (Sardinian Breeders Regional Association). Briefly, bulk milk samples were collected every two weeks from January to July and were stored at 1-4 C to be analysed. Fat, protein, casein and lactose contents of each farm bulk milk sample were determined using a FIL-IDF procedure FIL-IDF, International Dairy Federation (2000), with a Milkoscan-6000 (Foss Electric, Hillerød, Denmark) calibrated for sheep milk; the same equipment was used to obtain data of urea and NaCl concentration. SCC and TBC were determined using a Fossomatic360 instrument (Foss Electric, Hillerød, Denmark); milk pH was also measured.
The mid-infrared (MIR) analysis of milk, carried out in the same ARAS laboratory, with a MilkoScanFT6000 equipment (Foss Electric, Hillerød, Denmark), provided the estimation of the contents of three individual FA (C18:1t11, C18:3n3 and CLAc9,t11), and four groups of FA (saturated, SFA; unsaturated, UFA; monounsaturated, MUFA and polyunsaturated, PUFA). The milk content of these three individual FA and four classes of FA were obtained using previous prediction equations developed by Caredda et al. (2016). Concentrations of individual FA were expressed as grams per 100 g of total fatty acids methyl esters (FAME), whereas the concentrations of classes of FA were originally expressed as g per 100 g of milk, and then transformed (for the present work) in grams per 100 g of fat.
The composition of milk samples was summarised by the ARA laboratory and organised per month, for each dairy, starting from a dataset of 9,537,952 records.
In the preliminary data editing, inconsistent records (4 rotational records) were discarded. The monthly evolution of milk collected by all dairies, expressed as % of the annual overall, was regressed on the seven months under study to ascertain the dynamic of milk availability over the production season of PR.
A combined parameter, named yield, given by the ratio between the produced PR 24 h_cheese (kg), and litres of milk used, was analysed by using two different mathematical models. The first was aimed to compare the mean yield differences across years, months and their interaction. The following two factorial ANOVA model was used: where m was the overall mean; yr i was the i-th year (from 2015 to 2019); mth j was the j-th month (from January to July); yr i Ã mth j was the interaction factor between the i-th year and the j-th month and e ij was the random residual. The Tukey test was used for pairwise comparisons and the p-value significance threshold was set at a 0.05 level.
The second model was aimed to test yield for possible differences in the 28 dairies, considering the fat and the protein percentage equal in all dairies. These percentages are not, in fact, controlled by the single dairy. They depend on the quality of the milk that is conferred and processed. In consequence, possible yield differences among dairies can be attributed to the farm management. With this aim, the analysis of covariance (ANCOVA) was used to compare the dependent variable (the yield) with the fixed effect of a categorical variable (the 28 diaries) considering (or correcting for) the variability of other continuous variables, called covariates, (the fat and the protein percentage). The model used was: where, y i is the yield in one of the 28 diaries; b 0 is the overall intercept; d is the categorical variable representing the 28 dairies (i ¼ 1, … ,28 levels); fat (f) and protein (p) percentages are the two covariates and e is the random residual. The model produces 28 regression equations, which differ from each other only for the parameter b 1i (one different value for each dairy), being b 2 and b 3 the same for the 28 dairies. Finally, a one-factor ANOVA model was developed to test for possible monthly differences of chemical (fat, protein, casein, lactose, pH, urea, NaCl, SFA, UFA, MUFA, PUFA, C18:1t11, CLAc9,t11 and C18:3n3) and cellular composition (SCC and bacterial load) of milk. Also in this case, the Tukey test was used for pairwise comparisons and significance was declared at p < .05. Descriptive statistics (mean ± SD) were calculated and statistical analyses were performed using the SAS software (ver. 9.4, SAS Institute Inc., Cary, NC).

Milk collected and transformed in PR
Milk litres transformed in PR showed a fluctuating trend across the years, being the maximum close to 190 million litres (2016) and the minimum around 150 million litres (2019). The studied dairies transformed almost 70% of the total milk collected into PR, with this value ranging from 63.4% in 2019 to 77.7% in 2015. The remaining milk was usually transformed into Pecorino Sardo P.D.O. or other categories of soft and hard sheep cheeses. Monthly evolution of milk collected by dairies, expressed as percent of the annual overall, fits a quadratic model (Figure 1), with the maximum in April and minimum in July, which reflects the seasonality of the dairy sheep farming system mainly based on the pasture growth seasonality.
As showed in Figure 2, Mediterranean pastures reach medium productivity in autumn, low productivity in winter, high productivity in spring and scarce productivity in summer. Molle et al. (2008) described the pattern of the feeding strategy most adopted in Sardinia to manage pastures and integrate dairy sheep diets to maximise productivity and animal welfare.
Season has a strong effect both on the biomass availability and quality.

Cheese yield
The average yield of PR was 0.174 kg of 24 h_cheese per litre of milk, with a standard deviation of 0.017. The general output of the model indicated that the yield was different for year and month, with a significant interaction between these two factors. The yield in 2018 (0.176 kg of 24 h_cheese per litre of milk) was significantly higher than the one in 2016 (0.173; p < .05), whereas the other years did not differ from each other (0.175, 0.175 and 0.174 kg of 24 h_cheese per litre of milk for 2015, 2017 and 2019, respectively; p > .05). Regarding the months of production ( Figure  3), April, June and July showed significantly different yields from the other months. In particular, the lowest yield (i.e. more litres of milk needed per kg of cheese) was recorded in April, whereas the highest in June and July (i.e. less litres of milk per kg of cheese).
The PR yield is a function of the general percentage of bulk milk fat and protein processed in the 28 dairies (see below) and can be expressed by the following linear regression model with both categorical and continuous variables: where, b 0i d is different for each of the i-dairies (i ¼ 1, … .,28) with mean value 0.0757 and SD 0.0059. This relationship highlighted a higher value of protein than fat in cheese yielding, with a ratio between them equal to 1.16. This value was slightly lower than the one reported by Pirisi et al. (1994)

Percent of milk collected
Months (1 = January; 7 = July)) Figure 1. The monthly evolution of milk collected in 5 years by the 36 dairies member of the Pecorino Romano Consortium, expressed as percent of the annual overall (intercept not significantly different from zero).

Milk quality
The descriptive statistics of milk physicochemical traits (fat, protein, casein, lactose, SCC, TBC, pH, NaCl and C18:1t11, C18:3n3, CLAc9,t11, SFA, UFA, MUFA and PUFA) of the 28 dairies producing PR are summarised in Table 1. Average fat and protein contents were around 6.5 and 5.6%, respectively, whereas lactose content was about 4.8% and SCC was higher than 1 million per ml. As far as the FAs composition is concerned, SFA represented more than 68% of the total FA, whereas MUFA and PUFA accounted for about the 25 and 6%, respectively.
Monthly evolution of fat (Figure 4) was associated to the quantity of the milk delivered, confirming the concentration effect common to all dairy species (Pulina et al. 2003). Lowest fat content was observed in April (5.85%), whereas the highest one in July (7.52%). The evolution of C18:1t11 (vaccenic acid), C18:3n3 (alpha-linolenic acid) and CLAc9,t11 (rumenic acid) was coherent with the maturation trend of Mediterranean pasture grass, being these FA clearly related to C18:3n3 content in grass that is maximum in April and deceases rapidly until July (Dewhurstet al. 2006;Mele 2009).
On the other hand, the evolution of FAs classes showed in Figure 5 highlights the interaction between pasture grass quality and energetic metabolism of ewes. PUFAs, as reported above, are associated with the C18:3n3 content in the fresh herbages (Cabiddu et al. 2005;Mel 'uchov a et al. 2008). On the other hand, higher SFAs in April and their decrease in the next months, followed by the specular trend of UFA and in particular of MUFA, can be associated to the losses of weight occurring in lactating sheep in the last part of lactation (Gallagher et al. 1966), caused by the worsening of grass quality (Molle et al. 2008) and the stress caused by the high temperature (Seijan et al. 2017). Indeed, in these conditions, an increase in body reserve mobilisation can occur, with consequences on the milk FA composition. In particular, an increase in C18:1cis-9 concentration (MUFA are composed in large part by this FA) and a decrease of de novo synthetised FA (representing a great proportion of the SFA), have been related to the body fat mobilisation (Moallem et al. 2000). The association of body fat metabolisation with the feed restriction was pointed out in a recent work on Sarda dairy sheep (Correddu et al. 2021).
Proteins content evolution ( Figure 6) was quite different from that of fat in the second and third month of lactation: the relative higher value found in February and March, compared to January, April and    May, reflected the energy balance of the ewes, which is mainly positive; on the contrary, the rapid increase in the last two months reflects the concentration effect, as seen for fat content. Monthly evolution of casein percentage overlays the one of protein; milk urea was higher in winter and early spring, according to the soluble nitrogen and protein contents of grass (Pulina et al. 2006;Molle et al. 2008). The milk urea curve usually follows the availability and quality of the pastures, with high values during spring months (Molle et al. 2008). In the last decade, a strong decrease in Sardinia's milk urea values in pastures has also been achieved, leading to a major decrease in the nitrogen excretion in the urine of the ewes reared in Sardinia (Nudda et al. 2020); this was mainly due to the use of urea as a nutritional protein feeding index (Cannas et al. 1998;Giovanetti et al. 2019) . Milk urea is also strongly and negatively correlated to SCC in Sarda sheep milk (Nudda et al. 2020) and is reduced by dietary concentration of tannins (Molle et al. 2008;Correddu et al. 2019). Milk pH monthly evolution shows a significant narrowed wavy trend (Figure 6), remaining in the range of normality for all lactation period (Park et al. 2007).
The milk monthly trends of lactose, NaCl and SCC (Figure 7), showed the classical involution of the mammary gland: decrease of principal osmotic component at lobular level, twined with the increase of endogenous cellular component, partially caused by re-modulation of the secretive tissue (Nguyen and Neville 1998). Finally, bacterial count remained constantly low across the seasons, except for the last month (July) probably because of the combined effects of high temperatures and the less strict sanitary routine protocol at the end of the productive season.
The quality fingerprint of the Sardinian sheep milk represents the complex combination of: (a) lambing seasonality, with the majority of the multiparous ewes delivering in late autumn and the primiparous and the remainders multiparous (usually the 10-15% of flock) in late winter (Sitzia et al. 2015); (b) the lactation curve and the pattern of milk delivery, which in good feeding conditions, follow the Wood model (Cappio-Borlino et al. 1995); (c) the trends of pasture grass quality and quantity, that in Mediterranean areas reach the highest values in April (Molle et al. 2008) quality trends that, in Mediterranean pastures usually reach their lower quality in late-spring, due to the increase of maturity stage of the plants (high NDF and low CP). To compensate pasture quality variability and scarcity, sheep farms usually integrate animal diets with grain and mixes concentrates from both local and international markets. It has to be noticed that milk FA composition is also related to the herbage quality and to the intake of glucogenic and lipogenic carbohydrates; elevated intake of FA precursors from high digestible forages is associated with increases in de novo FA (Woolpert et al. 2017). Thus, an increase in de novo FA is expected in winter and early-spring with high herbage availability and intake, whereas a decrease is expected late-spring and summer because the lower pasture fibre quality.

Conclusions
This work is the biggest picture population-based of the Sarda sheep milk quality, its monthly evolution and of its yield rate for the PR production. The yield in PR depends, for more than half of its variability, on milk fat and protein contents, but it is also influenced, by 43% of the variability, by the technologies adopted by the dairies.
The monthly trends of Sardinian sheep milk quality and PR yield reflected the combination among animal factors (lactation curve), environmental constraints (pasture grass and climate seasonal evolutions) and technical choices made by the shepherds (dynamics of seasonal lambing and feeding strategies). The data here reported, because based on large majority of Sardinian sheep population, represent the reference values of Sarda dairy sheep performances regarding milk quality and its ability to be transformed into PR cheese. These data could be a relevant basis to set up future grids of milk payments based on quality standards at regional and dairy levels, to formulate administrative policies on dairy sector aimed to improve milk quality of Sardinian sheep milk, and moreover for the milk destined to PR production and diversification. Data showed in the present study will be also relevant to support basic and applicative research on sheep milk quality of Sarda dairy sheep.
A couple of take-home messages can be drawn by the results obtained in this work. First, pasture-based system of Sardinian dairy sheep industry leads to a strong seasonality not only in milk yield, but in cheese yield rate too causing scheduling problems for the cheese factories. Different lambing programs, which provide for the continuous production of milk throughout the year by splitting the flock into two breed-season groups, and feeding strategies that are able to absorb the quantitative and qualitative variability of pasture, are mature technologies that could  Figure 7. Monthly evolution (means ± SE of 5 years) of milk traits: lactose (%), NaCl (mg/100 mL), Somatic cell count (Â1000/mL; geometric mean) and bacterial load (Â1000/mL; geom. mean); different letters indicate significant differences (p < .05). be easily implemented at farm level for a majority of Mediterranean dairy sheep enterprises. Second, the high dependence level of cheese yield on technologies adopted at the specific cheese factory, demonstrates that the final transformation value of milk in PR cheese depends less than an half of its variability on milk quality, suggesting that the milk market have to be thought as unique at regional level and opening wide rooms to improve the competitivity among transformers.

Ethical approval
Because the work did not involve animals, it was not necessary to obtain permission from the ethics committee. All the used data were provided by the Pecorino Romano Consortium.