Advertisement

Opportunities for genomic selection of cheese-making traits in Montbéliarde cows

Open AccessPublished:April 18, 2022DOI:https://doi.org/10.3168/jds.2021-21558

      ABSTRACT

      As part of the From'MIR project, traits related to the composition and cheese-making properties (CMP) of milk were predicted from 6.6 million mid-infrared spectra taken from 410,622 Montbéliarde cows (19,862 with genotypes). Genome-wide association studies of imputed whole-genome sequences highlighted candidate SNPs that were then added to the EuroG10K BeadChip, which is routinely used in genomic selection. In the present study, we (1) assessed the reliability of single-step genomic BLUP breeding values (ssEBVs) for cheese yields, coagulation traits, and casein and calcium content generated from test-day records of the first 3 lactations, (2) estimated realized genetic trends for these traits over the last decade, and (3) simulated different cheese-making breeding objectives and estimated the responses for CMP as well as for other traits currently selected in the Montbéliarde breed. To estimate the reliability of ssEBVs, the available data were split into 2 independent training and validation sets that respectively contained cows with the oldest and the most recent lactation data. The training set included 155,961 cows (12,850 with genotypes) and was used to predict ssEBVs of 2,125 genotyped cows in the validation set. We first tested 4 models that included either lactation (LACT) or test-day (TD) records from the first (1) or the first 3 (3) lactations, giving equal weight to all 50K SNP effects. Mean reliabilities were 61%, 62%, 63%, and 64% for the LACT1, LACT3, TD1, and TD3 models, respectively. Using the most accurate model (TD3), we then compared the reliabilities of 3 scenarios with: SNPs from the Illumina BovineSNP50 BeadChip only, equally weighted (50K); 50K SNPs plus additional candidate SNPs, equally weighted (50K+); and 50K and candidate SNPs with additional weight given to 7 to 14 candidate SNPs, depending on the trait (CAND). The 50K+ and CAND scenarios led to similar mean reliabilities (67%) and both outperformed the 50K scenario (64%), whereas the CAND scenario generated the less biased ssEBVs. To assess genetic trends, SNP effects were estimated with a single-step GBLUP based on the TD3 model and the 50K scenario applied to the whole population (2.6 million performance records from 190,261 cows and 423,348 animals in the pedigree, of which 21,874 were genotyped) and then applied to 50K genotypes of 21,171 males and 311,761 females. We detected a positive genetic trend for all CMP during the last decade, probably due to selection for an increase in milk protein and fat content in Montbéliarde cows. Finally, we compared the selection responses to 3 different breeding objectives: the current Montbéliarde total merit index (TMI) and 2 alternative scenarios that gave a weight of 70% to TMI and the remaining 30% to either milk casein content (TMI-COMP) or a combination of 3 CMP (TMI-Cheese). The TMI-Cheese scenario yielded the best responses for all the CMP analyzed, whereas values in the TMI-COMP scenario were intermediate, with a slight effect on other traits currently included in TMI. Based on these results, a program of genomic evaluation for CMP predicted from mid-infrared spectra was designed and implemented for the Montbéliarde breed.

      Key words

      INTRODUCTION

      More than 36% of total cow milk is processed into cheese products worldwide (
      • International Dairy Federation
      Dairy's global impact.
      ) and this proportion has increased by 23% during the last decade. The milk processing industry thus stands to gain a great deal economically from the improvement of milk cheese-making properties (CMP). To assess CMP, various laboratory methods have been developed (reviewed in
      • Bittante G.
      • Penasa M.
      • Cecchinato A.
      Invited review: Genetics and modeling of milk coagulation properties.
      ), but they are costly and time consuming, and thus difficult to implement on a large scale. Mid-infrared (MIR) spectrometry has been proposed as an alternative method for the prediction of various milk characteristics, including cheese yield or coagulation properties (
      • De Marchi M.
      • Toffanin V.
      • Cassandro M.
      • Penasa M.
      Invited review: Mid-infrared spectroscopy as phenotyping tool for milk traits.
      ). Mid-infrared technology is cheap and routinely used. With this approach, CMP can be indirectly measured on a large scale and could therefore be incorporated into the breeding objective for dairy cattle.
      The From'MIR project, initiated in 2015, aims to analyze CMP and milk composition traits predicted from MIR spectra in Montbéliarde cows from the Franche-Comté region; the milk from these cows is mainly used to produce a variety of high-quality cheeses. As part of the work of From'MIR, MIR prediction equations were developed (
      • El Jabri M.
      • Sanchez M.P.
      • Trossat P.
      • Laithier C.
      • Wolf V.
      • Grosperrin P.
      • Beuvier E.
      • Rolet-Répécaud O.
      • Gavoye S.
      • Gauzere Y.
      • Belysheva O.
      • Notz E.
      • Boichard D.
      • Delacroix-Buchet A.
      Comparison of Bayesian and partial least squares regression methods for mid-infrared prediction of cheese-making properties in Montbéliarde cows.
      ) and a set of 8 CMP traits (3 laboratory cheese yields and 5 coagulation traits) was predicted with relatively high accuracy from 6.6 million MIR spectra taken from 410,622 Montbéliarde cows (19,862 with genotypes). In prior studies, we estimated genetic parameters and conducted GWAS on whole-genome sequences for CMP and milk composition traits (proteins, fatty acids, and minerals). These studies revealed medium-to-high heritabilities for CMP traits and high genetic correlations among CMP traits and between CMP and some milk composition traits (
      • Sanchez M.P.
      • El Jabri M.
      • Minéry S.
      • Wolf V.
      • Beuvier E.
      • Laithier C.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Genetic parameters for cheese-making properties and milk composition predicted from mid-infrared spectra in a large dataset of Montbéliarde cows.
      ), as well as several genomic regions (QTL) with large effects on these traits (
      • Sanchez M.P.
      • Ramayo-Caldas Y.
      • Wolf V.
      • Laithier C.
      • El Jabri M.
      • Michenet A.
      • Boussaha M.
      • Taussat S.
      • Fritz S.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows.
      ). Candidate causative variants were selected from the sequence-based GWAS and incorporated into the EuroG10K BeadChip, which was developed by the EuroGenomics consortium (https://www.eurogenomics.com/actualites/the-eurog-md:-a-unique-genotyping-microarray-for-cattle-.html) and is routinely used in genomic selection.
      The objectives of the present study were to (1) estimate the reliability of single-step GBLUP breeding values for cheese yields, coagulation traits, and casein and calcium content in milk using different data sets and models, (2) estimate realized genetic trends for these traits in recent years, and (3) simulate different breeding objectives to estimate responses for CMP traits as well as for other traits currently selected in the Montbéliarde breed.

      MATERIALS AND METHODS

      Reference Population: Animals, Genotypes, and Phenotypes

      The original data set is described in detail in prior studies (
      • Sanchez M.P.
      • El Jabri M.
      • Minéry S.
      • Wolf V.
      • Beuvier E.
      • Laithier C.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Genetic parameters for cheese-making properties and milk composition predicted from mid-infrared spectra in a large dataset of Montbéliarde cows.
      ,
      • Sanchez M.P.
      • Ramayo-Caldas Y.
      • Wolf V.
      • Laithier C.
      • El Jabri M.
      • Michenet A.
      • Boussaha M.
      • Taussat S.
      • Fritz S.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows.
      ). It consisted of 4,877,908 MIR spectra of milk samples from 311,613 Montbéliarde cows originating from 3,229 herds of the Franche-Comté region. To ensure that the data set was homogeneous in the present study, we retained cows whose age at first calving ranged between 22 and 44 mo and who had a complete first lactation and at least 3 test-day records. After filtering, the data set included 2,869,353 test-day records from 191,532 cows.
      For this study, we did not perform any experiment on animals; thus, no ethical approval was required.
      Of these cows, 19,564 were genotyped for the purpose of genomic selection using either the BovineSNP50 (50K, 6,498 cows) or the EuroG10K BeadChip (Illumina Inc., 13,066 cows). The latter chip contains generic SNPs as well as a research add-on for causal or predictive SNPs for traits of interest in cattle, including candidate causative SNPs identified for milk composition (
      • Boichard D.
      • Govignon-Gion A.
      • Larroque H.
      • Maroteau C.
      • Palhiere I.
      • Tosser-Klop G.
      • Rupp R.
      • Sanchez M.P.
      • Brochard M.
      • Amigues Y.
      • Boscher M.Y.
      • Leveziel H.
      Genetic determinism of milk composition in fatty acids and proteins in ruminants, and selection potential.
      ;
      • Sanchez M.P.
      • Govignon-Gion A.
      • Croiseau P.
      • Fritz S.
      • Hozé C.
      • Miranda G.
      • Martin P.
      • Barbat-Leterrier A.
      • Letaïef R.
      • Rocha D.
      • Brochard M.
      • Boussaha M.
      • Boichard D.
      Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle.
      ) whose effects were later confirmed (
      • Sanchez M.P.
      • Wolf V.
      • El Jabri M.
      • Beuvier E.
      • Rolet-Répécaud O.
      • Gaüzère Y.
      • Minéry S.
      • Brochard M.
      • Michenet A.
      • Taussat S.
      • Barbat-Leterrier A.
      • Delacroix-Buchet A.
      • Laithier C.
      • Fritz S.
      • Boichard D.
      Short communication: Confirmation of candidate causative variants on milk composition and cheesemaking properties in Montbéliarde cows.
      ). Missing genotypes were imputed using FImpute software (
      • Sargolzaei M.
      • Chesnais J.
      • Schenkel F.
      A new approach for efficient genotype imputation using information from relatives.
      ) for the 47,794 autosomal SNPs [41,492 (50K) and 6,852 custom SNPs (50K+)] that passed all quality control filters (individual call rate >95%; SNP call rate >90%; minor allele frequency >1%; genotype frequencies in Hardy-Weinberg equilibrium with P > 10−4).
      The following traits were investigated in this study:
      • 3 laboratory cheese yields: fresh curd yield (CYFRESH) = 100 × (g of curd/g of milk), curd yield in dry matter (CYDM) = 100 × (1 – g of DM whey/g of DM milk), and curd yield in protein and fat (CYFAT-PROT) = (PC + FC) × (g of milk/g of curd), where PC = protein content and FC = fat content;
      • 5 coagulation traits for pressed cooked cheese (PCC) and soft cheese (SC): curd organization index (K10/RCTPCC and K10/RCTSC), where RCT is rennet coagulation time to 0.5 firm index (FI) and K10 is time to obtain 10 FI, curd firmness at RCT (aPCC and aSC), and at 2 times RCT (a2SC); and
      • total casein (ΣCN) and Ca contents in milk.
      The MIR prediction accuracy, as estimated by the value of R2 in a validation population, varied between 0.54 (CYFAT-PROT) and 0.98 (ΣCN) depending on the trait (Table 1). All CMP and milk composition traits were predicted from MIR spectra using equations developed in 3 different projects: From'MIR for CMP traits (
      • El Jabri M.
      • Sanchez M.P.
      • Trossat P.
      • Laithier C.
      • Wolf V.
      • Grosperrin P.
      • Beuvier E.
      • Rolet-Répécaud O.
      • Gavoye S.
      • Gauzere Y.
      • Belysheva O.
      • Notz E.
      • Boichard D.
      • Delacroix-Buchet A.
      Comparison of Bayesian and partial least squares regression methods for mid-infrared prediction of cheese-making properties in Montbéliarde cows.
      ), PhénoFinLait for caseins (
      • Ferrand M.
      • Miranda G.
      • Guisnel S.
      • Larroque H.
      • Leray O.
      • Lahalle F.
      • Brochard M.
      • Martin P.
      Determination of protein composition in milk by mid-infrared spectrometry.
      ), and Optimir for Ca (
      • Soyeurt H.
      • Bruwier D.
      • Romnee J.
      • Gengler N.
      • Bertozzi C.
      • Veselko D.
      • Dardenne P.
      Potential estimation of major mineral contents in cow milk using mid-infrared spectrometry.
      ;
      • Gengler N.
      • Soyeurt H.
      • Dehareng F.
      • Bastin C.
      • Colinet F.
      • Hammami H.
      • Vanrobays M.-L.
      • Lainé A.
      • Vanderick S.
      • Grelet C.
      • Vanlierde A.
      • Froidmont E.
      • Dardenne P.
      Capitalizing on fine milk composition for breeding and management of dairy cows.
      ).
      Table 1Features of cheese-making properties and milk composition traits
      Trait
      For pressed cooked cheese (PCC) and soft cheese (SC).
      Description and unitMeanSDR2 val
      Accuracy of mid-infrared prediction equations (R2 val).
      h2 TD model
      Heritability (h2) and genetic standard deviation (σG) estimates for the test-day (TD; Sanchez et al., 2018a) and lactation (LACT) models.
      h2 LACT model
      Heritability (h2) and genetic standard deviation (σG) estimates for the test-day (TD; Sanchez et al., 2018a) and lactation (LACT) models.
      σG
      Heritability (h2) and genetic standard deviation (σG) estimates for the test-day (TD; Sanchez et al., 2018a) and lactation (LACT) models.
      CYFRESH100 × (g of curd/g of milk), %37.67.50.820.380.703.2
      CYDM100 × (g DM of curd/g DM of milk), %66.95.00.890.390.722.2
      CYFAT-PROT(g of milk fat + g of milk protein)/kg of curd, g/kg187.621.70.540.370.708.2
      aPCCCurd firmness at rennet coagulation time (RCT), firm index (FI)19.12.60.760.420.721.2
      K10/RCTPCCTime to obtain 10 FI from RCT (K10)/RCT, min/min0.370.100.680.470.710.044
      aSCCurd firmness at RCT, FI19.42.70.760.450.731.3
      a2SCCurd firmness at 2 times RCT, FI23.22.10.690.480.721.0
      K10/RCTSCTime to obtain 10 FI from RCT (K10)/RCT, min/min0.360.110.720.470.720.049
      ΣCNTotal caseins, g/100 g of protein80.91.30.980.460.730.12
      CaCalcium, mg/kg of milk1,161950.820.500.7552.6
      1 For pressed cooked cheese (PCC) and soft cheese (SC).
      2 Accuracy of mid-infrared prediction equations (R2 val).
      3 Heritability (h2) and genetic standard deviation (σG) estimates for the test-day (TD;
      • Sanchez M.P.
      • El Jabri M.
      • Minéry S.
      • Wolf V.
      • Beuvier E.
      • Laithier C.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Genetic parameters for cheese-making properties and milk composition predicted from mid-infrared spectra in a large dataset of Montbéliarde cows.
      ) and lactation (LACT) models.

      Training and Validation Data Sets

      To estimate the accuracy of genomic predictions, the data set was split into 2, so that training and validation sets were as independent as possible (Figure 1). Primiparous cows with genotypes were included in the training population (Training) if they had calved before May 2016 and if they had at least 3 TD records (13,272 cows). The validation population (VAL) was composed of genotyped primiparous cows who had calved after August 2016 and had at least 5 TD records (4,032 cows). All performance data of the Training cows that had been recorded after August 2016 were removed (Figure 1a). Out of the 642 bulls with genotyped daughters in the whole data set, 180 had daughters in both populations. Including these families in the analysis would result in increased accuracy but also increased bias. Therefore, we considered different bull families in each data set: we kept cows (genotyped or not) in the training set if at least 40% of the daughters of a bull were Training cows, and we removed all half-sisters of these cows in the VAL set. For all other bulls (with less than 40% of daughters in Training), we kept the genotyped daughters in the VAL set and removed all half-sisters (genotyped or not) of these cows from the Training set. We selected 155,961 cows (12,850 with genotypes and 143,111 without genotypes) from 3,278 bulls (including 496 bulls with genotyped daughters) for use in Training and 2,125 genotyped cows from 146 other bulls for use in the VAL set (Figure 1b).
      Figure thumbnail gr1
      Figure 1Selection of cows for training and validation sets.

      Single-Step Genomic Evaluation

      Single-step estimated breeding values (ssEBVs) were calculated using the hybrid single-step model proposed by
      • Fernando R.
      • Cheng H.
      • Golden B.
      • Garrick D.
      Computational strategies for alternative single-step Bayesian regression models with large numbers of genotyped and non-genotyped animals.
      and implemented in HSSGBLUP software (

      Tribout, T., V. Ducrocq, and D. Boichard. 2020. HSSGBLUP: A single-step SNP BLUP genomic evaluation software adapted to large livestock populations. In ICQG, Brisbane, Australia.

      ). Briefly, this method combines information on phenotypes, pedigrees, and genotypes to predict genomic breeding values for all the animals of the pedigree (genotyped or not) and directly provides marker effect estimates. Coherence between pedigree and genotypes was obtained by fitting a J vector (
      • Hsu W.L.
      • Garrick D.J.
      • Fernando R.L.
      The accuracy and bias of single-step genomic prediction for populations under selection.
      ) for each unknown parent group.

      Models and Scenarios

      The ssEBVs were estimated using 2 categories of models that respectively considered the mean performance over the course of a lactation (LACT) or the individual test-day records (TD) of each cow. Both LACT and TD models were applied to the first lactation (L1) as well as to the first 3 lactations (L3). In total, 4 different models were tested: LACT1 (1 observation per cow in L1), LACT3 (up to 3 observations per cow, 1.8 on average in L1–L3), TD1 (with an average of 7.1 observations per cow in L1), and TD3 (with an average of 13.3 observations per cow in L1–L3). The effects and number of observations included in each model are described in Table 2. The contemporary group was the herd × year combination for LACT and the herd × test-day for TD. The other fixed effects were the year × month of calving and the age of the cow at first calving. In the TD models, the lactation curve was fitted by the effect of DIM (one class every 10 d, within parity). In the LACT models, the effect of the number of TD records was included as a proxy of average DIM. Variances of the random effects of each model were estimated using Wombat software (
      • Meyer K.
      WOMBAT - A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML).
      ) with a pedigree relationship matrix. More details can be found in
      • Sanchez M.P.
      • El Jabri M.
      • Minéry S.
      • Wolf V.
      • Beuvier E.
      • Laithier C.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Genetic parameters for cheese-making properties and milk composition predicted from mid-infrared spectra in a large dataset of Montbéliarde cows.
      .
      Table 2Features of models and training populations tested
      NameLACT1LACT3TD1TD3
      PhenotypeMean phenotype over lactation (LACT) 1Mean phenotype over lactations 1 to 3Test-day (TD) records in lactation 1Test-day records in lactations 1 to 3
      Model
      a, p, and e are animal, permanent environment, and residual random effects with variances σa2, σp2, and σe2 estimated with Wombat software (Meyer, 2007), respectively; X and Z are incidence matrices for fixed and random effects, respectively.
      y = Xb + Za + ey = Xb + Za + Zp + ey = Xb + Za + Zp + ey = Xb + Za + Zp + e
      Fixed effects bHerd × year Month × year calving Age at first calving No. of TD recordsHerd × year Month × year calving Age at calving No. of TD records × parityHerd × TD Month × year calving Age at first calving DIMHerd × TD Month × year calving Age at calving DIM × parity
      Training population
       No. of cows with phenotypes
      The number of cows with phenotypes differed slightly among models due to the filters applied to the corresponding data set.
      155,479155,919155,961155,961
       No. of TD records155,479287,8001,115,0002,070,000
       No. of animals in pedigree413,787422,857414,919423,348
       No. of animals in pedigree with genotypes21,76421,87421,74721,874
      1 a, p, and e are animal, permanent environment, and residual random effects with variances σa2, σp2, and σe2 estimated with Wombat software (
      • Meyer K.
      WOMBAT - A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML).
      ), respectively; X and Z are incidence matrices for fixed and random effects, respectively.
      2 The number of cows with phenotypes differed slightly among models due to the filters applied to the corresponding data set.
      In addition to the polygenic models without genotyping information (LACT1-BLUP and TD1-BLUP), the LACT1, LACT3, TD1, and TD3 models were first tested using the 50K scenario (LACT1–50K, LACT3–50K, TD1–50K, and TD3–50K, respectively), and included all polymorphic SNPs of the BovineSNP50 chip (i.e., 41,942 SNPs). From these, we selected the model that generated the most accurate ssEBVs (TD3) and tested 2 additional scenarios in which we added SNPs from the custom part of the EuroG10K chip (i.e., considering a total of 47,794 SNPs). In the TD3–50K+ scenario, we used constant effect variances for all 47,794 SNPs, whereas in the TD3-CAND scenario, we established specific variances for 7 to 14 candidate SNPs, depending on the trait (Table 3). These SNPs had been identified in previous studies as the best candidates for QTL with large effects on milk composition and CMP traits (
      • Sanchez M.P.
      • Govignon-Gion A.
      • Croiseau P.
      • Fritz S.
      • Hozé C.
      • Miranda G.
      • Martin P.
      • Barbat-Leterrier A.
      • Letaïef R.
      • Rocha D.
      • Brochard M.
      • Boussaha M.
      • Boichard D.
      Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle.
      ,
      • Sanchez M.P.
      • Wolf V.
      • El Jabri M.
      • Beuvier E.
      • Rolet-Répécaud O.
      • Gaüzère Y.
      • Minéry S.
      • Brochard M.
      • Michenet A.
      • Taussat S.
      • Barbat-Leterrier A.
      • Delacroix-Buchet A.
      • Laithier C.
      • Fritz S.
      • Boichard D.
      Short communication: Confirmation of candidate causative variants on milk composition and cheesemaking properties in Montbéliarde cows.
      ,
      • Sanchez M.P.
      • Ramayo-Caldas Y.
      • Wolf V.
      • Laithier C.
      • El Jabri M.
      • Michenet A.
      • Boussaha M.
      • Taussat S.
      • Fritz S.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows.
      ). In the TD3-CAND scenario, we calculated the proportion of genetic variance for each SNP and each trait as follows: %σa2=100(2p(1p)α2σa2), with σa2 being the additive genetic variance, p and (1 − p) being the allelic frequencies, and α being the allelic substitution effect estimated in GWAS. The cumulative effects of the candidate SNPs explained from 17.5% to 58.4% of the total additive genetic variance of the traits, with the remaining part (41.6% to 82.5%) being equally distributed among all other SNPs (47,780 to 47,787 depending on the trait).
      Table 3Candidate SNPs of the 18 QTL detected for milk composition and cheese-making traits
      MAF: minor allele frequency; UTR: untranslated region.
      QTLBTASNP IDMAF% of additive genetic variance associated with each trait × SNP combinationNo. of traits affectedFunctional annotation
      CYFRESH
      For each SNP and trait, the % of genetic variance was calculated as follows: %σa2=100(2p(1−p)α2σa2), with σa2 being the additive genetic variance, p and (1 − p), allelic frequencies, and α, the allelic substitution effect estimated from GWAS; see description of traits in Table 1.
      CYDMCYFAT-PROTaPCCK10/RCT PCCaSCa2SCK10/RCT SCΣCNCa
      11rs1329653710.450.751.10.483Intron SLC37A1
      22rs1101231050.790.550.940.580.514Intergenic
      32rs1339443760.370.391.32Missense ALPL
      45rs2112105690.121.31.81.13Intron MGST1
      55rs5258807460.042.82.71.23.90.413.52.61.95.49Upstream GRAMD4
      66rs1368007730.732.53.04.42.32.04.30.737Downstream SLC34A2
      76rs437030160.407.77.28.72.342.74.335.842.11.51.310Missense CSN3
      87rs290244140.081.53.02.91.02.23.16Upstream ENSBTAG00000004032
      911rs1100662290.468.310.423.010.42.13.02.83.316.29Missense PAEP
      1014rs1090355860.3312.716.42.83.43.22.32.54.28Upstream GPT
      1116rs1090330260.080.450.450.260.580.580.470.570.528Intergenic
      1217rs4485010710.061.91.21.00.584Upstream BRI3BP
      1319rs1360670460.321.41.72Upstream FASN
      1420rs3857448460.070.411.30.562.40.872.32.27Intron ANKH
      1522rs416424780.261.92.20.730.362.00.901.31.60.789Intron FAM19A4
      1624rs1094782900.210.460.480.550.620.410.486Intron LMAN1
      1726rs412556930.460.500.592Missense SCD
      1827rs2086752760.411.82.61.735′UTR GPAT4
      Total (% σa2)42.946.545.922.857.118.349.358.434.217.5
      No. of QTL1411119810912117
      1 MAF: minor allele frequency; UTR: untranslated region.
      2 For each SNP and trait, the % of genetic variance was calculated as follows: %σa2=100(2p(1p)α2σa2), with σa2 being the additive genetic variance, p and (1 − p), allelic frequencies, and α, the allelic substitution effect estimated from GWAS; see description of traits in Table 1.

      Reliability and Bias

      The reliability and bias of each model and scenario were then estimated using the validation population. The ssEBVs were estimated by applying the Training genomic prediction equations to the genotypes of the 2,125 validation cows. Using GENEKIT software (
      • Ducrocq V.
      Genekit, BLUP software. June 2011 version.
      ), we adjusted the phenotypes of primiparous VAL cows for nongenetic effects estimated with the TD1 model (described in Table 2) applied to the complete Training + VAL data set, but ignoring marker information. To compare the accuracies of the different models, we calculated, for each trait, the correlation coefficients (r) between means of individual performances (corrected for nongenetic effects using the TD1 model) and ssEBVs estimated with the LACT1–50K (r LACT1–50K), LACT3–50K (r LACT3–50K), TD1–50K (r TD1–50K), or TD3–50K (r TD3–50K) model. Similarly, for the TD3 model, the 3 scenarios were compared by calculating correlations between mean individual performances (corrected for nongenetic effects using the TD1 BLUP model) and ssEBVs estimated with the 50K (r TD3–50K), 50K+ (r TD3–50K+), or CAND (r TD3-CAND) scenario. Reliabilities were calculated by dividing the squared correlations (r2 LACT1–50K, r2 LACT3–50K, r2 TD1–50K, r2 TD3–50K, r2 TD3–50K+, or r2 TD3-CAND) by the heritability of the corresponding trait estimated using Wombat software (
      • Meyer K.
      WOMBAT - A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML).
      ) with the LACT1 model; these are reported in Table 1. Bias was estimated as the deviation from 1 of the slope of the regression of corrected performances on ssEBVs.

      Estimation of Genetic Trends

      To estimate genetic trends in CMP and milk composition traits in Montbéliarde cattle, we applied the TD3–50K model to the whole From'MIR data set. This was equivalent to the data set used when the TD3–50K model was applied to the training population (423,348 animals in the pedigree, of which 21,874 with genotypes for 41,942 SNPs), except that all phenotypes were kept (i.e., 2.6 million TD records from 190,261 cows).
      We then applied the genomic prediction equations obtained from the whole data set to all genotyped Montbéliarde animals belonging to the Umotest breeding company, including 17,859 males born between 2005 and 2018 and 256,861 females born between 2008 and 2018. Genetic trends were assessed separately for males and females by averaging ssEBVs per year of birth.

      Breeding Simulation with Different Breeding Goals

      Finally, we conducted a breeding simulation to investigate the effect of incorporating cheese-making traits into the breeding goal of the population. We compared 3 different breeding goals: the total merit index (TMI), which is the current breeding goal in French Montbéliarde cattle, and 2 alternative indexes (TMI-COMP and TMI-Cheese), defined as follows:
      • 1)
        TMI = 0.45 Milk + 0.145 Udder-Health + 0.18 REPROD + 0.05 LONG + 0.05 Speed + 0.125 MORPH, with weights expressed in genetic standard deviations. This is a combination of
        • the production index Milk = 14 PY + 2 FY + 4.2 PC + 1.3 FC, which combines breeding values for milk protein yield (PY), fat yield (FY), protein content (PC), and fat content (FC);
        • the udder health index Udder-Health = 0.6 Cell + 0.4 MAST, which combines breeding values for somatic cell score (Cell) and clinical mastitis (MAST);
        • the fertility index REPROD = 0.5 FERC + 0.25 FERH + 0.25 CALV-AI, which combines breeding values for conception rate of cows (FERC) and heifers (FERH) and the interval between calving and first artificial insemination (CALV-AI);
        • functional longevity (LONG);
        • speed of milking (Speed);
        • the overall morphology index MORPH = 0.4 Udder + 0.3 Body + 0.15 Legs + 0.10 Rump + 0.05 Muscularity, with Udder, Body, Legs, Rump, and Muscularity being the aggregate index for udder conformation, body capacity, feet and legs, rump morphology, and muscularity, respectively;
      • 2)
        TMI-COMP = 0.7 TMI + 0.3 ΣCN;
      • 3)
        TMI-Cheese = 0.7 TMI + 0.1 (CYDM + aPCC − K10/RCTPCC).
      We simulated selection by truncation on bulls based on the breeding strategy applied by the Umotest breeding company. For the traits investigated in the present study, we considered the ssEBVs of bulls estimated with the TD3–50K model with the whole data set, whereas for all other traits routinely included in genomic evaluations, we considered official genomic indexes. To ensure a sufficient number of bulls per year of birth, only the 14,389 bulls born between 2009 and 2017 were considered (1,022 to 2,210 bulls per year of birth). Within each year of birth, we sorted bulls by decreasing TMI, TMI-COMP, or TMI-Cheese value, calculated using the formulas described above. In a given year, we selected the 80 bulls with the highest TMI, TMI-COMP, or TMI-Cheese indexes among all the candidates (1,600 candidates on average). We then calculated the responses to selection for the ssEBVs estimated in the present study (10 milk CMP and composition traits) as well as for the 8 official indexes described earlier: TMI, Udder-Health, REPROD, PY, FY, PC, FC, and milk yield (MY).
      For each breeding goal scenario, we estimated the annual responses to selection by calculating the selection differential of each ssEBV or index (i.e., the difference between the mean of the 80 best bulls and the mean of all candidates). The yearly responses, expressed in genetic standard deviations, were then averaged over the whole period (2009–2017).

      RESULTS

      Reliability of Single-Step Genomic Evaluation

      Estimates of heritability coefficients and genetic standard deviations are in Table 1. The heritability coefficients for the lactation traits were used to estimate the reliability of the evaluations.
      To estimate the reliability of single-step genomic evaluation for milk cheese-making traits, we removed these phenotypes from the validation population and calculated breeding values (ssEBVs). We then assessed the reliabilities of each model × trait combination in the validation population.
      As described in the Materials and Methods, we first compared 6 models, without genotyping data (LACT1-BLUP and TD1-BLUP) or applied to the 50k genotypes (LACT1–50K, LACT3–50K, TD1–50K, and TD3–50K). The correlations, calculated between ssEBVs and adjusted performance, are presented in Supplemental Table S1 (https://doi.org/10.6084/m9.figshare.19477925). All genomic models largely outperformed the polygenic model. Mean correlation values of all traits for LACT1 and TD1 models were, respectively, 0.37 and 0.38 for polygenic models and 0.66 and 0.67 for genomic models. In the following part, only genomic models will be compared.
      Regardless of the model tested and the trait analyzed, all reliabilities were higher than 59% (Figure 2). We found slightly higher reliabilities for milk coagulation and composition traits, but overall, reliabilities were similar for all traits. For a given category of model, based on performance data averaged over the course of lactation (LACT) or individual test-day records (TD), reliabilities were slightly higher for models that included data from the second and third lactations together with those from the first (+1 point on average). In all cases, the TD models (TD1–50K and TD3–50K) gave more reliable ssEBVs than the LACT models (LACT1–50K and LACT3–50K) by about 2 points on average. Across all traits, average reliabilities were 61%, 62%, 63%, and 64% for the LACT1–50K, LACT3–50K, TD1–50K, and TD3–50K models, respectively.
      Figure thumbnail gr2
      Figure 2Reliabilities estimated in the validation population (n = 2,150) with the 4 tested models. LACT = lactation; TD = test-day. CYFRESH = fresh curd yield; CYDM = curd yield in dry matter; CYFAT-PROT = curd yield in protein and fat; K10/RCTPCC curd organization index for pressed cooked cheese, where RCT is rennet coagulation time to 0.5 firm index (FI) and K10 is time to obtain 10 FI; aPCC = curd firmness at RCT for pressed cooked cheese; K10/RCTSC = curd organization index for soft cheese; aSC = curd firmness at RCT for soft cheese; a2SC = curd firmness at 2 times RCT; ΣCN = total casein; Ca = Ca content.
      Because the TD3 model gave the most reliable results for all traits, we applied this model to the training population to test 3 scenarios that differed in the density or weighting (or both) of SNPs (TD3–50K, TD3–50K+, and TD3-CAND). We then compared the reliabilities and biases of ssEBVS calculated with the 3 scenarios in the validation population. The inclusion of SNPs from the custom part of the EuroG10K chip (the TD3–50K+ scenario: 5,852 additional SNPs compared with the TD3–50K scenario) further improved the reliability of ssEBVs (Figure 3). Across all traits, the average gain was about 2 points: from +2.7 to +3.4 points for cheese yields and from +1.1 to +1.9 points for coagulation parameters and casein and calcium content. The TD3-CAND scenario, which assigned higher weights to candidate causal SNPs, provided only a marginal gain of 0.25 points on average, with differences depending on the trait (from +0.2 to +1.1 points for aPCC, ΣCN, aSC, a2SC, CYDM, CYFRESH, CYFAT-PROT, and Ca versus from −0.1 and −0.7 points for K10/RCTSC and K10/RCTPCC, respectively). Instead, when we evaluated bias in the scenarios, based on the slope of the regression line of corrected performance on ssEBVs, we found different results (Figure 4). For all traits and all 3 scenarios, slopes were less than 1, indicating bias in the estimation of ssEBVs. The TD3-CAND scenario generated the least biased ssEBVs, whereas the other 2 scenarios gave equivalent results, with the TD3–50K scenario slightly outperforming the TD3–50K+ scenario.
      Figure thumbnail gr3
      Figure 3Reliabilities estimated in the validation population (n = 2,150) with the 3 tested scenarios. TD = test day; CAND = weighted candidate SNPs.
      Figure thumbnail gr4
      Figure 4Biases estimated in the validation population (n = 2,150) with the 3 tested scenarios. TD = test day; CAND = weighted candidate SNPs.

      Genetic Trends in Cheese-Making Traits

      Genetic trends in CMP traits were estimated by calculating the mean genetic values (ssEBVs) of males and females per year of birth, and are shown in Figure 5. Over the period 2005–2018, we observed regular improvement in ssEBVs for all of the cheese-making criteria and milk composition traits analyzed in this study (i.e., an increase in ssEBVs for all traits except for 2 traits for which reduction is desirable: CYFAT-PROT, which is inversely proportional to curd weight, and K10/RCT, which is inversely proportional to the curd organization speed). In the Comté PDO region, the average genetic level of Montbéliarde bulls and cows with respect to cheese yields and coagulation parameters is therefore higher today than in 2005. When expressed as genetic standard deviations (σg), the traits with the most marked genetic improvement over this period were coagulation parameters and total caseins in milk (between 0.5 and 0.6 σg). The genetic trend was instead more moderate for cheese yields (between 0.2 and 0.4 σg) and calcium content in milk (0.1 σg).
      Figure thumbnail gr5
      Figure 5Genetic trends in genetic standard deviation (σa) estimated for bulls born between 2005 and 2018 and cows born between 2008 and 2018 for milk cheese-making and composition traits. CYFRESH = fresh curd yield; CYDM = curd yield in dry matter; CYFAT-PROT = curd yield in protein and fat; K10/RCTPCC curd organization index for pressed cooked cheese, where RCT is rennet coagulation time to 0.5 firm index (FI) and K10 is time to obtain 10 FI; aPCC = curd firmness at RCT for pressed cooked cheese; K10/RCTSC = curd organization index for soft cheese; aSC = curd firmness at RCT for soft cheese; a2SC = curd firmness at 2 times RCT; ΣCN = total casein; Ca = Ca content.

      Effect of Different Breeding Strategies on Cheese-Making Traits and Other Traits

      Annual responses to selection, averaged over the 2009–2017 period, were calculated with the current breeding objective (TMI) and 2 alternative scenarios that gave a weight of 70% to TMI and 30% to either casein content in milk (TMI-COMP) or to a combination of 3 cheese-making traits (TMI-Cheese). As expected, the addition of new criteria to the breeding goal had an effect on the breeding values of all investigated traits (Figure 6). For almost all traits, the TMI-COMP scenario yielded values that were intermediate between those of the TMI and TMI-Cheese scenarios. When we examined the traits incorporated in the current breeding goal (i.e., noncheese-making traits), we found that both alternative scenarios resulted in lower responses than TMI for PY, FY, MY, and Udder-Health, a similar response for REPROD, and higher responses for PC and FC. However, the reductions in gain per year were relatively limited (0.1 σg for PY, FY, and Udder-Health and 0.2 σg for Milk) considering the high weight given to casein content or cheese-making traits in the alternative scenarios. Instead, the additional emphasis on selection on cheese-making traits in bulls in TMI-COMP, and a fortiori in TMI-Cheese, led to much higher response levels for all cheese-making traits included in this study (i.e., all cheese yields and coagulation parameters measured on pressed cooked and SC). For CYDM, for example, the annual gain was almost doubled with the TMI-COMP scenario (+0.68 σg) and almost tripled with the TMI-Cheese scenario (+0.95 σg) compared with the TMI scenario (+0.38 σg). Gains for other cheese-making traits were similar. Milk composition traits were also significantly improved, to a similar extent, in both alternative scenarios. Total casein content increased annually by 0.52, 0.86, and 0.85 σg, whereas calcium content increased by 0.38, 0.50, and 0.52 σg with the TMI, TMI-COMP, and TMI-Cheese scenarios, respectively.
      Figure thumbnail gr6
      Figure 6Responses to selection for actual total merit index without cheese-making properties (TMI), 70% TMI + 30% casein content (TMI-COMP), and 70% TMI + 30% cheese-making properties (TMI-Cheese) on milk cheese-making and composition traits and other traits (in genetic standard deviations). CYFRESH = fresh curd yield; CYDM = curd yield in dry matter; CYFAT-PROT = curd yield in protein and fat; K10/RCTPCC curd organization index for pressed cooked cheese, where RCT is rennet coagulation time to 0.5 firm index (FI) and K10 is time to obtain 10 FI; aPCC = curd firmness at RCT for pressed cooked cheese; K10/RCTSC = curd organization index for soft cheese; aSC = curd firmness at RCT for soft cheese; a2SC = curd firmness at 2 times RCT. PC = protein content; FC = fat content; FY = fat yield; PY = protein yield; MY = milk yield; REPROD = reproduction.

      DISCUSSION

      The results presented in this study are very encouraging for the implementation of single-step genomic evaluation of CMP in the Montbéliarde breed. Breeding values with high reliability were obtained and favorable genetic trends in CMP traits were identified over the last 13 yr. However, we showed that amending the current breeding goal to include a combination of 3 CMPs could result in more rapid improvement to all CMP traits, with only a limited effect on other traits currently selected in the Montbéliarde breed.
      By testing different models and scenarios, we showed that all genomic models strongly outperformed the polygenic model (with a reliability over 0.59 vs. less than 0.25) and we were able to identify which genomic model gave the best results. The model that generated the most reliable breeding values was the test-day model that considered all the observations recorded once a month during the first 3 lactations of the cows. Compared with the other models tested (lactation models with performance data averaged over the lactation period or the test-day model for the first lactation only), the best-performing model integrated the largest amount of information and was more accurate in modeling phenotypes. Furthermore, by incorporating SNPs from the custom part of the EuroG10K BeadChip, we included those that were particularly concentrated in the genomic regions affecting many important traits (QTL) in cattle, which probably explained the gain in reliability obtained with these scenarios (50K+ and CAND scenarios) compared with the scenario based only on 50K SNPs. Indeed, in addition to having the candidate variants identified by the PhenoFinlait and From'MIR projects, the custom part of the EuroG10K chip is also particularly enriched in SNPs from QTL regions associated with milk production traits (
      • Boichard D.
      • Boussaha M.
      • Capitan A.
      • Rocha D.
      • Hozé C.
      • Sanchez M.P.
      • Tribout T.
      • Letaief R.
      • Croiseau P.
      • Grohs C.
      • Li W.
      • Harland C.
      • Charlier C.
      • Lund M.S.
      • Sahana G.
      • Georges M.
      • Barbier S.
      • Coppieters W.
      • Fritz S.
      • Guldbrandtsen B.
      Experience from large scale use of the EuroGenomics custom SNP chip in cattle. In 11th WCGALP. Auckland, New Zealand.
      ). Among these QTL regions are, for example, the regions of the casein or PAEP genes, which have been widely investigated in recent decades for their strong effects on milk composition and cheese-making traits in many dairy cattle breeds (
      • Martin P.
      • Szymanowska M.
      • Zwierzchowski L.
      • Leroux C.
      The impact of genetic polymorphisms on the protein composition of ruminant milks.
      ;
      • Ganai N.A.
      • Bovenhuis H.
      • van Arendonk J.A.
      • Visker M.H.
      Novel polymorphisms in the bovine beta-lactoglobulin gene and their effects on beta-lactoglobulin protein concentration in milk.
      ). Scenarios with custom SNPs are mostly beneficial to cheese yields, which are highly genetically correlated with milk fat content (
      • Sanchez M.P.
      • El Jabri M.
      • Minéry S.
      • Wolf V.
      • Beuvier E.
      • Laithier C.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Genetic parameters for cheese-making properties and milk composition predicted from mid-infrared spectra in a large dataset of Montbéliarde cows.
      ) and thus particularly affected by the DGAT1 gene region (
      • Sanchez M.P.
      • Ramayo-Caldas Y.
      • Wolf V.
      • Laithier C.
      • El Jabri M.
      • Michenet A.
      • Boussaha M.
      • Taussat S.
      • Fritz S.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows.
      ). This region contains numerous SNPs that are featured in the research add-on of the EuroG10K chip (153 SNPs in the 50K+ scenario versus 25 SNPs in the 50K scenario). As a consequence, the average absolute values of SNP effects estimated in this region were higher in the 50K scenario than in the 50K+ or CAND scenarios (e.g., 0.052, 0.021, and 0.019 points for cheese yield in DM, respectively); this was due to the fact that, in the 50K+ and CAND scenarios, the effect of the chromosomal region was distributed over a larger number of SNPs, each with smaller individual effects. As expected, in the CAND scenario, the candidate variant with 16% of the genetic variance of the trait had the largest estimated effect (0.75 points), whereas the other 152 SNPs in the region had a very small estimated effect. Regardless of the weights attributed to the SNPs, the relative SNP enrichment of this region in the 50K+ and CAND scenarios, and the consequent increase in LD with causal mutations, probably helped to better capture the effects of the causal mutations and thus increase the reliabilities of breeding values. The inclusion of weighted candidate mutations (CAND scenario) resulted in a somewhat similar result in terms of reliability, but with a different distribution of effects: the strong effect of the candidate causal mutation led to reductions in the effects of nearby SNPs. However, the slight increase in reliability in the CAND scenario was also accompanied by a reduction in bias of genetic values. This was likely a reflection of the fact that the estimated effect of a causal variant is expected to be more stable over generations than the estimated effects of neutral markers, which depend on recombinations and distance to the reference population and are thus more prone to inflation.
      The reliabilities of the breeding values we obtained (66% on average for all traits with the best model) were high in spite of a medium-size reference population (3,278 sires with a total of 155,961 daughters having MIR information, including 12,850 cows with genotypes in the training population). This result is explained by both the heritability of the traits and the number of repeated records per cow, resulting in highly informative phenotypes per cow, explain these reliabilities. They are consistent with the theoretical accuracy expected given the heritability estimates of these traits and the size of the reference population included in the present study (
      • Goddard M.
      • Hayes B.
      Mapping genes for complex traits in domestic animals and their use in breeding programmes.
      ). With this level of reliability, comparable to that obtained for other currently selected traits (e.g., milk production traits), it becomes possible to implement a program of genomic evaluation of cheese-making traits predicted from MIR spectra. It is important to note here that our division of the reference population into 2 independent subpopulations, although necessary to give rigorous reliability estimates, had the effect of greatly reducing the size of the population used to estimate SNP effects (12,850 cows out of the 19,564 cows in the reference population). The inclusion of all cows with phenotypes and genotypes, including the new ones generated since the time of the study (around +40,000/yr), will further improve the prediction accuracy of genomic values.
      Another point to consider is that CMP measurements are indirect MIR predictions and all our results (reliability but also genetic trend; next paragraph) assume that these predictions are robust over the whole data set (i.e., over herds and years, beyond fat and protein contents). Even if it is probably the case for the best-predicted traits, our reliability results may be optimistic and a periodic validation of the CMP prediction equations will be useful to avoid any potential drift between real CMP values and predictions.
      The favorable genetic trends for cheese yields and coagulation parameters we observed over the last 13 yr are likely due to strong genetic correlations between milk protein and fat content and cheese-making traits (
      • Sanchez M.P.
      • El Jabri M.
      • Minéry S.
      • Wolf V.
      • Beuvier E.
      • Laithier C.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Genetic parameters for cheese-making properties and milk composition predicted from mid-infrared spectra in a large dataset of Montbéliarde cows.
      ), the latter having probably been selected indirectly together with the former, which is included in the TMI in the Montbéliarde breed. We found a particularly strong response for total casein content, which represents about 80% of total proteins in milk. This result suggests that selection on protein content, which has a relatively important economic weight in the Montbéliarde TMI, has led to an improvement of CMP and, in particular, coagulation parameters, which present the strongest genetic correlations with protein content (
      • Sanchez M.P.
      • El Jabri M.
      • Minéry S.
      • Wolf V.
      • Beuvier E.
      • Laithier C.
      • Delacroix-Buchet A.
      • Brochard M.
      • Boichard D.
      Genetic parameters for cheese-making properties and milk composition predicted from mid-infrared spectra in a large dataset of Montbéliarde cows.
      ). These results are consistent with the increase detected in the frequencies of some milk protein variants that are associated with better cheese-making abilities in French dairy cattle breeds (
      • Sanchez M.P.
      • Fritz S.
      • Patry C.
      • Delacroix-Buchet A.
      • Boichard D.
      Frequencies of milk protein variants and haplotypes estimated from genotypes of more than 1 million bulls and cows of 12 French cattle breeds.
      ).
      To enhance this positive trend, we show here that a more significant genetic gain in cheese-making traits could be possible with only a limited effect on the other traits currently considered in the breeding objective for Montbéliarde cattle. This could be accomplished by amending the breeding objective to include 3 traits (cheese yield in DM and 2 coagulation traits for PCC), which were identified as the most relevant by a group of experts, composed of farmers, cheesemakers, and ripeners in the Franche-Comté region.
      Despite the potential economic benefit of improving CMP, to the best of our knowledge, very few countries currently include CMP in their breeding objective. The United States calculates a so-called cheese index (Dollars Cheese Merit Index, CM$), but like the current TMI in Montbéliarde, it simply gives a greater weight to protein content (
      • VanRaden P.M.
      Invited review: Selection on net merit to improve lifetime profit.
      ). In Belgium, the ProFARMilk project (2011–2017) investigated the technological properties of milk, as predicted by MIR spectrometry, with respect to its ability to be processed into cheese, yogurt, and butter (

      Colinet, F. G., T. Troch, O. Abbas, V. Baeten, F. Dehareng, E. Froidmont, H. Soyeurt, P. Dardenne, M. Sindic, and N. Gengler. 2013. Potentiel d'utilisation de la spectrométrie moyen infrarouge pour prédire le rendement fromager du lait et étudier sa variabilité génétique. In Renc. Rech. Rum. Vol. 20, Paris, France.

      ). This project led to the implementation of a tool for monitoring milk processing abilities in the framework of milk recording in Wallonia, with the initial goal of discarding noncoagulating milk and eventually enabling the selection of these criteria (
      • AWE
      La transformation fermière du lait.
      ). Italy, which processes more than 80% of its milk into cheese, has for many years used indicators of cheese aptitude (e.g., casein content or lactodynamic behavior of milk measured with a Formagraph) in the milk payment scale in the Parmigiano-Reggiano PDO region (
      • Malacarne M.
      • Formaggioni P.
      • Franceschi P.
      • Summer A.
      Seasonal variations of milk quality in Parmigiano-Reggiano cheese manufacture on a period of 10 years.
      ). In addition, some coagulation parameters predicted by MIR spectrometry have been routinely recorded since 2011 (

      De Marchi, M., M. Penasa, F. Tiezzi, V. Toffanin, and M. Cassandro. 2012. Prediction of milk coagulation properties by Fourier Transform Mid-Infrared Spectroscopy (FTMIR) for genetic purposes, herd management and dairy profitability. Pages 47–53 in International Strategies and New Developments in Milk Analysis. VI ICAR Reference Laboratory Network Meeting. Vol. ICAR Technical Series No. 16, Cork, Ireland.

      ;
      • Pretto D.
      • Lopez-Villalobos N.
      • Penasa M.
      • Cassandro M.
      Genetic response for milk production traits, somatic cell score, acidity and coagulation properties in Italian Holstein-Friesian population under current and alternative selection indices and breeding objectives.
      ). Finally, in the Veneto region, an index (Cheese Aptitude Index) that combines coagulation speed (RCT) and firmness (a30), with equivalent weights, has been published periodically since January 2012 for Holstein bulls (http://www.intermizoo.com/research/cheese-aptitude-index).
      In France, traits associated with technological properties for cheese-making are not included in the current breeding objective of any dairy cattle breed. This study demonstrates, though, that their inclusion could have an obvious economic benefit for the dairy sector (e.g., by increasing cheese yields). However, it is worth noting that we have attributed a relatively high weight to cheese-making traits in the scenarios tested in the present study, to emphasize and illustrate the potential effect of selection. In practice, lower weight is given to CMP, ensuring positive gains for both MY and CMP. Selection on cheese-making traits will require further study to determine the optimal economic weights for these traits with respect to all the traits of the breeding goal. In addition, before selecting for cheese yields and milk coagulation traits, it will be necessary to assess the effects of such selection on the sensory quality of ripened cheeses. For this purpose, we have implemented a pilot genomic evaluation program in the Montbéliarde breed that makes indexes available for these traits, which will make it possible to control genetic trends in CMP with respect to the sensory quality of cheeses.

      CONCLUSIONS

      In this study, we investigated the potential for selecting milk technological traits, as predicted from MIR spectra, in the Montbéliarde dairy cattle breed. We showed it was possible to obtain reliable breeding values for cheese yields and coagulation parameters and identified the model that produced the most reliable and least biased results. The application of genomic prediction equations to all genotyped bulls and cows revealed favorable genetic trends in cheese-making traits over the last 13 yr. Finally, we found that incorporating cheese-making traits into the breeding goal could lead to more rapid improvement in these traits, with only a limited effect on other traits currently selected in this breed. Further investigation is needed to determine how best to integrate cheese-making traits into the breeding objective and to study the effect of selecting these traits on the quality of ripened cheeses. Regardless, all of these results are very encouraging for the implementation of single-step genomic evaluation of cheese-making parameters in the Montbéliarde breed, which could subsequently be extended to other dairy cattle breeds.

      ACKNOWLEDGMENTS

      This study was funded by the French Ministry of Agriculture, Agro-food, and Forests; the French Dairy Interbranch Organization (CNIEL); the Regional Union of Protected Designation cheeses of Franche-Comté (URFAC); and the Regional Council of Bourgogne-Franche-Comté, as part of the project FROM'MIR. The authors have not stated any conflicts of interest.

      REFERENCES

        • AWE
        La transformation fermière du lait.
        • Bittante G.
        • Penasa M.
        • Cecchinato A.
        Invited review: Genetics and modeling of milk coagulation properties.
        J. Dairy Sci. 2012; 95 (23021752): 6843-6870
        • Boichard D.
        • Boussaha M.
        • Capitan A.
        • Rocha D.
        • Hozé C.
        • Sanchez M.P.
        • Tribout T.
        • Letaief R.
        • Croiseau P.
        • Grohs C.
        • Li W.
        • Harland C.
        • Charlier C.
        • Lund M.S.
        • Sahana G.
        • Georges M.
        • Barbier S.
        • Coppieters W.
        • Fritz S.
        • Guldbrandtsen B.
        Experience from large scale use of the EuroGenomics custom SNP chip in cattle. In 11th WCGALP. Auckland, New Zealand.
        Molecular Genetics. 2018; 4: 675
        • Boichard D.
        • Govignon-Gion A.
        • Larroque H.
        • Maroteau C.
        • Palhiere I.
        • Tosser-Klop G.
        • Rupp R.
        • Sanchez M.P.
        • Brochard M.
        • Amigues Y.
        • Boscher M.Y.
        • Leveziel H.
        Genetic determinism of milk composition in fatty acids and proteins in ruminants, and selection potential.
        INRA Prod. Anim. 2014; 27: 283-298
      1. Colinet, F. G., T. Troch, O. Abbas, V. Baeten, F. Dehareng, E. Froidmont, H. Soyeurt, P. Dardenne, M. Sindic, and N. Gengler. 2013. Potentiel d'utilisation de la spectrométrie moyen infrarouge pour prédire le rendement fromager du lait et étudier sa variabilité génétique. In Renc. Rech. Rum. Vol. 20, Paris, France.

      2. De Marchi, M., M. Penasa, F. Tiezzi, V. Toffanin, and M. Cassandro. 2012. Prediction of milk coagulation properties by Fourier Transform Mid-Infrared Spectroscopy (FTMIR) for genetic purposes, herd management and dairy profitability. Pages 47–53 in International Strategies and New Developments in Milk Analysis. VI ICAR Reference Laboratory Network Meeting. Vol. ICAR Technical Series No. 16, Cork, Ireland.

        • De Marchi M.
        • Toffanin V.
        • Cassandro M.
        • Penasa M.
        Invited review: Mid-infrared spectroscopy as phenotyping tool for milk traits.
        J. Dairy Sci. 2014; 97 (24440251): 1171-1186
        • Ducrocq V.
        Genekit, BLUP software. June 2011 version.
        INRA GABI, 1998
        • El Jabri M.
        • Sanchez M.P.
        • Trossat P.
        • Laithier C.
        • Wolf V.
        • Grosperrin P.
        • Beuvier E.
        • Rolet-Répécaud O.
        • Gavoye S.
        • Gauzere Y.
        • Belysheva O.
        • Notz E.
        • Boichard D.
        • Delacroix-Buchet A.
        Comparison of Bayesian and partial least squares regression methods for mid-infrared prediction of cheese-making properties in Montbéliarde cows.
        J. Dairy Sci. 2019; 102 (31178172): 6943-6958
        • Fernando R.
        • Cheng H.
        • Golden B.
        • Garrick D.
        Computational strategies for alternative single-step Bayesian regression models with large numbers of genotyped and non-genotyped animals.
        Genet. Sel. Evol. 2016; 48 (27931187): 96
        • Ferrand M.
        • Miranda G.
        • Guisnel S.
        • Larroque H.
        • Leray O.
        • Lahalle F.
        • Brochard M.
        • Martin P.
        Determination of protein composition in milk by mid-infrared spectrometry.
        in: Proc. International Strategies and New Developments in Milk Analysis. VI ICAR Reference Laboratory Network Meeting, 2012: 41-45
        • Ganai N.A.
        • Bovenhuis H.
        • van Arendonk J.A.
        • Visker M.H.
        Novel polymorphisms in the bovine beta-lactoglobulin gene and their effects on beta-lactoglobulin protein concentration in milk.
        Anim. Genet. 2009; 40 (19032698): 127-133
        • Gengler N.
        • Soyeurt H.
        • Dehareng F.
        • Bastin C.
        • Colinet F.
        • Hammami H.
        • Vanrobays M.-L.
        • Lainé A.
        • Vanderick S.
        • Grelet C.
        • Vanlierde A.
        • Froidmont E.
        • Dardenne P.
        Capitalizing on fine milk composition for breeding and management of dairy cows.
        J. Dairy Sci. 2016; 99 (26778306): 4071-4079
        • Goddard M.
        • Hayes B.
        Mapping genes for complex traits in domestic animals and their use in breeding programmes.
        Nat. Rev. Genet. 2009; 10: 381-391
        • Hsu W.L.
        • Garrick D.J.
        • Fernando R.L.
        The accuracy and bias of single-step genomic prediction for populations under selection.
        G3 (Bethesda). 2017; 7 (28642364): 2685-2694
        • International Dairy Federation
        Dairy's global impact.
        https://fil-idf.org/dairys-global-impact
        Date: 2016
        Date accessed: November 10, 2021
        • Malacarne M.
        • Formaggioni P.
        • Franceschi P.
        • Summer A.
        Seasonal variations of milk quality in Parmigiano-Reggiano cheese manufacture on a period of 10 years.
        Sci. Tecn. Latt. Cas. 2004; 55: 63-77
        • Martin P.
        • Szymanowska M.
        • Zwierzchowski L.
        • Leroux C.
        The impact of genetic polymorphisms on the protein composition of ruminant milks.
        Reprod. Nutr. Dev. 2002; 42 (12537255): 433-459
        • Meyer K.
        WOMBAT - A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML).
        J. Zhejiang Univ. Sci. B. 2007; 8 (17973343): 815-821
        • Pretto D.
        • Lopez-Villalobos N.
        • Penasa M.
        • Cassandro M.
        Genetic response for milk production traits, somatic cell score, acidity and coagulation properties in Italian Holstein-Friesian population under current and alternative selection indices and breeding objectives.
        Livest. Sci. 2012; 150: 59-66
        • Sanchez M.P.
        • El Jabri M.
        • Minéry S.
        • Wolf V.
        • Beuvier E.
        • Laithier C.
        • Delacroix-Buchet A.
        • Brochard M.
        • Boichard D.
        Genetic parameters for cheese-making properties and milk composition predicted from mid-infrared spectra in a large dataset of Montbéliarde cows.
        J. Dairy Sci. 2018; 101 (30197141): 10048-10061
        • Sanchez M.P.
        • Fritz S.
        • Patry C.
        • Delacroix-Buchet A.
        • Boichard D.
        Frequencies of milk protein variants and haplotypes estimated from genotypes of more than 1 million bulls and cows of 12 French cattle breeds.
        J. Dairy Sci. 2020; 103 (32773310): 9124-9141
        • Sanchez M.P.
        • Govignon-Gion A.
        • Croiseau P.
        • Fritz S.
        • Hozé C.
        • Miranda G.
        • Martin P.
        • Barbat-Leterrier A.
        • Letaïef R.
        • Rocha D.
        • Brochard M.
        • Boussaha M.
        • Boichard D.
        Within-breed and multi-breed GWAS on imputed whole-genome sequence variants reveal candidate mutations affecting milk protein composition in dairy cattle.
        Genet. Sel. Evol. 2017; 49 (28923017): 68
        • Sanchez M.P.
        • Ramayo-Caldas Y.
        • Wolf V.
        • Laithier C.
        • El Jabri M.
        • Michenet A.
        • Boussaha M.
        • Taussat S.
        • Fritz S.
        • Delacroix-Buchet A.
        • Brochard M.
        • Boichard D.
        Sequence-based GWAS, network and pathway analyses reveal genes co-associated with milk cheese-making properties and milk composition in Montbéliarde cows.
        Genet. Sel. Evol. 2019; 51 (31262251): 34
        • Sanchez M.P.
        • Wolf V.
        • El Jabri M.
        • Beuvier E.
        • Rolet-Répécaud O.
        • Gaüzère Y.
        • Minéry S.
        • Brochard M.
        • Michenet A.
        • Taussat S.
        • Barbat-Leterrier A.
        • Delacroix-Buchet A.
        • Laithier C.
        • Fritz S.
        • Boichard D.
        Short communication: Confirmation of candidate causative variants on milk composition and cheesemaking properties in Montbéliarde cows.
        J. Dairy Sci. 2018; 101 (30219425): 10076-10081
        • Sargolzaei M.
        • Chesnais J.
        • Schenkel F.
        A new approach for efficient genotype imputation using information from relatives.
        BMC Genomics. 2014; 15 (24935670): 478
        • Soyeurt H.
        • Bruwier D.
        • Romnee J.
        • Gengler N.
        • Bertozzi C.
        • Veselko D.
        • Dardenne P.
        Potential estimation of major mineral contents in cow milk using mid-infrared spectrometry.
        J. Dairy Sci. 2009; 92 (19447976): 2444-2454
      3. Tribout, T., V. Ducrocq, and D. Boichard. 2020. HSSGBLUP: A single-step SNP BLUP genomic evaluation software adapted to large livestock populations. In ICQG, Brisbane, Australia.

        • VanRaden P.M.
        Invited review: Selection on net merit to improve lifetime profit.
        J. Dairy Sci. 2004; 87 (15377590): 3125-3131