A comprehensive integration of factors affecting vitamin B 12 concentration in milk of Holstein cows: genetic variability, milk productivity, animal characteristics, and feeding management

: Daily vitamin B 12 (VB 12 ) requirements of humans can naturally be fulfilled by animal product consumption, especially products from ruminants because of bacteria dwelling in their rumen. Indeed, only bacteria can synthesize this vitamin. Milk is hence an excellent source of VB 12 . This cross-sectional study was undertaken to unravel factors, such as genetic variation, diet and cow characteristics, and milk production, explaining the large variation in milk B 12 concentration among cows by using an integrative approach. Milk samples from 2 consecutive milkings were collected from 3,533 Canadian Holstein cows (1,239 first, 932 s


Short Communication Genetics
The list of standard abbreviations for JDSC is available at adsa.org/jdsc-abbreviations-24.Nonstandard abbreviations are available in the Notes.
Abstract: Daily vitamin B 12 (VB 12 ) requirements of humans can naturally be fulfilled by animal product consumption, especially products from ruminants because of bacteria dwelling in their rumen.Indeed, only bacteria can synthesize this vitamin.Milk is hence an excellent source of VB 12 .This cross-sectional study was undertaken to unravel factors, such as genetic variation, diet and cow characteristics, and milk production, explaining the large variation in milk B 12 concentration among cows by using an integrative approach.Milk samples from 2 consecutive milkings were collected from 3,533 Canadian Holstein cows (1,239 first, 932 s, and 1,362 third and more lactations) located in 99 herds with various feeding management.For genetic variation analysis purpose, pedigrees were traced back for 3 complete generations for each sire and dam.A total of 10,021 identities were used in the subsequent genetic analyses.Milk VB 12 averaged 4.2 ng/mL with a range between 0.7 and 9.0 ng/mL.Dietary fiber (NDF from forage, dietary NDF, ADF, and lignin) increased and dietary components related to energy (non-fiber carbohydrate, starch, net energy of lactation, and percentage of concentrate) decreased VB 12 in milk.Milk VB 12 varied with days in milk, with a similar pattern as milk fat and protein concentration lactation curves.Milk VB 12 increased as age at calving increased.When disregarding the herd variance, heritability value was 0.37, meaning that milk VB 12 can be modified by genetic selection.The final model including factors related to the diet, animal characteristics and milk productivity, and genetic variation explained 79% (pseudo-R 2 ) of the milk VB 12 variation.When excluding the random effect of the cow, i.e., excluding the animal and genetic relationships, the pseudo-R 2 dropped to 43%, reinforcing the importance of genetic variation in explaining milk VB 12 variation.To our knowledge, the present study is the most comprehensive evaluation of factors affecting milk VB 12 variation including the greatest number of cows from various lactation stages.
V itamin B 12 (VB 12 ) is essential for animals and humans.In humans, this vitamin is important for preventing neuropathy, megaloblastic anemia, and gastrointestinal problems (Combs and McClung, 2022).Vitamin B 12 differs from other B vitamins as it is only synthesized by bacteria when cobalt is available (Martens et al., 2002).Hence, as the rumen contains a myriad of bacteria synthesizing VB 12 , ruminant products are an excellent source of this vitamin for humans.Indeed, cow milk is a good dependable and natural source of VB 12 (Walther et al., 2022); according to the National Academy of Sciences (1998), a 250-mL glass of milk provides 46% of the recommended dietary allowance of VB 12 for humans above 13 years old.A recent study has shown that milk VB 12 concentration was highly variable among herds; a 250-mL glass of milk fulfilled from 28 to 61% of the recommended dietary allowance (Duplessis et al., 2019).Given this wide variability and the importance of VB 12 for humans, some authors studied factors explaining its variation in milk.It has been previously observed that genetic effects (Rutten et al., 2013;Duplessis et al., 2016;Gebreyesus et al., 2021), milk productivity and cow characteristics (Duplessis et al., 2019), rumen microbiota (Franco-Lopez et al., 2020), cow breed (Anthony et al., 1951;Duplessis et al., 2018), cow nutrition (Duplessis et al., 2016;2019), rumen fermentation parameters (Duplessis et al., 2021), and seasons (Scott et al., 1984;Duplessis et al., 2016) have an impact on milk VB 12 .It shows that the regulation of milk VB 12 concentration is multifactorial.
Increasing the knowledge on factors affecting milk VB 12 variation is the first step to optimize and stabilize its concentration, hence strengthening the perception that milk is an excellent source of VB 12 for consumers.This is particularly important as milk consumption has constantly declined over the last decades (Stewart et al., 2021).Thus, the primary objective of this study was to build a model that integrates several factors (genetic and non-genetic) affecting milk VB 12 from a large data set of Canadian Holstein cows.Another rationale was to evaluate the heritability (h 2 ; i.e., the proportion of variation caused by genetic variation) of milk VB 12 of Holstein dairy cows from a large data set containing animals with a broad range of parities and DIM located in several herds.
All procedures were approved by the Animal Care Committee from Université Laval, QC, Canada following the guidelines of the Canadian Council on Animal Care (2009).Animals and data collection were previously described by Duplessis et al. (2019).Briefly, 100 commercial dairy herds located in Québec, Canada were involved in this cross-sectional study.To participate, it was required that Holstein was the main breed of the herd, that performance data were recorded through DHI agency, and that cows were milked twice daily.Herds that fulfilled these requirements were contacted by phone and participation was on a voluntary basis.Cows were housed in tie (n = 98) or free (n = 2) stalls and herd size ranged from 16 to 113 lactating Holstein cows.Three different feeding systems were used: 1) TMR (n = 31); 2) Automatic component feeding (ACF; n = 49); and 3) Manual component feeding (MCF;n = 20).
Herds were visited during 3 consecutive milkings (i.e., morning and evening milkings the first day and morning milking the second day of the visit).Each ingredient fed to the cows was sampled during the first morning milking visit and stored at −20°C until further analysis.Daily quantity of each ingredient given was recorded as previously reported (Duplessis et al., 2019).During evening milking of the first day and morning milking of the second day, milk yield was recorded, and milk samples were separately collected for each cow using calibrated in-line milk meters throughout the milking.A total of 4,340 Holstein cows (1,484 first, 1,093 s, and 1,763 third and more lactations) from 3 to 592 DIM were milk sampled.One aliquot of milk per milking was preserved with bronopol for major component analyses and another one, without bronopol preservative, was kept at −20°C until VB 12 analysis.Characteristics of cows (i.e., parity, DIM) and registration numbers were obtained from Lactanet (Canadian Network for Dairy Excellence, Sainte-Anne-de-Bellevue, QC, Canada).Body weight was estimated using heart girth circumference measurement (Yan et al., 2009).
Milk samples with bronopol from evening and morning milkings were immediately sent to the Lactanet laboratory for milk fat, protein, and lactose concentration analyses by mid-infrared spectroscopy using MilkoScanTM FT 6000 (Foss, Hillerød, Denmark).Milk VB 12 of evening and morning milkings were analyzed in duplicate by radioassay using a commercial kit (SimulTrac B 12 / Folate-S; MP Biomedicals, Solon, OH Duplessis et al. (2015)).The inter-assay CV was 3.1%.Daily milk VB 12 was calculated as follows: (evening milk VB 12 concentration × evening milk yield + morning milk VB 12 concentration × morning milk yield)/daily milk yield.Vitamin B 12 was obtained in pg/mL and transformed into ng/ mL to facilitate numerical stability.
Feed samples were dried at 55°C for 48 h, ground at 1 mm, and sent at SGS Canada (Guelph, ON, Canada) for nutrient concentrations based on AOAC International (2005) as previously described (Fadul-Pacheco et al., 2017).Percentage of inclusion of each ingredient in the diet was computed on a DM basis.Nutrient composition of the diet was then calculated by multiplying the inclusion percentage of each ingredient on a DM basis by its respective nutrient composition.
After merging the Lactanet registration number file of the 100 dairy herds with sample collection data set file based on cow identification numbers, there were 4,672 records in the data set.Data edition was then needed to remove irrelevant records.After the edition for duplicate, mismatched or invalid registration number, there were 4,407 records remaining.When merging with the daily milk VB 12 file, 278 records were excluded due to missing milk VB 12 data.Only animals with age at calving between 21 and 96 mo were kept (158 records excluded).There were insufficient numbers of animals in parities 7 and greater to warrant retaining them and they were dropped from the analyses (n = 89).Records with no milk production information (n = 27) as well as milk VB 12 greater than 9.0 ng/mL were deleted (n = 29) as they were considered as outliers or coming from unhealthy cows (Duplessis et al., 2018).Records were dropped because animal identification was inconsistent from the registration number file of Lactanet and the daily milk VB 12 file (n = 90).From the latter, one herd was dropped because all records were inconsistent (n = 34).All records from cows which were not pure Holstein or one of their parents were not pure Holstein were removed from the data set (n = 127).Finally, records of cows with more than 430 DIM were deleted (n = 76).The final data set had 3,533 records (1,239 first, 932 s, and 1,362 third and more lactations), one record per cow from 99 herds.From these, pedigrees (ancestors) were traced back for 3 complete generations from each of their sires and dams using the Holstein breed pedigree files maintained by Lactanet.This produced a total of 12,980 animal identities.Non-informative animals were pruned out; i.e., animals who had both parents unknown, no actual data records themselves and only one progeny.There remained 10,021 identities for use in the subsequent genetic analysis.
Descriptive statistics were obtained using Proc UNIVARIATE of SAS (version 9.4; SAS Institute, 2012).The current assessment uses a subset of cows from the study of Duplessis et al. (2019) in which the authors concluded that milk production, diet, and cow characteristics explained 52% of the milk VB 12 variation.However, they did not include genetic relationships among animals in their model, and more exclusion criteria as described above were applied in the current work leading to a suppression of 807 cows relative to the data set of Duplessis et al. (2019).These authors used principal component (PC) analysis to overcome multicollinearity among independent variables related to dietary composition and ingredients.Then, they used PC as independent variables in a multivariable model aiming to study the relationship between daily milk VB 12 as dependent variable.Thus, a similar approach as Duplessis et al. (2019) with slightly different independent variables was used in the current paper to also include the impact of genetic variation on VB 12 in the model for a better integration.Briefly, 26 independent continuous variables pertaining to diet composition [NE L (Mcal/kg of DM), CP, RUP, starch, ADF, aNDF, NDF from forage, lignin, NFC, fat (% of DM), and Co (mg/kg of DM)], diet ingredients [corn silage, chopped or baled mixed grass-legume silage, hay, concentrate, corn, other cereals (i.e., oat, barley and wheat), soy products, commercial protein supplement, commercial energy supplement, and fat supplement (% of DM)], and milk production (kg/day) and composition (%; fat, protein, and lactose concentrations) were used in the PC analysis using Proc FACTOR of SAS.To be retained for further analysis, a PC should explain at least 5% of the total variance.A VARIMAX rotation was applied to the PC to facilitate their interpretation.To be part of a PC, a variable should have the highest coefficient among the obtained PC.Then, each PC can be considered as an independent variable in a multiple regression analysis without multicollinearity among the PC.From the pedigree (3 generations back on both the sire and dam pathways of each sampled cow), inbreeding values, which are associated with a loss of genetic diversity (Meuwissen et al., 2020), were calculated according to Quaas (1976).
The model used in WOMBAT (Meyer, 2007) to estimate genetic parameters in addition to the PC and other independent factors was: Where Y ijklmpstuw is the dependent variable, i.e., daily VB 12 in milk collected from the wth cow from the ith feeding system (FS) in the jth month of sampling, kth herd nested within the ith FS in the jth month, which cow was sampled when she was at her lth DIM at mth parity and calved at the pth age at calving (AaC) in months nested within parity, with inbreeding level s with tth BW and uth PC; • µ is the overall mean, • FS i = the fixed effect of FS (TMR, ACF, and MCF), • month = the fixed effect of sampling month, • herd ijk = the random effect of the kth herd nested within the ith FS and the jth month of sampling [k = 1 to 99; herd ijk ~N(0, σ2 )], • DIM l = the fixed effect of DIM of the cow when the sample was collected (l = 1 to 430), • e (−0.05*DIM) l is to take into account the Wilmink curve model (Wilmink, 1987), • parity m = the fixed effect of parity of the sampled cows (m = 1 to 6), • AaC mp = the fixed effect of AaC (P = 21 to 96 mo) nested within parity, • inbreeding s = the fixed effect of the inbreeding level category (increment of 1%; s = 1 to 15%), • BW t = the fixed effect of estimated cow BW, • PC u = the fixed effect of PC (u = 1 to 6), • animal ijklmpstuw = the random effect of the cow to include animal and genetic relationships ~N(0, Aσ 2 ), where A is the numerator relationship matrix, and • e ijklmpstuw is the residual error term.
Proc REG of SAS was used to compute pseudo-R 2 from predicted milk VB 12 from the final model and observed values.The predicted values were obtained from WOMBAT.The residuals were visually assessed for normality and homoscedasticity.Results were considered significant when P ≤ 0.05.Heritability was computed as per Poulsen et al. (2015) and Gebreyesus et al. (2021), including or not the herd variance.
On average, cows in the current assessment yielded 32.4 kg/d of milk (Table 1).Milk VB 12 ranged from 0.68 to 8.96 ng/mL and averaged 4.18 ± 1.45 ng/mL (Table 1).A typical ration had a forage: concentrate ratio averaging 67:33 with 15.2% of crude protein, 1.57Mcal/kg of net energy of lactation, 21.8% of ADF, 37.4% of NDF, 40.7% of NFC, and 0.60 mg/kg of Co on a DM basis.Among the 10,021 identities used for subsequent genetic analyses, inbreeding values ranged between 0 and 15%; 83% of the animals had an inbreeding value between 0 and 1%, 10% between 1 and 2% and the remaining 7% of animals had inbreeding coefficients above 2%.Current inbreeding values were lower than reported for Canadian Holstein cattle born from 1980 to 2004 (3.20%) (Sewalem et al., 2006).
There were 6 PC that fulfilled the requirement of explaining at least 5% of the total variance.These 6 PC were then retained for further analysis and altogether explained 59% of the total variance; 22, 11, and 8% of the variance were explained by the PC1, PC2, and PC3, respectively.The first PC was negatively associated with dietary fiber (NDF from forage, NDF, ADF, and lignin) and positively with indicators of dietary energy (NFC, starch, NE L , and percentage of concentrate).Lactation performance was grouped into the PC2 (negatively with milk yield and positively with milk fat and protein concentrations).The third PC was positively related to dietary protein (CP and RUP) and Co concentrations.The fourth PC was negatively associated with mixed silage as bale, whereas it was positively associated with chopped mixed silage stored in silo.Dietary fat and percentages of commercial fat inclusion were positively associated while percentage of corn silage was negatively related with PC5.The last PC, PC6, had a positive relationship with percentage of soya products.
The novelty of the current assessment is to allow integrating the impact of genetic variation, animal characteristics and milk productivity, and diets on milk VB 12 in the same multivariable statistical model.It is noteworthy that the independent variables that affected VB 12 in milk are similar to results reported by Duplessis et al. (2019); the current assessment using a subset of their cows and a comparable methodology.Nevertheless, the random effect of the cow to include animal and genetic relationships and the fixed effect of inbreeding were not previously added.
Table 2 depicts which independent variables had a significant relationship with daily milk VB 12 .Milk VB 12 differed according to DIM and the exponential term of DIM (P < 0.001), i.e., the effect of DIM was well-modeled with the Wilmink regression.Age at calving and estimated BW were positively related with VB 12 in milk (P < 0.001).Out of 6 PC, 2 were retained in the final model because they reached significance.The one pertaining to dietary The h 2 was 0.27 ± 0.04 when considering the herd variance and 0.37 ± 0.05 without the herd variance (Table 2).To our knowledge, the current report is the most comprehensive assessment on genetic variation in milk VB 12 , including more than 3,500 Canadian Holstein cows at all lactation stages and parities.The first report demonstrating a genetic variation among animals for milk VB 12 used 544 primiparous Dutch Holstein cows between 76 and 282 DIM (Rutten et al., 2013).Other previous genetic evaluations of VB 12 in milk included 343 early-lactating Canadian Holstein (Duples-sis et al., 2016) and 341 Danish Holstein cows (Gebreyesus et al., 2021).In the literature, 2 different h 2 calculations have been found (Rutten et al., 2013;Poulsen et al., 2015), considering or not the herd variance.We obtained a greater value when the herd variance was not considered as per Poulsen et al. (2015).Our results are in agreement with previous experiments reporting an h 2 of 0.37 when disregarding the herd variance in the calculation (Rutten et al., 2013;Gebreyesus et al., 2021).Moreover, Duplessis et al. (2016) estimated the h 2 value to 0.23 when the herd variance was included in the calculation, similar to the 0.27 of the current study.These studies showing moderate h 2 for milk VB 12 indicate that there is a significant additive genetic variability, which could be exploited to increase milk VB 12 supplied to the consumer.Our results reported for the first time that inbreeding did not have an impact on milk VB 12 , although current inbreeding values were lower than previously reported in Canadian Holstein.
The final WOMBAT model including significant independent variables with the random effect of the cow to include animal and genetic relationships explained 79% of the milk VB 12 variability (pseudo-R 2 of 79%).Not including the animal and genetic relationships gave a pseudo-R 2 of 43% (i.e., animal and genetic explained about 36% of the variation) and further excluding DIM and e (−0.05*DIM) gave a pseudo-R 2 of 27% (i.e., DIM and e (−0.05*DIM) explained about 16% of the variation).Interestingly, PC2, which was the third largest source of the variation, only contributed to 6% of the variation, and PC1, 1% of the variation.Duplessis et al. (2019) concluded that 52% of the milk VB 12 variability was explained by milk production, cow characteristics, and diet composition, indicating a substantial improvement by including the animal and genetic relationships.
In summary, these results highlight that several factors affect milk VB 12 .With a h 2 of 0.37 when not considering the herd variance, this strengthens the idea that VB 12 in milk can be increased by genetic selection.The impact of other factors such as ruminal   According to the calculation of Poulsen et al. (2015). 5 Variance components: σ 2 a = genetic variation of the animal (0.421 ± 0.065); σ 2 herd = herd variation (0.384 ± 0.063); σ 2 e = residual variation (0.739 ± 0.053).
microbiome and breed of the cow on VB 12 can be further explored in an integrated model from data collected in several herds.

1PC1=
Negatively associated with dietary fiber (NDF from forage, NDF, ADF, and lignin) and positively with indicators of dietary energy (non-fiber carbohydrate, starch, net energy of lactation, and percentage of concentrate).2 PC2 = Negatively associated with milk yield and positively with milk fat and protein concentrations.3 Sum of the three variance components: genetic variation of the animal, herd variation, and residual variation.

Table 1 .
Duplessis et al. (2019)2019)Vitamin B 12Variation Descriptive statistics regarding the 3,533 Holstein dairy cows located in 99 herds included in the analysis Obtained from evening and morning milk vitamin B 12 concentration samples weighted according to respective milk yield.fiber and dietary indicators of energy (PC1) was negatively associated and the one related to lactation performance (PC2) was positively related to milk VB 12 (P < 0.001).Results indicate that milk VB 12 increased as dietary ADF, NDF, NDF from forages, and lignin increased, whereas an inverse relationship was obtained for dietary NE L , starch, NFC and concentrate percentages as similarly observed byDuplessis et al. (2016;2019).vandenOeveretal.(2021)observedthatmilk VB 12 was greater when cows ingested silage instead of hay as the sole forage.Our results did not support this finding, possibly because all herds fed their cows with a blend of 2 or more different forages in contrary with the study of van den Oever et al. (2021).The relationship with PC2 showed that milk yield was negatively and milk fat and protein concentrations were positively related with milk VB 12 .Hence, cows producing low milk yield with high milk fat and protein concentrations will have greater milk VB 12 than the opposite.Results on the relationship between DIM, exponential term of DIM, milk performance, and milk VB 12 are similar to the ones reported inDuplessis et al. (2019), in which the lactation curve of VB 12 reached the nadir at 50-55 DIM.Results suggest that as AaC increased, milk VB 12 also increased.There was no significant effect of parity over and above the effect of age at calving (P = 0.43); parity was therefore dropped and age at calving was kept in the model.It is noteworthy that age at calving is closely related to parity.Similarly, in the used byDuplessis et al. (2019), milk VB 12 increased as parity increased.

Table 2 .
Duplessis et al. | Cow Milk Vitamin B 12Variation Parameter estimates included in the multivariable model and heritability of daily milk vitamin B 12 concentration (ng/mL) of 3,533 Canadian Holstein cows located in 99 herds