An attempt to identify milk protein fraction genotypes using unsupervised and supervised near-infrared spectroscopy methods

Abstract

The aim was to evaluate near-infrared spectroscopy (NIRS) potential to discriminate among β-casein (CN), κ-CN and β-lactoglobulin (LG) genotypes to be used as an authentication method. A total of 168 milk samples with known genetic information for β-CN, κ-CN and β-LG were collected at the same farm and paired with the NIRS spectrum. Spectra were evaluated with an unsupervised method (principal component analysis, PCA) and a supervised method (partial least squares-discriminant analysis, PLS-DA). For the PLS-DA, data were split into a train (75%) and a test set (25%), and the variable in projection >1 criterion was applied to select informative wavelengths. Results obtained confirmed that milk quality was similar among genetic variants. For the PCA, the observed variance explained by the first two principal components was 94%, but samples were not clustered by their genotypes of β-CN (i.e. A1A2, A2A2), κ-CN (i.e. AA, AB, AE, BB, BE) and β-LG (i.e. AA, AB, BB). The best accuracy for the PLS-DA models was reached by β-CN (train and test set, 64%), followed by β-LG (train set, 56%; test set, 52%) and κ-CN (train set, 41%; test set, 36%). In conclusion, the PCA on milk spectra was not able to cluster β-CN, κ-CN and β-LG genotypes, but the PLS-DA models revealed promising results for β-CN and β-LG. It could be interesting to increase the number of samples to equilibrate genetic variants and to apply a sampling selection method before discarding the applicability of NIRS as an authentication method.

Highlights

• Near-infrared spectroscopy discriminates β-casein (CN) and β-lactoglobulin (LG) more accurately than κ-CN genotypes.

• Scarce cluster ability of β-CN, κ-CN and β-LG genotypes with principal component analysis (PCA).

• Partial least squares-discriminant analysis (PLS-DA) moderately discriminates β-CN and β-LG genotypes.

• PLS-DA lowly discriminates κ-CN genotypes.

Keywords:

• A2 milk
• genetic polymorphism
• NIRS
• protein variants

Introduction

Nowadays, there is an interest in producing A2 milk to fulfil consumers’ demands because it is considered a healthier alternative to conventional A1 milk (Bodnár et al. 2018). The genetic variant of β-casein (CN) for milk A2 codes for the amino acid proline instead of histidine at position 67. This change prevents β-CN hydrolysis at this position and avoids peptide β-casomorphin-7 (β-CM-7) release after digestion (Thiruvengadam et al. 2021). The β-CM-7 has been associated with worsening gastrointestinal symptoms, cognitive traits and increasing intestinal transit time (Jianqin et al. 2016). Moreover, β-CN genotype variants impact milk production, protein yield and fat percentage (Bovenhuis et al. 1992). Genotypes of κ-CN affect protein yield and content (Bovenhuis et al. 1992; Tsiaras et al. 2005). Genetic variants of β-lactoglobulin (LG) have been associated with changes in whey composition and properties (Tsiaras et al. 2005; Rutten et al. 2011), milk production, protein yield and fat percentage (Bovenhuis et al. 1992). Thus, all these genetic variations impact milk quality and technological properties of milk (Čítek et al. 2021).

Genetic variants are usually identified through genetic tests (Caroli et al. 2009) which are highly accurate but expensive, labour-intensive and not feasible for large-scale screening (Xiao et al. 2022). Near-infrared spectroscopy (NIRS) is considered as a quick method, objective, non-destructive with the sample, free of chemical reagents, low cost and environmentally friendly analysis. It is routinely used in the dairy industry to predict milk quality and represents an interesting approach to detecting food adulteration (dos Santos Pereira et al. 2020). Few studies have attempted to apply principal component analysis (PCA; unsupervised method) or/and partial least squares-discriminant analysis (PLS-DA; supervised method) to infrared spectra to identify protein and whey milk fraction genotypes. Those studies used mid-infrared (MIR) spectroscopy instead of NIRS, and focused on β-CN and β-LG genetic variants. Daniloski et al. (2022) reported that the PCA plot failed to discriminate among β-CN genetic variants. On the other hand, Rutten et al. (2011) and Xiao et al. (2022) obtained a 74% for β-LG and 96% for β-CN accuracy, respectively, when applying a PLS-DA.

Therefore, this study aimed to evaluate NIRS potential to discriminate among β-CN, κ-CN and β-LG genetic variants.

Materials and methods

Milk samples analysis

A total of 168 individual milk (2 × 50 mL) samples from Holstein-Friesian cows were collected in plastic tubes with preservative bronopol (Broad Spectrum Micro-tabs II, D&F Control Systems, San Ramon, CA) in June 2022 from a commercial farm located in the area of Barcelona, Spain (42°00′11.2″N 2°13′16.8″E). This farm has been selecting milk A2 and has the specific genotype information for β-CN (i.e. A1A2, A2A2), κ-CN (i.e. AA, AB, AE, BB, BE) and β-LG (i.e. AA, AB, BB) provided by the Frisian Federation of Catalonia (FEFRIC, Barcelona, Spain). One aliquot (50 mL) was analysed by the Interprofessional Dairy Association of Catalonia (ALLIC) for gross composition (fat, protein, lactose and milk urea nitrogen (MUN)) by MilkoScan (FOSS, Hillerød, Denmark) and somatic cell count (SCC) by Fossomatic (FOSS, Hillerød, Denmark). As described by Manuelian et al. (2021), SCC was transformed into SCS by applying the Wiggans and Shook (1987) equation.

The other aliquot was analysed with a NIRSystems 5000 spectrophotometer (FOSS, Hillerød, Denmark) equipped with a scanning monochromator working from 1100 to 2500 nm every 2 nm at the Agriculture and Animal Production laboratory of the Universitat Autònoma de Barcelona (Barcelona, Spain). Before scanning the samples, they were heated at 40 °C for 10 min in a water bath (Aparatos Científicos, J.P. Selecta SA, Barcelona, Spain) and homogenised by gently rotating the bottles four times. Then, it was transferred into a quartz glass (diameter 4.8 cm and height 3.9 cm) where a gold reflector (0.5 mm path length) was placed. All samples were scanned in duplicate, manually shaking gently between each scan, and the average spectrum of each sample was used. The absorbance was recorded as a log (1/transflectance). Each spectrum was then matched with the corresponding milk protein genotype information.

Chemometric analysis

The chemometric analysis was conducted following Manuelian et al. (2021) procedure with R version 4.2.0 (R Core Team 2022). Briefly, ‘stats’ and ‘ggbiplot’ (Vu 2011) packages were used to perform the PCA. The ‘DiscriMiner’ package (Sanchez 2013) was applied to perform the PLS-DA for two components after removing the low signal-to-noise ratio wavelength, splitting the dataset into a train set (75% of the observations) and a test set (25% of the observations), and selecting the wavelengths with a variable importance in projection (VIP) score >1. The PLS-DA model performance was assessed by calculating the error rate and accuracy.

Statistical analysis

The sample size was calculated with G*Power software ver. 3.1.9.6 (Faul et al. 2007, 2009; Heinrich Heine Universität Düsseldorf, Düsseldorf, Germany). The normality of the traits was evaluated with PROC UNIVARIATE of SAS ver. 9.4 (SAS Institute Inc., Cary, NC). A multigene approach was applied to evaluate milk quality through a PROC GLM of SAS ver. 9.4 as suggested by Bovenhuis et al. (1992). The model included cow β-CN, κ-CN and β-LG genotypes, stage of lactation (i.e. class 1–4, being the first class from 3 to 100 DIM; class 2 from 101 to 200 DIM; class 3 from 201 to 300 DIM; and class 4 from >300 DIM), parity (i.e. primiparous and multiparous) and stage of lactation × parity as fixed effects. Results are reported as least squares means (LSM) and multiple comparisons were performed using Tukey’s adjustment when necessary. Significance was declared at p < .05 unless otherwise stated.

Results and discussion

Database description and statistics

Average fat and protein content (Table 1) are in line with Sola-Larrañaga and Navarro-Blasco (2009) who evaluated 348 bulk milk from herds located in northern Spain (Navarra region). Moreover, milk quality variability agreed with Franzoi et al. (2020), with Italian Holstein-Friesian individual records with similar MUN content, lower lactose content (4.76%), and greater fat (4.08%) and protein (3.31%) content, and SCS (2.59).

Table 1. Descriptive statistics of milk quality traits.

Trait^a	N	Mean	SD	Minimum	Maximum	CV, %
Milk composition
Fat, %	162	3.28	0.65	1.74	5.32	20
Protein, %	164	3.15	0.30	2.38	3.02	10
Lactose, %	161	5.18	0.28	4.31	5.84	5
MUN, mg/L	162	193.34	53.56	48.00	369.00	28
SCS, score	164	1.76	1.63	−0.32	7.69	93

^aMUN: milk urine nitrogen; SCS: somatic cell score.

For β-CN, the greater number of cows A2A2 (62.5%) than A1A2 (37.5%; Table 2) is linked to the selection strategy of the farm. For κ-CN, genotypes AB and AA were the most frequent ones (33.9% and 31.0%, respectively); and AE and BE were the least frequent ones (9.5% and 8.3%, respectively; Table 2). For β-LG, AA was the most frequent (51.2%) and BB was the least frequent (8.3%; Table 2). The unbalanced frequency of the genetic variants for κ-CN and β-LG could be attributed to favouring A2A2 for β-CN as indicated by Comin et al. (2008) and supports the multigene statistical approach (Bovenhuis et al. 1992). In agreement with our results, Comin et al. (2008) also reported a greater frequency of AA and AB for κ-CN when evaluating 1042 Italian Holstein cows. Some authors also observed that BE and EE genotypes for κ-CN are infrequent (Comin et al. 2008; Gai et al. 2021). In line with our results, the Gai et al. (2021) review indicated that genetic variant A is more frequent than B for β-LG in Holstein-Friesian cows.

Table 2. Least squares means (±SE) of milk quality¹ traits for the main effects.

Fixed effect^2	n	Fat, %	Protein, %	Lactose, %	MUN, mg/L	SCS
β-CN genotypes
A1A2	63	3.36 ± 0.10	3.20 ± 0.04	5.17 ± 0.04	200.3 ± 8.2	1.79 ± 0.23
A2A2	105	3.27 ± 0.10	3.16 ± 0.04	5.23 ± 0.04	196.2 ± 8.2	1.89 ± 0.23
κ-CN genotypes
AA	52	3.41 ± 0.11	3.16 ± 0.04^ab	5.11 ± 0.05	187.4 ± 8.9	2.103 ± 0.246
AB	57	3.30 ± 0.11	3.30 ± 0.04^a	5.14 ± 0.05	192.6 ± 9.1	2.099 ± 0.252
AE	16	3.36 ± 0.18	3.15 ± 0.07^ab	5.29 ± 0.08	206.9 ± 14.5	1.792 ± 0.404
BB	29	3.43 ± 0.14	3.27 ± 0.05^ab	5.20 ± 0.06	199.0 ± 11.1	1.492 ± 0.307
BE	14	3.10 ± 0.19	3.03 ± 0.08^b	5.25 ± 0.08	205.2 ± 15.8	1.732 ± 0.439
β-LG genotypes
AA	86	3.20 ± 0.08	3.13 ± 0.03	5.218 ± 0.031	193.4 ± 6.2	1.71 ± 0.17
AB	68	3.39 ± 0.09	3.12 ± 0.04	5.223 ± 0.037	212.4 ± 7.3	1.65 ± 0.20
BB	14	3.36 ± 0.18	3.29 ± 0.07	5.155 ± 0.076	188.9 ± 14.3	2.17 ± 0.40
Lactation stage
3–100 d	49	3.33 ± 0.12^ab	2.99 ± 0.05^c	5.21 ± 0.05	175.7 ± 9.4^b	1.64 ± 0.26
101–200 d	37	3.03 ± 0.12^b	3.12 ± 0.05^bc	5.26 ± 0.05	196.1 ± 10.0^ab	1.72 ± 0.28
201–300 d	51	3.23 ± 0.10^b	3.20 ± 0.04^b	5.16 ± 0.04	194.7 ± 8.2^ab	1.92 ± 0.23
>300 d	31	3.69 ± 0.13^a	3.42 ± 0.05^a	5.17 ± 0.05	226. ± 10.6^a	2.08 ± 0.30
Parity
Primiparous	87	3.24 ± 0.09	3.21 ± 0.04	5.30 ± 0.04^a	200.7 ± 7.3	1.19 ± 0.21^b
Multiparous	81	3.40 ± 0.09	3.15 ± 0.04	5.10 ± 0.04^b	195.8 ± 7.6	2.50 ± 0.20^a

¹ MUN: milk urine nitrogen; SCS: somatic cell score.

² β-CN: β-casein; κ-CN: κ-casein; β-LG: β-lactoglobulin.

^a,bSuperscript letters within trait and effect are significantly different (p < .05).

In the present study, milk quality parameter estimates (Table 2) did not differ among β-CN and β-LG variants. Whereas some authors have reported higher protein yield and lower fat percentage in A2 than A1 variants for β-CN, others have not identified differences in protein and fat percentage (Gai et al. 2021). Moreover, Cendron et al. (2021) also obtained similar fat, protein, MUN and SCS between A1A2 and A2A2 genotypes of β-CN in Italian Holstein. Nevertheless, Albarella et al. (2020) found significant differences in protein percentage in β-CN genotypes in the autochthonous agerolese cattle breed. In line with our results, Botaro et al. (2008) also described a similar milk composition among the β-LG genetic variants (i.e. AA, AB, BB) in Holstein cows. However, Cendron et al. (2021) reported greater fat percentage, MUN and SCS in genetic β-LG variant BB than AA. The lower protein content for BE than AB variants of κ-CN (p = .031) agreed with Cendron et al. (2021) results; however, we cannot discard that the observed differences were related to the lower number of cows with the BE variant. Although Heck et al. (2009) and Cendron et al. (2021) also reported similar milk fat percentage among κ-CN genotypes as we did, Cendron et al. (2021) reported significantly lower protein and MUN, and a greater SCS for AA than AB variants. In agreement with our results, Albarella et al. (2020) did not find differences in milk from the agerolese breed for AB and BB variants. Differences observed across lactation and parities are in agreement with the literature, and we are not going to discuss them further.

PCA and PLS-DA

The two first principal components (PC1 = 78.0%; PC2 = 16.3%) of the PCA applied to the NIR spectra explained 94.3% of the observed variance, but were unable to cluster β-CN, κ-CN and β-LG genotypes (Figure 1). Similar results were reported by Daniloski et al. (2022) for β-CN where overlapping clusters were observed when evaluating 114 milk samples from Australian Holstein-Friesian. However, those authors improved the PCA model’s performance by selecting only five samples of each variant that were separated enough to avoid overlaps. We are not aware of any other study trying to cluster samples based on protein fraction genetics applying PCA to infrared spectra.

Figure 1. Standardised plot of first (PC1) and second (PC2) principal components for β-casein (β-CN), κ-casein (κ-CN) and β-lactoglobulin (β-LG). The ellipse represents 95% confidence interval.

The average spectra of the samples presented a peak around 1200 nm, a slightly broad peak between 1380 nm and 1550 nm, a small peak around 1750 nm, and two broad peaks between 1860 and 2220 nm and between 2300 and 2498 nm (Figure 2). For β-CN, VIP was >1 from 1100 nm to 1116 nm; 1190 nm to 1222 nm; 1692 nm to 1774 nm; 1888 nm to 2002 nm; 2408 nm; 2410 nm; 2414 nm to 2422 nm; and 2430 nm to 2490 nm (Figure 2). For κ-CN, the VIP was >1 from 1100 nm to 1138 nm; 1888 nm to 2018 nm; and 2380 nm to 2498 nm (Figure 2). For β-LG the, VIP was >1 from 1708 nm to 1736 nm; 1880 nm to 2084 nm; and 2382 nm to 2498 nm (Figure 2). Those VIP peaks identified in all three fractions around 1940 nm have been linked to water (Coppa et al. 2010); around 2000 nm and 2400 nm to protein and fat content, respectively (Mohamed et al. 2021); and around 2466 nm to protein (Núñez-Sánchez et al. 2016). Moreover, those VIP peaks in β-CN and κ-CN fractions around 1700 nm have been associated with fat and protein content (Mohamed et al. 2021). Additionally, the VIP peak in β-CN fraction around 1200 nm has been associated with fat content (Núñez-Sánchez et al. 2016).

Figure 2. The absorbance of milk samples (blue), and the importance of the variable for projection (VIP) for β-casein, κ-casein and β-lactoglobulin for two components after removing signal-to-noise ratio wavelength. The red straight line indicates the threshold VIP = 1.

The best PLS-DA model based on the accuracy of the test set was observed for β-CN (64%), followed by β-LG (53%) and κ-CN (36%; Table 3). In contrast with our results, Xiao et al. (2022) obtained a greater accuracy in discriminating A1 and A2 milk in the test set (96%) using 754 MIR spectra of Chinese Holstein cows reared on the same farm. Although Daniloski et al. (2022) observed that the PLS-DA analysis with a subset of samples (30 out of 114) was able to separate milk samples collected in the same farm into different groups based on their β-CN variants, they stressed that this could be purely by chance due to the selected samples included in the final model. Similarly, Rutten et al. (2011) also reported a better accuracy for β-LG variants in the test set (74%) using MIR in 3826 milk samples from Dutch Holstein-Friesian cows collected in different farms. To the best of our knowledge, there are no studies regarding κ-CN genotype discrimination analysis.

Table 3. Partial least squares-discriminant analysis performance.

Fractions	Train set (n = 126)	Test set (n = 42)
β-CN genotypes
Error rate	0.36	0.36
Accuracy	0.64	0.64
κ-CN genotypes
Error rate	0.59	0.64
Accuracy	0.41	0.36
β-LG genotypes
Error rate	0.44	0.48
Accuracy	0.56	0.52

β-CN: β-casein; κ-CN: κ-casein; β-LG: β-lactoglobulin.

Conclusions

Our results demonstrated the difficulties of using NIRS to cluster milk samples with an unsupervised method such as the PCA based on their genotypes for β-CN, κ-CN and β-LG when milk quality is similar. The supervised method PLS-DA led to slightly better results depending on the fractions evaluated, revealing promising results for β-CN and β-LG. The best accuracy was achieved for β-CN, reaching up to 64% when dealing with A1A2 and A2A2 milk. The accuracy decreased for β-LG, but it still was above 50% considering the three variants. The worst results were obtained for κ-CN, which included five variants. We cannot discard that the decrease in the accuracy of the models could be related to the unbalanced number of samples for the different variants of each fraction. Further research trying to equilibrate variants’ representativeness within each fraction could be advisable before discarding NIRS as an authentication method. Moreover, it could be interesting to apply sample selection techniques to identify samples with different genetic variants that are farther in the PCA to also improve the models.

Ethics statement

Procedures adopted in the present study do not fall into the scope of an animal ethics evaluation.

Acknowledgements

The authors thank the Interprofessional Dairy Association of Catalonia (ALLIC) (Barcelona, Spain) for providing the samples.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data presented in this study are available free of charge for any user upon reasonable request from the corresponding authors.

Additional information

Funding

This research was funded by the Ministry of Science and Innovation of Spain with the Project PID2019-110752RB.I00. C.L. Manuelian is currently a postdoctoral researcher funded with a María Zambrano Grant from the Spanish Ministry of Universities (funded by European Union-Next Generation EU; MZ2021-86).