Milk somatic cell count-derived traits as new indicators to monitor udder health in dairy buffaloes

Abstract Milk somatic cell count (SCC), an indicator of udder health and milk hygiene, has not been implemented either in milk payment systems or in the current selection index of Italian buffalo so far. As a matter of fact, there is room to improve udder health through genetic and management strategies in this species. Repeated milk SCC test-day (TD) records on the same animal are useful to derive novel phenotypes on lactation basis by using somatic cell score (SCS) and thus to disclose buffaloes with udder issues to be monitored. In this study, sources of variation of SCC-derived traits in buffalo milk were investigated and their effect on milk yield and composition traits was estimated. Mean SCS (SCS150), standard deviation of SCS (SCS_SD150) and severity (SEV150, ratio between the number of TD SCC >200,000 cells/mL and total number of TD available in the first 150 DIM) in the first 150 days in milk were calculated using TD data of 45,312 lactations from 35,623 buffaloes. Sources of variation of such traits were investigated through a linear model. Both SCS150 and SEV150 increased with parity, whereas SCS_SD150 decreased (p < .001). Subsequently, the three traits were separately included in the model as explanatory variables to estimate their effect on milk yield and composition traits. Results showed that milk yield and lactose content were lower in animals with high and variable SCS in the first 150 days of lactation (p < .001). This study opens the debate on the development of an udder health index for the Italian buffalo. HIGHLIGHTS So far, somatic cell count has not been considered in payment systems of buffalo milk in Italy. There is need to improve mastitis resistance and udder health in dairy buffaloes. Traits derived from test-day records can be used for management and breeding purposes in buffaloes.


Introduction
Both milk somatic cell count (SCC) and somatic cell score (SCS) are used worldwide in milk-recording systems as quality and udder health indicators in livestock (Finocchiaro 2018;Franzoi et al. 2020;Costa, Negrini, Campanile et al. 2020). In many countries, regulations on the hygienic and sanitary characteristics of milk apply an SCC threshold to bovine milk intended to direct human consumption. For instance, in Europe, New Zealand, Canada, and Australia the limit is 400,000 cells/mL, and in the USA it reaches 750,000 cells/mL (Norman et al. 2011;Sharma et al. 2011). As regards buffalo milk, the European law only regulates bacterial count, which should not exceed 500,000 colony-forming unit/mL at 30 C when milk is intended to raw-milk cheese manufacture (EU Regulation no. 853/20042004. Italy is one of the major Mediterranean countries rearing buffaloes, together with Turkey and Egypt (Borghese and Moioli 2002). In these countries, buffalo milk is predominantly used to manufacture cheeses and other dairy products, like mozzarella, ghee, and domiati.
In buffaloes under the official milk recording program, SCC is routinely assessed in individual milk.
Despite high SCC is associated with impaired coagulation ability of milk in buffalo (Guccione 2013) like in cattle (Franzoi et al., 2020), still, most of the payment systems of cheese industries do not account for this trait. Test-day (TD) data allow the monitoring of udder health at both individual and farm level across days in milk (DIM); however, the current selection index of Italian Mediterranean buffalo accounts for udder morphology (24%), feet and legs (20%), 270-days milk yield (21%), 270-days fat percentage (15%) and 270days protein percentage (20%), but not for SCC (Associazione Nazionale Allevatori Specie Bufalina 2019).
It is well established that mastitis is often subclinical in buffalo, that is, latent inflammation with no signs, and frequently may become chronic (Sharma and Sindhu 2007;Guccione 2013). As regards causative pathogens, it has been observed that in half of the total mastitis cases there are usually more than one pathogen isolated and involved in the infection and inflammation (International Dairy Federation 2008). Some studies have investigated sources of variation of milk SCC in Anatolian (Şahin et al. 2017, 2019) and Marrah buffalo (Singh and Ludri 2001). Despite his, there is still uncertainty on the level of SCC that may be used to identify mastitis in buffalo and there is a lack of knowledge on udder health monitoring among farmers and milking personnel (International Dairy Federation 2008). In this context, repeated and objective measurements of milk performances, including SCC, may disclose buffaloes with udder issues to be monitored, managed with specific cares (i.e. milked separately) and possibly treated. Recently,  investigated milk SCC in the Italian buffalo population and reported that 28.4% of first-parity buffaloes had lactation average SCC !200,000 cells/mL, and in sixth-parity buffaloes, this percentage was even greater (39%). Thus, it is clear that the use of a single cut-off for SCC could have some limitations, but for routine practice, a single threshold would be much easier to be adopted (Guccione 2013). Some studies on bovine milk have proposed different phenotypes based on TD SCC or SCS to catch more variation of the trait and identify patterns and trends within lactation. For example, Green et al. (2004) reported that the maximum level of SCC observed during lactation and variation of TD SCC were more efficient in predicting clinical mastitis than lactation means SCC. Koeck et al. (2012) suggested that the percentage of TD SCC >200,000 and >500,000 cells/mL were informative in Canadian Holsteins. All the aforementioned studies have demonstrated that patterns and variability of SCS rather than the single TD or the lactation average are generally very informative on cow udder health and also on the type and severity of mastitis, particularly in the first half of lactation. Furthermore, it is important to recall that indicator traits like SCC must show genetic and phenotypic variation to be included in breeding programs and to identify the 'best' individuals in the population to be mated. Given this background, a similar approach is potentially exploitable in buffalo to investigate novel phenotypes for both managerial and breeding purposes to be used worldwide. Therefore, the present study aimed to investigate phenotypic variation of SCC-derived traits in buffalo milk and estimate their effects on milk production and composition in the first 150 DIM.

Data
A total of 1,414,449 TD records of 106,388 Italian water buffaloes (Bubalus bubalis) were provided by the Italian Breeders Association (Rome, Italy) for the period January 2013 to December 2017. All herds underwent official milk recording procedures according to International Committee for Animal Recording (2017). For each TD, information on buffalo ID, farm, birth date, date of calving, parity, DIM, milk yield (MY, kg/ day), fat (FP, %), protein (PP, %), lactose (LP, %) and SCC (cells/mL) was available. Composition traits (FP, PP and LP) were determined by mid-infrared spectroscopy using MilkoScan TM FT6000 (Foss Electric A/S, Hillerød, Denmark) according to ISO 21543 (2006) and SCC was assessed using Fossomatic TM FC (Foss Electric A/S, Hillerød, Denmark) according to ISO 13366 (2006). All the analyses were performed in the milk laboratory of the Breeders Association of Campania region (Benevento, Italy). According to Associazione Italiana Allevatori (2020) and International Committee for Animal Recording (2017), lactating buffaloes were sampled every 4 or 6 weeks during the lactation.
All the farms were located in the area of production of the Protected Designation of Origin mozzarella cheese, where lactating buffaloes are usually housed in free stalls with access to the paddock, are fed with total mixed ration, and are typically divided into different groups according to the stage of lactation.

Editing
Buffaloes were required to be in parity 1 to 6 and to have the first calving between 24 and 55 months of age. Animals calved between 2013 and 2017. Only TD with SCC from 1000 to 10,000,000 cells/mL were retained. As regards MY, FP, PP and LP, values outside the range mean ± 3 standard deviations (SD) were treated as missing. In the window from 5 to 400 DIM, the median of SCC was 70,000, 88,000, 94,000, 97,000, 99,000 and 96,500 cells/mL for parity 1, 2, 3, 4, 5 and 6, respectively, and to achieve normality, SCC was converted to SCS using the formula of Ali and Shook (1980): SCS ¼ 3 þ log 2 (SCC/100,000).
Subsequently, the focus was put on the first 150 DIM and thus TD above this limit were discarded and not considered for further investigation. Each lactation was required to have !5 TD SCC from 5 to 150 DIM, with the first available TD SCC within the first 35 DIM. The interval between two consecutive TD SCC was around 30 days. Finally, 232,669 TD from 45,312 lactations (5-150 DIM) of 35,623 buffaloes in 324 herds remained for subsequent statistical analysis.

Definition of SCC-derived traits
A description of the SCC-derived traits is provided in Table 1. In particular, mean (SCS 150 ) and SD (SCS_SD 150 ) of SCS in the first 150 DIM were calculated for each lactation as follows: where x i is the TD SCS, x is the lactation mean of all TD SCS within the lactation and n is the total number of TD SCS within the lactation. The severity (SEV 150 ) was calculated as follows: where n H is the number of TD SCC >200,000 cells/mL in the first 150 DIM and n T is the total number of TD SCC available in the first 150 DIM for each lactation. The threshold of 200,000 cells/mL was identified based on the existing literature (Dhakal 2006;Moroni et al. 2006;Tripaldi et al. 2010;Guccione 2013;. For each lactation, both mean and SD of MY, FP, PP and LP were calculated in a similar manner using the available TD (n ! 3) in the first 150 DIM.
Pearson correlations of SCS 150 and SCS_SD 150 with other traits were estimated, as well as Spearman correlations of SEV 150 with other traits. In fact, considering the way SEV 150 was calculated, ranking correlations were more appropriate. All statistical analyses were performed using the SAS software v.9.4 (SAS Institute Inc., Cary, NC).

Sources of variation of SCC-derived traits
The SCS 150 , SCS_SD 150 and SEV 150 were analysed according to the following random intercept multilevel linear model: where Y ijklmn is SCS 150 , SCS_SD 150 or SEV 150 ; l is the population mean; P i is the fixed effect of the ith parity (i ¼ 1 to 6); S j is the fixed effect of the jth season of calving (j ¼ winter -December to February, spring -March to May, summer -June to August, autumn -September to November); Y k is the fixed effect of the kth year of birth (k ¼ 2003 to 2015); (P Â S) ij is the fixed interaction effect between parity and season of calving; H l is the random intercept effect of the lth herd (l ¼ 1 to 324); B m (H l ) is the random intercept effect of the mth buffalo (m ¼ 1 to 35,623) nested within the lth herd; and e ijklmn is the random residual. The HPMIXED procedure of SAS software v.9.4 (SAS Institute Inc., Cary, NC) was used for the analysis of variance and the Bonferroni multiple comparisons post-hoc test was used to test if least squares means (LSM) of fixed effects differed significantly (p .05).

Effects of SCC-derived traits on milk yield and composition
In order to investigate the impact of SCS 150 , SCS_SD 150 and SEV 150 on the variability of mean and SD of MY and composition traits, SCC-derived traits were categorised and separately included as explanatory variables in a mixed linear model along with the same effects described in previous model. In particular, 5 classes of SCS 150 , 5 classes of SCS_SD 150 and 4 classes of SEV 150 were defined ( Table 2). The HPMIXED procedure of SAS software v.9.4 (SAS Institute Inc., Cary, NC) was used and the Bonferroni multiple Ratio between the number of test-day somatic cell count >200,000 cells/mL and the total number of test-day somatic cell count (!5) in the first 150 days in milk. a SCS 150 ¼ mean of test-day SCS of the first 150 days in milk; SCS_SD 150 ¼ SD of test-day SCS of the first 150 days in milk; SEV 150 ¼ severity, calculated as ratio between the number of test-day SCC >200,000 cells/mL and the total number of test-day SCC in the first 150 days in milk.
comparison post-hoc test was used to test if LSM of the levels of SCS 150 , SCS_SD 150 or SEV 150 differed significantly (p .05).

Data overview
The first TD SCC (<35 DIM) was >200,000 cells/mL in 15.3% of lactations, likely due to the high SCC which normally occurs immediately after calving.  observed that 65.98% of first-parity buffaloes had at least one TD SCC !200,000 cells/mL in the lactation, and this percentage increased up to 72.12% in fifth-parity animals. In the same study, authors reported that a relevant proportion of animals had high levels of SCC for almost the entire lactation; in particular, 28.36% first-parity and 39.05% sixth-parity buffaloes had average SCC !200,000 cells/mL in the lactation. As in other dairy species, authors are aware that the physiological baseline of SCC could change across lactating animals and also across herds. Therefore, the use of a fixed threshold (e.g. 200,000 cells/mL) may be questionable; nevertheless, it is important to highlight that the overall aim of farmers and the dairy industry is to reduce SCC or SCS to optimise milk technological traits and improve health status. Therefore, a potential selection index to improve udder health in buffalo should define new indicators able to account for both the average and the variability of SCC.
Descriptive statistics of SCC-derived traits, MY and composition traits are presented in Table 3. Average SCS 150 (2.80) corresponded to a SCC of 87,055 cells/ mL and, together with average SCS_SD 150 (1.25), mirrored findings reported by Bobbo et al. (2018) in firstcalving Holstein cows. In fact, those authors reported average SCS 150 and SCS_SD (5-305 DIM) of 2.66 and 1.29, respectively. The SD of SCS in the first 150 DIM may be an efficient indicator of anomalous values and peak occurrence, and thus can be useful to identify buffaloes with unstable (undesired) SCS, regardless of the SCC level. Among the SCC-derived traits, SEV 150 exhibited the greatest coefficient of variation (CV ;  Table 3). Again, these results agreed with the findings of Bobbo et al. (2018) in dairy cattle, where severity averaged 0.14 and had the greatest CV. About 47% of lactations did not show TD SCC >200,000 cells/mL in the first 150 DIM (class 'zero' of SEV 150 ; Table 2), meaning that the remaining 53% had at least 1 TD SCC >200,000 cells/mL. According to the threshold of 200,000 cells/mL proposed by Tripaldi et al. (2010), this means that there was at least one udder inflammation event in 53% of lactations. Milk yield averaged 10.75 kg/day and showed a CV of 23.14%, and means of FP, PP and LP were 7.67, 4.60 and 4.82%, respectively (Table 3). Mean FP and LP had the highest (11.77%) and the lowest CV (3.61%), in agreement with findings on buffalo and bovine milk-based on whole lactation (Tripaldi et al. 2010;Costa, Lopez-Villalobos, Visentin, et al. 2019;Costa, Negrini, De Marchi et al. 2020). Comparison with the literature is difficult for the SD of MY, FP, PP and LP, due to the  lack of large-scale studies that have investigated such traits in dairy species. The average SD of MY was 1.96 kg/day, with quite a high CV (47.91%). Concerning milk components, the lowest average SD was observed for LP (0.25%) and the greatest for FP (1.25%). Also, the greatest CV was calculated for LP (57.72%; Table 3). Parity-specific descriptive statistics are provided in Supplemental Table S1.

Correlations
The strongest correlation (0.85; p < .05) was estimated between SCS 150 and SEV 150 , whereas weak associations were assessed between SCS 150 and SCS_SD 150 (0.15; p < .05), and SCS_SD 150 and SEV 150 (0.28; p < .05, Table  4). In Holstein cows, Bobbo et al. (2018) estimated a slightly lower correlation of 0.70 between SCS 150 and severity, and slightly stronger associations of 0.20 between SCS 150 and SCS_SD, and 0.42 between severity and SCS_SD compared with our findings, likely due to differences in the temporal window (DIM) and species. Correlations suggest that each SCC-derived trait can provide additional information on the udder health status of animals. Correlations of SCC-derived traits with mean and SD of MY, FP, PP and LP were generally very weak (<0.10; Table 4); indeed, the strongest correlation (À0.22; p < .05) was observed between SCS 150 and LP. The inverse relationship between LP and SCS has been widely described in the literature for dairy species (Fox et al. 2015;Costa, Lopez-Villalobos, Sneddon, et al. 2019) and indicates that LP tends to drop in presence of high SCC (udder inflammation). Barbosa et al. (2019) estimated a similar correlation (À0.23) in Brazilian Murrah buffalo. In the study of  correlations were calculated within subsets created on the basis of average lactation SCC. In buffaloes with average lactation SCC <200,000 cells/mL the correlation between SCS and LP was À0.27 and the magnitude increased up to À0.35 in buffaloes with average SCC !500,000 cells/mL. Overall, both SCS 150 and SEV 150 were weakly negatively associated with means of MY, FP and PP, whereas SCS_SD 150 was weakly positively associated with FP and PP (Table 4). Such findings agreed with , except for the correlation between SCS 150 and MY which was slightly weaker in our study compared with  who used TD data rather than mean MY and SCS. Supporting this, Moroni et al. (2006) reported a weak correlation between SCS and MY in Italian water buffaloes. Furthermore, the same authors observed a slightly increasing trend of the magnitude of the correlation between MY and SCS at different levels of SCC, that is, buffaloes with average lactation SCC <200,000 (À0.10), !200,000 (À0.15), !300,000 (À0.17), !400,000 (À0.18) and !500,000 (À0.19).

Factors affecting SCC-derived traits
All fixed effects included in the statistical model significantly explained the variation of the SCC-derived traits, with the only exception of the interaction between parity and season of calving for SCS_SD 150 . In particular, p-values of fixed effects were <.001 for SCS 150 and SCS_SD 150 and <.01 for SEV 150 . The proportion of variance explained by the random effects indicated that the variation of the traits was more attributable to the herd effect rather than to the buffalo within-herd effect (Table 5); this could suggest that preventive rather than curative efforts should be allocated at the farm level and thus the key factor for a short-term efficient reduction of SCC is the management.
The LSM of SCS 150 , SCS_SD 150 and SEV 150 for the fixed effects of parity and season of calving are shown in Table 5. While SCS 150 and SEV 150 were the greatest in parity 5 and 6, and the lowest in parity 1, the maximum value of SCS_SD 150 (1.44 ± 0.02) was estimated for parity 1. This suggested that variability of SCS in buffalo decreases with parity, that is when the average SCS increases. As regards the season of calving, results highlighted that buffaloes calving in spring and summer had significantly lower SCS 150 (2.73 ± 0.05 and 2.76 ± 0.05, respectively) than buffaloes calving in Table 4. Correlations a (p .05) between somatic cell countderived traits b and mean and SD of milk yield and composition traits. À0.01 ns 0.14 À0.01 ns a Pearson correlations for SCS 150 and SCS_SD 150 and Spearman rank correlations for SEV 150 are presented. Superscript 'ns' indicates not significant correlations. b SCS 150 ¼ mean of test-day somatic cell score of the first 150 days in milk; SCS_SD 150 ¼ SD of test-day somatic cell score of the first 150 days in milk; SEV 150 ¼ severity, calculated as ratio between the number of test-day somatic cell count >200,000 cells/mL and the total number of test-day somatic cell count in the first 150 days in milk.
autumn and winter (Table 5). The opposite was observed for SCS_SD 150 , with the greatest values in spring and summer and the lowest in autumn and winter (Table 5). Finally, despite significance, the variation of SEV 150 across the four seasons was negligible ( Table 5). The trend of the SCC-derived traits across the year of birth of buffaloes is depicted in Figure 1. Overall, the results suggested that the average level of SCS is increasing in the population while the variability of SCS is slightly decreasing (Figure 1).

Impact of SCC-derived traits on milk yield and composition
Overall, MY and all composition traits were significantly affected by SCC-derived traits (p < .05), with only a few exceptions. The LSM of milk traits for the fixed effect of SCS 150 (5 classes) are depicted in Figure  2. Mean MY gradually decreased from class 1 (10.54 ± 0.09) to class 5 (9.83 ± 0.09), oppositely to MY SD, which exhibited the greatest estimates in classes 4 and 5 and the lowest in class 1 (Figure 2). The increase of SCS 150 resulted in a progressive reduction of mean LP from 4.89 ± 0.01 (class 1) to 4.72 ± 0.01 (class 5) and to a linear increase of PP. On the other hand, mean FP increased from class 1 to 3 of SCS 150 (Figure 2). These results generally agreed with previous findings of Cer on-Muñoz et al. (2002) and Tripaldi et al. (2010) in buffalo, andCosta, Egger-Danner, et al. (2019) Bansal et al. (2007) found LP to have the best discrimination ability for udder inflammation in 101 Murrah buffaloes. In cattle, Ebrahimie et al. (2018) observed similar findings using machine learning algorithms; in particular, lower LP in combination with lower MY was a pattern discovered by one of the decision trees built by the model. In this study, results suggest also that synthesis of fat and protein and thus FP and PP were not impaired by high SCC and/or udder inflammation. In other words, considering the concentration effect, the reduction of MY seems to be  responsible for the (apparent) increase of FP and PP in milk, as previously observed by Tripaldi et al. (2010) in buffalo. In support of this, SD of both FP and PP did not show a specific trend across classes of SCS 150 , whereas LP SD progressively increased. The LSM of mean MY and LP were the greatest in class 1 of SCS_SD 150 (10.42 ± 0.09 and 4.83 ± 0.01, respectively) and the lowest in class 5 (10.16 ± 0.09 and 4.78 ± 0.01, respectively; Figure 3). Mean PP slightly increased in a linear manner moving from class 1 (4.58 ± 0.01) to class 5 (4.61 ± 0.01). In general, it is important to highlight that a large variability of SCS is undesired in dairy species and thus animals with low SCS 150 and low SCS_SD 150 should be preferred as they exhibit optimal resistance to udder infections and prompt immune response. This is supported by findings of the present study; in fact, mean MY and LP were the greatest in class 1 of SCS_SD 150 , where their SD was minimum (Figure 3). In general, the SD of all milk traits was affected by SCS_SD 150 . Anyway, it is worth stressing that mastitis in buffaloes is often chronic or subclinical and thus mean and/or SD of SCS may not be exhaustively informative.
The LSM of mean MY and LP showed that these traits were affected similarly by SEV 150, with the greatest means estimated for the 'zero' class and the lowest for the 'high' class ( Figure 4). In other words, the greatest MY was estimated for buffaloes with all TD SCC <200,000 cells/mL during the first 150 DIM. On the contrary, buffaloes with more than half of TD SCC Figure 2. Least squares means of mean and SD of yield and composition traits of buffalo milk for the fixed effect of the mean somatic cell score in the first 150 days in milk (SCS 150 , 5 classes: from 1 ¼ low to 5 ¼ high). Different letters indicate significantly different means (p .05).
>200,000 cells/mL were the less productive and were characterised by the lowest milk LP. Mean FP and PP were less easy to be interpreted; for instance, mean PP in 'low' and 'high' classes were not significantly different, as well as mean PP in 'zero' and in 'intermediate' classes. Overall, for all traits, the SD was the lowest in class 'zero' but did not follow the expected linear pattern in 'low', 'intermediate' and 'high' classes ( Figure 4).

Conclusions
In the present study, milk SCC-derived traits were investigated in a large dataset of Italian water buffaloes to promote informative novel phenotypes for breeding and management purposes, and to quantify the impact of different levels of SCS on milk production and quality. Results showed that milk traits were impaired by high and variable SCS in the first 150 DIM. A negative association was observed between SCS 150 and LP, opening the debate on the validation of such trait for monitoring buffalo udder health in combination with other features. Finally, we demonstrated that high SCS 150 are translated into less MY. In perspective, further studies will validate these udder health indicators with clinical mastitis data and milk spectra to exploit routinely available information for genetic purposes and develop an udder health index for the dairy buffalo. Moreover, the effect of high SCS on the sensory characteristics of cheeses may also be investigated. Such findings on Italian buffalo may be extended to other populations worldwide and may open the debate on the need of implementing a specific recording scheme for udder diseases in the near future. In fact, the improvement of udder health/mastitis resistance at the population level usually takes a long time to be observed, particularly if artificial insemination is scarcely adopted and if productive life of lactating animals is long like in dairy buffalo.

Ethical approval
Data and milk spectra used in this study were collected during the routine milk recording, therefore the authorisation of the animal welfare and ethics committee was not required.