Introduction
The buffaloes in Serbia belong to the Mediterranean buffalo breed descended from the Arni river buffalo. They came to Europe and thus to Serbia in the 6th and 7th century (Moioli et al., 2001; Zhang et al., 2020). The largest number of buffaloes in Serbia is in southwestern Serbia in the Raška region, in the municipalities of Novi Pazar, Sjenica and Tutin.
The buffaloes in Serbia were not given special consideration during selection, they were treated as the most primitive group of cattle, which resulted in lower buffalo milk production in Serbia compared to buffalo milk of Mediterranean breed in other countries (Joksimović, 1968; Perišić et al., 2015). However, governmental programs to protect genetic resources supported buffaloes and milk production.
Buffalo milk is characterized by high levels of fat, lactose, protein, and ash contents. Due to its composition, buffalo milk is used in Serbia to produce a variety of dairy products such as white brined cheese and of fat-based dairy products like butter and kajmak, a traditional Serbian product (Mitić et al., 1987; Becskei et al., 2020). Cheeses are made both from raw milk and heat-treated (95 °C/10 min) milk remaining after kajmak production. These products are very important in the abovementioned region but to the best of our knowledge have not been investigated in detail, as well as buffalo milk used for them.
The casein fractions play a fundamental role in the production of cheese through their differences in primary structure among the ruminant species and breeds of dairy interest.
Casein fraction of buffalo milk, as in other ruminant species, is composed of αs1-casein (αs1-CN), β-casein (β-CN), αs2-casein (αs2-CN) and k-casein (k-CN), encoded by four major genes specifically expressed in mammary gland epithelial cells during lactation: CSN1S1, CSN2, CSN1S2, and CSN3, respectively (Balteanu et al., 2013). The study of milk protein variants has proven to be a useful tool for associating genetics to productive traits, as has already happened with dairy cattle.
The αs1-CN of buffalo milk consists of a single polypeptide chain with 199 residues and shows high homology to variant B of the αs1-CN of cow milk (D’Ambrosio et al., 2008). Genetic variation of casein genes in buffaloes ( Bubalus bubalis) has been poorly studied so far (Balteanu et al., 2013; Zicarelli et al., 2020). Only two genetic variants have been designated as buffalo αs1-CN A and αs1-CN B, which have a single silent amino acid substitution Leu178(A)→Ser178(B) with allele frequencies of 0.31 and 0.69, respectively (Chianese et al., 2009; El Salam and El-Shibiny, 2011). Balteanu et al. (2013) were the first to find and describe the exon skipping phenomenon and a novel gene variant with an allele frequency of 0.18 in the CSN1S1 locus of Romanian buffalo.
One of the main processes during cheese production is the two-stage rennet coagulation, which depends on numerous factors and significantly contributes to gel properties as a base for the final product structure. It is known that the primary phase presents enzymatic hydrolysis of k-CN, while the secondary phase is the aggregation of destabilized casein micelles in the presence of Ca2+ contributing to gel formation. Factors influencing the coagulation parameters and gel properties include ionic strength, casein concentration, milk heat treatment, coagulation temperature, milk pH, amount of rennet used and calcium concentration (Dalgleish, 1983; Van Hooydonk et al., 1987; Nájera et al., 2003).
The dynamic oscillatory rheological measurements have been widely used to monitor the kinetics of milk rennet coagulation and the properties of the obtained gel (Zoon et al., 1988; Lucey, 2002; Lucey et al., 2003). The dynamic measurements show that during gel formation and aging, the modulus of elasticity (G') and modulus of viscosity (G") increase with time, while the ratio of modulus G" to modulus G', tan δ, decreases very quickly when the moduli start to increase, and it usually becomes constant during the coagulation process as a result of the transition from a fluid into a gel system (Zoon et al., 1988). The rate of the modulus of elasticity increase depends on different factors, such as the casein content, rennet concentration, milk pH, temperature, etc. According to Dimassi et al. (2005), casein micelle aggregation, which occurs in the nonenzymatic phase, can be characterized by rennet coagulation time (RCT), aggregation rate (AR) and gel firmness (GF) obtained based on dynamic oscillatory rheological measurements.
The genetic variant distribution of the caseins and major whey proteins, β-lactoglobulin (β-LG) and α-lactalbumin (α-LA), contribute to the difference in the milk coagulation properties (Čítek et al., 2020; Sabour et al. 1996; Bonfatti et al., 2010, 2012; Zicarelli et al., 2021., La Gatta et al., 2021). Thus, the αs1-CN B alleles in milk contribute to better curd yield efficiency. Furthermore, in the cattle breed, milk with the αs1-CN B/k-CN B alleles show the best coagulation efficiency (Caroli et al., 2009). McLean (1984) indicates a relationship between a higher concentration of αs1-CN and a concomitant decrease in κ-casein in cow milk, which affects the process of rennet coagulation.
Genetic variants in the Italian buffalo milk proteins and the role of αs1-CN and β-CN genes, in the coagulation properties affecting the RCT and the setting time of the curd have been described by different authors (Cecchinato et al., 2012; Dettori et al., 2009; Chianese et al., 2009; Zicarelli et al., 2021). Bonfatti et al. (2012) confirmed the existence of CSN1S1 and CSN3 genetic variants in Mediterranean water buffalo and assumed that the B allele at CSN1S1 was associated with increased RCT and curd firming time and with weaker curds compared with allele A. However, Bonfatti et al. (2013a) pointed out that the effects of genotype on rennet coagulation properties were not entirely due to existing differences in milk protein composition related to variants in the CSN1S1 and CSN3 genes, but primarily to total casein and k-CN content.
Milk composition obtained from buffaloes raised in Serbia, particularly casein gene variants, has not been studied much. Hence, the aim of this study was to analyze the CSN1S1 472G>C variant and its potential influence on rennet coagulation properties that are crucial for cheese production. Rennet coagulation properties were studied on raw and high heat-treated milk, usually used for cheese production in Serbia. Data obtained in this study will contribute to the limited knowledge on buffalo milk properties originating from Serbia and could help in profiling the selection process direction to improve its technological properties and get more economic processing.
Materials and methods
The Ethical Permission for the use of animals for the research implementation was issued by the Ministry of Agriculture, Forestry and Water Management of the Republic of Serbia, number 323-07-10974/2022-05. All rules related to animal welfare were followed.
Genotyping
Genetic analysis was performed on blood samples from a total of 30 female buffaloes of the Mediterranean buffalo breed raised in Novi Pazar and Tutin (Serbia). DNA was extracted by the phenol-chloroform extraction method.
CSN1S1 (Gene ID: 102396531) 472G>C was genotyped by Sanger sequencing method in house, on ABI PRISM 3100 Genetic Analyzer (Applied Biosystems, USA), using BigDye® Terminator v3.1 Cycle Sequencing Kit protocol. PCR reaction was performed using the following primer sequences: forward primer I5CZS1-F: 5’ -ACT TAG CAA GGA GAT AAT GCA AGA A-3’ and reverse primer E7BCZS1-R: 5’ - CTC AGT TGA TTC ACT CCC AAC ATC-3’ were used for amplifying genomic region from intron 5 to exon 7 of CSN1S1. PCR conditions were: 94 °C for 3’, than 35 cycles of 94 °C for 1’, 58 °C for 1’ , 72 °C for 1’ and final step of 72 °C for 7’ (Balteanu et. al., 2013). Genomic DNA (100 ng) was amplified in a total volume of 20 μL reaction mixture containing 100 ng of each primer, 200 μmol/L deoxynucleoside triphosphates; 1 U of DreamTaq polymerase in 10× DreamTaq buffer (Thermo Fisher Scientific Inc.).
Detection of genotypes was conducted using Sequencing Analysis Software V4.0 (Applied Biosystems™, USA). Genotypes were coded and included into statistical files for analysis.
Sampling of buffalo milk
After completion of genotype analysis, raw milk samples were collected weekly for three consecutive weeks in February 2023 from buffaloes of each genotype in the same lactation period. The milk samples were immediately cooled and analyzed.
Buffalo milk composition analysis
The composition of raw buffalo milk samples was determined by the following methods: total solids by the standard drying method at 102±2 °C (IDF, 1987); fat content according to the Gerber method (IDF, 1981); nitrogen content by the Kjeldahl method (AOAC, 1990), while the protein content was calculated as the nitrogen content multiplied by 6.38. Milk pH was measured using a digital pH-meter (Mettler-Toledo GmbH, Switzerland).
Experimental design
The experiment design presented in Table 1 was used to examine the influence of genotypes of CSN1S1 472G>C on rennet coagulation and the gel properties including the fresh curd yield.
Three pooled milk samples were prepared, each representing a specific genotype (GG, GC, and CC). Each pooled sample, totaling 1 liter, had an identical protein content of 4.9 %.
The milk samples were defatted by centrifugation at 2,000 g for 30 min at a temperature of 40 °C (Eppendorf centrifuge model 5430; Eppendorf AG, Hamburg, Germany) and then the fat layer was removed. Defatted milk samples were divided into two batches: a batch of raw milk (without heat treatment) and one of high heat-treated milk (at 95 °C/10 min). The heat-treated milk was immediately cooled and stored in the refrigerator.
Before addition of rennet the milk was incubated at 40 °C for 60 min in order to reverse the dissociation of casein and calcium from casein micelles caused by cold storage (Hallen et al., 2007) and then was cooled to the appropriate rennet coagulation temperature (31 °C).
Table 1. Buffalo milk samples and rennet coagulation parameters used in the experimental design
Milk coagulation was initiated by the addition of rennet powder (Caglificio Clerici spa, Cadorago, Italy, 94 % of chymosin) in the quantity of 0.003 % (w/v). Rennet dilution in distilled water (concentration 1.2 %) was prepared 15 min before the addition into the milk sample.
Rennet coagulation parameters were achieved by following: the adjustment of pH by adding 10 % lactic acid (Merck, Germany); the level of added CaCl2 by using 20 % solution of CaCl2 (Table 1).
Rennet coagulation parameters were achieved by using 10 % lactic acid (Merck, Germany) for pH adjustment and 20 % solution of CaCl2. The pH of high heat-treated buffalo milk was adjusted to 6.3 to enable rennet coagulation (Radovanovic et al., 2021).
Each milk sample prepared was divided in two loads for rheological measurements (40 mL) and a fresh curd yield determination (20 mL), respectively.
Sodium dodecyl sulfate polyacrylamide gel electrophoresis
Raw milk (GG1, GC1, CC1) and heat-treated defatted buffalo milk (GG2, GC2, CC2) were analyzed by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE) under reducing (SDS -R) and non-reducing conditions (SDS-NR) according to the method of Laemmli (1970) to study the influence of the CSN1S1 genotypes and heat treatment of milk on the protein fractions distribution.
Resolving and stacking gels were prepared according to the method of Anema and Klostermeyer (1997). The resolving gel contained 15 % acrylamide (2.6 % BIS) and 0.1 % SDS in Tris/HCl buffer at pH 8.85, while the stacking gel contained 4 % acrylamide and 0.1 % SDS in Tris/HCl buffer at pH 6.8.
Sample buffers (for SDS-R and SDS-NR) at pH 6.8 and the running buffer at pH 8.30 were prepared according to Miloradovic et al. (2015). Buffalo milk samples were skimmed by centrifugation at 600 g at 20 °C for 15 min respectively (centrifuge model 5430; Eppendorf AG, Hamburg, Germany). Centrifugal supernatants were diluted with the sample buffer to achieve a protein concentration of 3 μg/μL. The prepared samples were boiled for 5 min. The total sample quantity of 5 μL was loaded onto the gel.
SDS-PAGE was performed using a vertical twin-plate electrophoresis unit (TV200YK, Consort, Belgium). Electrophoresis was done under the following conditions: I const. = 80 mA, U max. 300V. Gel was run in the pH 8.30 running buffer for 75 minutes. Gel staining was performed in 0.23 % Coomassie Blue R-250 dye solution and destaining in a solution containing 8 % (v/v) acetic acid and 18 % (v/v) methanol (Miloradovic et al., 2015). Densitometric analyses of the SDS electrophoretograms were performed by using densitometric software ImageJ (National Institute of Health, Bethesda MD, USA).
Rheological properties of buffalo milk rennet coagulation
Rheological measurements of milk rennet coagulation and gel properties were performed by dynamic small amplitude oscillatory rheology (SAOR) using a Kinexus Pro rheometer (Malvern Instruments Ltd, Worcestershire, UK) equipped with a vane geometry tool according to Radovanovic et al. (2021).
Rheological measurement was started 4 min after the rennet addition into buffalo milk. The coagulation process was monitored at a constant temperature by measuring the modulus of elasticity, G' and the modulus of viscosity, G'', as well as the related phase angle δ= arc tan (G''/G'). The measurements were done over 60 min, at a fixed frequency of 1 Hz and 0.01 % shear strain. To examine gel firmness (GF), a subsequent frequency sweep test was performed in the range of 0.1-20 Hz at 0.01 % shear strain at the coagulation temperature. GF represents the value of G' at a frequency of 1.5 Hz. Rennet coagulation time (RCT) was estimated at the point when modulus of elasticity reached the value of G' ≥ 1 Pa (Dimassi et al., 2005). The setting time (ST) was determined at the point when the gel reached the value of G' ≥ 60 Pa. The G' value of 60 Pa for ST was defined by preliminary testing on the coagulation process performed by combining the rheological measurements (data not shown) and visual control of gel strength. Additionally, G’ at 60th min of rheological measurement was recorded as another indicator of the gel firmness (G'60).
Aggregation rate (AR), in mPa s-1, was obtained from the slope of the tangent at the point corresponding to the maximum of the first derivative of G'(t), according to Steffl et al. (1999). Based on the data obtained by rheological measurements, coagulation parameters (RCT, AR, GF) were calculated using Origin Pro 8.0 (OriginLab Corporation, Northampton, USA). Two replicate rheological measurements were carried out for each sample in each batch.
Fresh curd yield (Yf)
Fresh curd yield (Yf) was determined according to the method of Hallen et al. (2010) with minimal modifications. Defatted milk samples (10 mL) were placed in test tubes. Rennet solution (25 µL) was added to each sample, after which they were mixed and kept undisturbed in a water bath at coagulation temperature 31°C. After reaching ST (60Pa), the curd was vertically cut in four equally sized sections, using a four-edged knife. After another 30 min incubation the tubes were removed from the water bath. Pressing of the curd was simulated by centrifugation at 1.260 g for 20 min at 31 °C (Eppendorf centrifuge 5430, Eppendorf AG, Hamburg, Germany). The expelled whey was separated, weighed, and expressed as grams of whey per 100 g of milk sample. The fresh curd yield (Yf) was calculated as the weight difference between the initial milk sample and the expelled whey and expressed as grams of curd per 100 g milk. The test was performed in duplicate.
Statistical analysis
The one-way analysis of variance was calculated to test the influence of CSN1S1 472G>C variant on rennet coagulation parameters and gel properties of buffalo milk obtained by using STATISTICA 6.0 (StatSoft, USA) data analysis software. Mean comparisons of the parameters were performed using the LSD-test, with a significance level of 0.05.
Results and discussion
Genotyping
Genotyping results demonstrate the existence of three CSN1S1 472G > C genotypes: GG, GC, and CC in investigated buffalo population. The frequency of the genotypes is presented in Table 2.
In the tested sample of 30 buffaloes, the GG genotype was the most prevalent, with a frequency of 53.33 %. The genotype CC had the lowest frequency, at 10 % (Table 2).
Table 2. Genotype frequencies of CSN1S1 472G > C gene variant (n = 30)
Buffalo milk composition
The average protein content was 4.9 %, fat content was 9.40 %, and dry matter content was 20.10 %. The composition of buffalo milk samples fell within the range of variations in the chemical composition of buffalo milk reported in the literature (Tufarelli et al., 2008; El Salam and El-Shibini, 2011; Braun and Preuss, 2008; Radovanovic et al., 2021) for Italy, Serbia, Germany, Turkey, and other countries outside the European continent.
SDS-PAGE
By densitometric analysis under SDS-R conditions it was established that raw milk samples of all three genotypes have a similar percentage of total casein which amounts to about 72-75 % and the whey protein fractions (β-LG + α-LA) about 23-24 % (Table 3) which is in accordance with the literature data (Islam et al., 2014). In addition, high weight fractions, including immunoglobulins, IgM, IgA, and IgG (El Salam and El-Shibini, 2011) are visible in these raw milk samples being the most abundant in the CC genotype (Figure 1).
Densitometric analysis of casein fractions in the raw milk samples showed that the percentage of total αs-CN in total casein fractions was about 53.3 % (GG); 54.3 % (GC) and 52.5 % (CC). At the same time, the proportion of k-CN in total casein was 19.6 %, 15.11 % and 16.5 % for GG, GC and CC, respectively. The percentage of k-CN in total casein was slightly increased compared to literature data (Bonfatti et al., 2013b; Cosenza et al., 2011), especially for the GG genotype.
The high heat treatment of milk (90 °C/10 min) applied in the production of kajmak leads to strong changes in milk proteins, such as denaturation of whey proteins and formation of coaggregates of milk proteins (Corredig and Dalgleish, 1999; Anema and McKenna, 1996; Anema and Li, 2003) which has been proven in various examinations of the kajmak formation (Radovanovic et al., 2020; Radovanovic et al., 2022).
Based on the electrophoretogram (Figure 1) all three samples (genotypes) of heat-treated milk in SDS-NR conditions have a similar profile of protein fractions. In these samples, the β-LG fraction is missing, which means that β-LG participated in the formation and is completely bounded by disulfide bonds in high molecular weight complexes, coaggregates with some of them positioned at the entrance of the stacking gel (Figure 1).
At the same time, α-LA, also participates in these complexes, but is not completely bounded, so a certain amount of free α-LA remained. The incomplete aggregation of α-LA in all samples is probably the result of an insufficient duration of the heat treatment of milk necessary for the complete denaturation of α-LA (Corredig and Dalgleish, 1996; Anema, 2001; Anema and Li, 2003). In addition, the results of densitometric analysis of the samples under SDS-NR conditions showed that despite the formation of coaggregates it could be observed that a significant number of unbound caseins remained in all samples.
Figure 1. SDS electrophoretograms of row (GC1, CC1, GG1) and heat-treated (GC2, CC2, GG2) buffalo milk under reducing (SDS-R) and non-reducing (SDS-NR) condition
Table 3. Ratios of protein fractions (%) in raw buffalo milk
GG1, GC1 and CC1-raw buffalo milk of GG, GC and CC genotype respectively
Rheological properties of buffalo milk rennet coagulation
Based on the data from rheological measurements it can be concluded that there are significant statistical differences between the tested parameters of CSN1S1 472G > C genotypes of buffalo milk, which are expressed in both raw and heat-treated milk (Tables 4 and 5).
Table 4. Average values of rennet coagulation parameters of raw buffalo milk
*Values with different letters within the same column are significantly different (p<0.05)
GG1, GC1 and CC1-raw buffalo milk of GG, GC and CC genotype respectively
Table 5. Average values of rennet coagulation parameters of heat-treated buffalo milk
*Values with different letters within the same column are significantly different (p<0.05)
GG2, GC2 and CC2 - heat-treated (95 °C/10 min) buffalo milk of genotype GG, GC and CC respectively
When comparing the rennet coagulation parameters of raw milk, GG genotype had a very similar RCT to GC genotype, although GG has a higher content of k-CN, which improves the coagulation properties of milk and can be considered as a factor for accelerating the enzymatic phase of rennet coagulation (El Salam and El-Shibini, 2011). This was supported by our results with an exceptionally high value of the coagulation rate, AR (0.1232 Pa/s). On the other hand, compared with GC and GG genotypes, CC genotype had slightly less total αs-CN and a shorter RCT was expected (Bonfatti et al., 2013a), but a longer RCT was observed. With the electrophoretic technique mentioned above, it was not possible to determine the involvement of αs1-CN and αs2-CN because the peaks overlapped, and it is possible that this genetic variant (CC) had a higher involvement of αs1-casein, causing the prolongation of RCT.
Results (Table 4, Figure 2a) showed that the genotype GG was also characterized by high values for G'60 (212 Pa) and for GF (233 Pa), which were statistically significantly different (p<0.05) compared to the parameters of the variants GC and CC. The reason for the lower values for G' at the end of the measurement and for GF of the genotype CC could be the higher proportion of β-CN. According to St-Gelais and Hache (2005), milk with a relatively high proportion of β-CN is characterized by poorer coagulation properties.
Considering heat-treated milk, genotype CC was characterized by high values of the rennet coagulation parameters. Recorded data for G' at the end of the measurement (about 178 Pa) and especially for GF (about 198 Pa), were statistically significantly different from the parameters of the genotypes GC and GG (p<0.05). The evolution of the G’ during coagulation of the heat-treated milk was graphically represented in Figure 2b where it could be observed that the high heat treatment had the least influence on the dynamics of the G' of the milk genotype CC, while the strongest influence was recorded for the GC genotype.
GG1, GC1 and CC1-raw buffalo milk and GG2, GC2 and CC2 – heat-treated (95°C/10 min) buffalo milk of genotype GG, GC and CC respectively
Figure 2. Dynamics of the modulus of elasticity (G') during rennet coagulation of raw (a) and heat-treated buffalo milk (b)
Raw and heat-treated buffalo milk were characterized by similar renneting properties, regarding RCT and ST values. However, the other parameters were completely different (Tables 4 and 5). Highly heat-treated buffalo milk showed a much slower aggregation rate compared to raw milk, with lower values for GF and G'60, as well. An exception was the genotype CC which showed even higher values for these parameters compared to raw milk (Table 4, 5). Nevertheless, the results clearly suggest that highly heat-treated buffalo milk of these genotypes could be suitable for rennet coagulation if the pH of the milk is adjusted to 6.3.
These results are consistent with the conclusions of Radovanovic et al. (2021), but in the current study, we obtained lower values for ST and rather higher values for RCT, G'60, GF, and AR for both raw and heat-treated milk.
Bearing in mind that in Serbia it is common to use highly thermally processed cow milk which remains after the production of kajmak, we thus confirmed that buffalo milk could be used in the same way. These results could contribute to valorizing the usage of buffalo milk, which has economic importance for many people living in the rural regions in South-Western Serbia, where these animals are predominantly raised.
Fresh curd yield
The fresh curd yield (Yf) obtained from fresh and heat-treated buffalo milk is shown in Tables 6 and 7. Values of the amount of whey separated by centrifugation method, which simulates the pressing of the cheese curd, are also expressed in percentage. The highest yield was obtained with genotype CC in both, fresh and heat-treated milk. In general, for all genotypes, heat-treated milk had a higher curd yield than raw milk.
A particularly high Yf (about 72 %) was found in heat-treated milk of the CC genotype, probably due to the pronounced hydrophilicity of the milk protein coaggregates formed, the presence of which is confirmed by the results of SDS analysis under non-reducing conditions (Figure 1).
There are different methods and formulas for calculating the real and estimated curd yield, but they have not given satisfactory results so far (Zicarelli et al., 2020; Zicarelli et al., 2021). We have used the method of determining yield by measuring the whey separated from the fresh curd. This method can be used to compare the curd yield of different milk samples but does not correspond to the actual yield of cheese i.e. real curd yield, RCY (Zicarelli et al., 2020).
Table 6. Fresh curd yield and separated whey obtained from raw buffalo milk
*Values with different letters within the same column are significantly different (p<0.05)
Table 7. Fresh curd yield and separated whey obtained from heat-treated buffalo milk
*Values with different letters within the same column are significantly different (p<0.05)
The Yf values we obtained are generally quite high (Table 6, 7), compared to literature data where the mean fresh curd yield was 23.9, i.e., in the range of 14.5-34.1 % (Hallen et al., 2010; Othmane et al., 2002). The reason is probably that the pressure used was insufficient to pull the whey out from the strong buffalo milk gel. As mentioned earlier, buffalo milk is characterized by a very high protein, fat and mineral content (El-Salam and El-Shibiny, 2011), which contributes to a high firmness of the gel. Thus, it appears that the method commonly used to analyze curd yield and syneresis of cow milk gels (Hallen et al., 2010) is not suitable for determining the gel properties obtained from coagulation of buffalo milk. It is likely that such strength of the buffalo milk rennet gel may be overcome by applying a higher centrifugation force to the gel or a longer centrifugation time. These findings open a space for future research on buffalo milk.
Conclusion
The presence of three CSN1S1 472G > C genotypes (GG, GC, and CC), was determined in buffalo milk samples. The results of this study show that the genotypes affect the variation in rheological properties of rennet coagulation of buffalo milk.
All three genotypes showed high values for GF and G' at the 60th minute of rennet coagulation, which were also pronounced in highly heat-treated milk (95 °C/10 minutes). Moreover, it appears that highly heat-treated milk left over after the production of kajmak could be used for cheese production, although the influence of these genotypes on the choice of cheese type, cheese quality, and especially cheese yield should be further investigated.