Comparison of the genetic characteristics of directly measured and Fourier-transform mid-infrared-predicted bovine milk fatty acids and proteins

Fourier-transform mid-infrared (FT-MIR) spectroscopy is a high-throughput and inexpensive methodology used to evaluate concentrations of fat and protein in dairy cattle milk samples. The objective of this study was to compare the genetic characteristics of FT-MIR predicted fatty acids and individual milk proteins with those that had been measured directly using gas and liquid chromatography methods. The data used in this study was based on 2,005 milk samples collected from 706 Holstein-Friesian × Jersey animals that were managed in a seasonal, pasture-based dairy system, with milk samples collected across 2 consecutive seasons. Concentrations of fatty acids and protein fractions in milk samples were directly determined by gas chromatography and high-performance liquid chromatography, respectively. Models to predict each directly measured trait based on FT-MIR spectra were developed using partial least squares regression, with spectra from a random selection of half the cows used to train the models, and predictions for the remaining cows used as validation. Variance parameters for each trait and genetic correlations for each pair of measured/predicted traits were estimated from pedigree-based bivariate models using REML procedures. A genome-wide association study was undertaken using imputed whole-genome sequence, and quantitative trait loci (QTL) from directly measured traits were compared with QTL from the corresponding FT-MIR predicted traits. Cross-validation prediction accuracies based on partial least squares for individual and grouped fatty acids ranged from 0.18 to 0.65. Trait prediction accuracies in cross-validation for protein fractions were 0.53, 0.19, and 0.48 for α-casein, β-casein, and κ-casein, 0.31 for α-lactalbumin, 0.68 for β-lactoglobulin, and 0.36 for lactoferrin. Heritability estimates for directly measured traits ranged from 0.07 to 0.55 for fatty acids; and from 0.14 to 0.63 for individual milk proteins. For FT-MIR predicted traits, heritability estimates were mostly higher than for the corresponding measured traits, ranging from 0.14 to 0.46 for fatty acids, and from 0.30 to 0.70 for individual proteins. Genetic correlations between directly measured and FT-MIR predicted protein fractions were consistently above 0.75, with the exceptions of C18:0 and C18:3 cis -3, which had genetic correlations of 0.72 and 0.74, respectively. The GWAS identified trait QTL for fatty acids with likely candidates in the DGAT1 , CCDC57 , SCD , and GPAT4 genes. Notably, QTL for SCD were largely absent in the FT-MIR predicted traits, and QTL for GPAT4 were absent in directly measured traits. Similarly, for directly measured individual proteins, we identified QTL with likely candidates in the CSN1S1 , CSN3 , PAEP , and LTF genes, but the QTL for CSN3 and LTF were absent in the FT-MIR predicted traits. Our study indicates that genetic correlations between directly measured and FT-MIR predicted fatty acid and protein fractions are typically high, but that phenotypic variation in these traits may be underpinned by differing genetic architecture.


INTRODUCTION
Bovine milk is a rich source of dietary nutrients that are important to human health, including proteins, fats, carbohydrates, vitamins, and minerals. The concentrations of these components are determined by genetic factors such as breed and sire, as well as nongenetic factors related to the environment, stage of lactation, feed, and the nutritional status of the animal. Fats are important to human health due to the role they play in growth, development, hormone regulation, and inflammation management. In bovine milk, a typical fatty acid profile comprises about 70% saturated, 25% monounsaturated, and 5% polyunsaturated fatty acids.
Bovine milk is also a common source of protein, an important nutrient in the human diet because of the role it has in body maintenance, as well as the growth and repair of cells. However, the concentrations of casein and whey proteins in bovine milk differ to that of human milk, with bovine milk protein comprising approximately 80% casein and 20% whey proteins, whereas most of the protein in human milk represents whey proteins. These differences in protein composition are important because casein and whey proteins have different digestibilities and AA profiles. Moreover, the protein profiles have implications for cheese processing and the manufacture of casein supplements.
Fourier-transform mid-infrared (FT-MIR) spectroscopy is a method to determine the presence of specific chemical bonds in a composite substance such as milk, and is widely used in the dairy industry to characterize milk composition. The approach involves directing infrared light through a milk sample, leading to interactions between the infrared light and molecules in the milk that cause vibrations and rotational changes in molecular bonds, resulting in the differential absorption of the various infrared light wavelengths. From this process, a spectrum of absorbance values for light wavelengths across the mid-infrared range is generated, which can be used to predict a variety of traits. This is a high-throughput and inexpensive method for predicting milk composition from milk samples and is widely used to reliably quantify concentrations of fat and protein for dairy cattle. This methodology is also of interest for characterizing fat composition, casein, and whey proteins in milk because of the implications these milk components may have for human health and milk processability, and because the FT-MIR spectra are already available from routine milk testing.
Applications using FT-MIR spectral data to predict milk composition traits typically involve using a set of samples with directly measured trait values to develop a calibration equation based on the spectrum of absorbance values, using methods such as partial least squares (PLS) regression. The resulting calibration equation can then be applied to future samples to predict trait values as a linear combination of individual wavenumber absorbances from any milk sample with FT-MIR spectral data. The success of using FT-MIR data as a phenotyping tool relies on the strength of the phenotypic correlation between the directly measured trait and the FT-MIR predicted trait. However, the success of using an FT-MIR predicted trait in breeding programs is further dependent on the heritability of the predicted trait, and the genetic correlation between the directly measured and predicted trait.
Several GWAS have been conducted on fatty acids and protein fractions in bovine milk, across a range of genotype densities. This includes studies of directly measured fatty acids using 50k (Bouwman et al., 2011) or high-density (HD) genotypes (Buitenhuis et al., 2014;Palombo et al., 2018), and FT-MIR predicted fatty acids using 50k (Cruz et al., 2019;Iung et al., 2019;Freitas et al., 2020), HD (Olsen et al., 2017), or imputed whole-genome sequence  genotypes. Studies of directly measured protein fractions include those using 50k (Schopen et al., 2011;Pegolo et al., 2018) or HD (Buitenhuis et al., 2016;Zhou et al., 2019) genotypes, and studies of FT-MIR predicted protein fractions include those using imputed sequence genotypes (Sanchez et al., 2017b. Aside from differences in genotype density, the breed composition of animals in these studies also varies. In particular, studies of directly measured fatty acids include Dutch Holstein-Friesians (Bouwman et al., 2011), Danish Holsteins and Jerseys (Buitenhuis et al., 2014), and Italian Simmental and Holsteins (Palombo et al., 2018), whereas studies of FT-MIR predicted fatty acids include Holstein (Cruz et al., 2019;Iung et al., 2019;Freitas et al., 2020), Norwegian Red (Olsen et al., 2017), and Montbéliarde  cows. Studies of directly measured protein fractions in milk include Dutch Holstein-Friesians (Schopen et al., 2011), Italian Brown Swiss cows (Pegolo et al., 2018), and Danish Holsteins and Jerseys (Buitenhuis et al., 2016), whereas studies of FT-MIR predicted protein fractions include Montbéliarde, Normande, and Holstein cows (Sanchez et al., 2017b. Differences in genotype density and breed composition for GWAS conducted on directly measured and FT-MIR predicted fatty acid and protein traits make it difficult to compare QTL between studies. To date, as far as we are aware, there have been no GWAS that compare QTL for directly measured fatty acids and protein traits to QTL for the corresponding FT-MIR predicted traits within the same study population. The objective of this study was to compare the genetic characteristics of directly measured fatty acids and protein fractions to the same traits predicted from FT-MIR spectra. Calibration equations were developed using milk samples from New Zealand crossbred dairy cattle, and pedigree-based models were used to evaluate the (co)variance parameters of each directly measured trait and its corresponding FT-MIR predicted trait. To understand the underlying differences in the genetic architecture of directly measured and FT-MIR predicted traits, we conducted GWAS using imputed whole-genome sequence, and compared QTL from directly measured traits to QTL from the corresponding FT-MIR predicted traits. It was expected that the use of imputed whole-genome sequence genotypes from an F 2 study population would enhance our ability to identify trait QTL and candidate causative mutations, and that using the same data set to conduct GWAS across directly measured and FT-MIR predicted traits would be valuable for determining differences between QTL.

Ethics Statement
Animal ethics approval for the collection of data used in this study was granted by the Ruakura Animal Ethics Committee (Hamilton,New Zealand;approval numbers 4,232,4,621,and 10,174), according to the rules and guidelines outlined in the New Zealand Animal Welfare Act 1999.

Study Population/Animals and Milk Samples
Animals included in this study were from an F 2 design crossbreeding experiment with a half-sibling family structure, as previously described (Spelman et al., 2001;Berry et al., 2010). Briefly, 6 F 1 bulls were generated from reciprocal crosses of Holstein-Friesian and Jersey animals that were then mated to high genetic merit F 1 cows. This resulted in a herd of 850 F 2 female progeny, consisting of 2 cohorts produced over consecutive seasons, which were managed in a seasonal, pasture-based dairy system. Because of the phenotypic differences between milk composition for Friesian and Jersey animals, it was expected that the genetic variation exhibited in F 2 animals would typically be higher compared with what would be seen in a study of purebred animals, and that this could assist in the identification of trait QTL.
Measurements of FT-MIR spectra, and fatty acid and protein composition, were evaluated from second lactation milk samples collected at peak-, mid-and late-lactation in the 2003 to 2004 season for cohort 1, and the 2004 to 2005 season for cohort 2. Calving for each cohort took place over ~3 mo between July and October. Samples for each cohort representing peak milk were collected on a daily basis for these cows at 35 d postcalving, whereas mid-and late-lactation samples were collected at a fixed date across the herd within the season. A frequency distribution of the number of samples classified by DIM at the time of sampling has been provided in Appendix A1.
Concentrations of fatty acids were directly determined in milk fat samples by fatty acid methyl ester analysis using GC (MacGibbon and Reynolds, 2011), within 1 of up to 5 batches on a given sample collection day, and were expressed as grams per 100 g of total fat content. In this study, we report an analysis for 17 individual fatty acids and 6 fatty acid groups that were classified based on the degree of saturation and the length of the carbon chain, as follows: (1) SFA (no double bonds); (2) UFA (1 or more double bonds); (3) PUFA (2 or more double bonds); (4) short-chain fatty acids (SCFA; 4, 6, or 8 carbons); (5) mediumchain fatty acids (MCFA; 10, 12, or 14 carbons); and (6) long-chain fatty acids (LCFA; 18 carbons). Milk proteins were determined using HPLC, as described by Palmano and Elgar (2002), and were analyzed within 1 of up to 6 batches on a given sample collection day, and were expressed as grams per liter of total milk volume. Traits were assessed for deviation from normality by visual inspection of normal quantile plots and by evaluating asymmetry according to skewness. With the exception of lactoferrin, all directly measured traits were approximately normally distributed with absolute skewness values less than 1. For lactoferrin, log, square-, and cube-root transformations were applied to determine which transformation minimized skewness. A cube-root transformation was the most effective of those investigated for minimizing skewness and was applied to lactoferrin trait values for all downstream analyses. Frequency distributions of untransformed lactoferrin concentrations and lactoferrin concentrations after applying a cube-root transformation are provided in Appendix A2. Outliers for each fatty acid and protein trait were identified and removed if the trait value was more than 3 standard deviations from the mean for the corresponding season and stage of lactation (peak, mid, late). After removal of outliers, each trait was adjusted to remove batch effects, where batch effects were evaluated from a random effects model with batch nested within season and stage of lactation, using Nelder-Mead optimization as implemented in the lme4 package in R (Bates et al., 2015).
The same milk samples assessed for fatty acid and protein composition were also analyzed on a Foss MilkoScan FT6000 (Foss) instrument, to generate spectral records consisting of 1,060 wavenumbers across the range from 925.66 to 5,010.15 cm −1 . Spectral data from regions associated with low signal-to-noise ratios and poor sample measurement repeatability due to the water content in milk were excluded, according to the definitions by Tiplady et al. (2019). Specifically, the excluded low signal-to-noise regions were 649 to 970 cm −1 , 1,608 to 1,682 cm −1 , and ≥3,021 cm −1 . This resulted in 542 wavenumbers for use in the development of prediction equations. Outliers in the spectral data were identified using the methodology described in Tiplady et al. (2019). Briefly, the squared Mahalanobis distance between each spectral record and the average spectra were evaluated using the 542 wavenumbers identified as being outside low signal-to-noise regions. The distributions of Mahalanobis distance values for each season were compared and found to be similar, indicating that although the spectra were collected in 2 different seasons, the effect of instrument drift across time was likely to be small. Based on the lowest average information criterion, a logistic distribution with location and scale parameters of 541.7 and 27.3, respectively, had the best fit to the overall Mahalanobis distance values, and based on a P-value of 0.001, 18 outliers were identified and removed. In total, after outlier removal, we had 2,005 samples for 706 animals with FT-MIR spectra and either a fatty acid or protein composition result. Traits varied in the final number of records available for analysis, ranging from 1,686 to 1,977 records, and representing from 699 to 704 animals. The overall mean fat and protein concentrations as predicted from the Foss instrument calibration equation were 5.40 (SD = 0.70) and 3.98 (SD = 0.36), respectively.

Development and Validation of Calibration Equations
Phenotypic calibration equations for each fatty acid and protein fraction were evaluated within a cross-validation framework, whereby records for a random selection of half the animals were assigned to a training dataset, and the remaining records were assigned to a validation dataset. This ensured that validation was cow-independent in that none of the records for animals included in the training dataset were included in the validation dataset. Partial least squares models for each trait were developed using 542 spectral wavenumbers with the caret package in R (Kuhn et al., 2022), based on training data with 10 repeats of 10-fold cross-validation. In addition to the untreated spectra, several mathematical treatments of spectra were assessed using the mdatools package in R (Kucheryavskiy, 2020), in-cluding standard normal variate (SNV) transformation, multiplicative scatter correction, and first-order Savitzky-Golay derivative (Savitzky and Golay, 1964) treatments. First-derivative treatments were applied to untreated spectra and spectra after SNV or multiplicative scatter correction treatments using a range of window sizes, with up to 1 and 10 points either side. For each trait, the performance of the PLS model was assessed according to the coefficient of determination between actual and predicted phenotypic trait values in the validation dataset R cv 2 ( ) , and the relative prediction error (RPE) between actual and predicted trait values in the validation dataset (RPE cv ), as described by Lopez-Villalobos et al. (2014).

Genetic Parameters of Traits
Genetic (co)variances of each directly measured trait and its corresponding FT-MIR predicted trait were estimated using a pairwise bivariate repeated measures animal model in ASReml-R (Butler et al., 2009) based on a pedigree comprising 5,943 animals. The model was defined as follows: where y 1 is a vector of the directly measured fatty acid or protein fraction, y 2 is a vector of the corresponding FT-MIR predicted trait; X 1 , Z 1 , W 1 , X 2 , Z 2 , and W 2 are design matrices for the fixed, additive genetic and permanent environment effects, respectively, for y 1 and y 2 ; b 1 and b 2 are vectors of the fixed effect of DIM (represented as 35-day windows from the start of lactation) within season (2003,2004) for the directly measured and the FT-MIR predicted trait, respectively; u 1 and u 2 are vectors of random additive genetic effects for each trait; p 1 and p 2 are vectors of permanent environment effects for each trait; and e 1 and e 2 are vectors of residuals. The following (co)variance structure for each directly measured (y 1 ) and FT-MIR predicted (y 2 ) trait pair is assumed: of order corresponding to the length of the vector p, I e is an identity matrix of order corresponding to the length of the vector e, ⊗ is the Kronecker product. Additionally, G, C, and R are genetic, permanent environment and residual (co)variance matrices, respectively, and are defined as follows: residual covariances, respectively, where a and b ranged from 0.1 to 0.9 in increments of 0.1. Among models that converged for each pair of traits, genetic parameter estimates were highly consistent. For traits that had different solutions from different models, the model that minimized the squared sum of the difference between single-and multi-trait model heritability estimates was selected.

Genotypes and Imputation
Of the 706 animals with phenotypic data, 685 were genotyped on Illumina BovineHD (HD; n = 12; ~777k SNP) or Illumina BovineSNP50k (50k; n = 685; ~53k SNP) panels, or were genotyped on both. The resultant genotypes were imputed to sequence density as part of a wider set of 153,357 animals, as described previously (Jivanji et al., 2019;Tiplady et al., 2021). Briefly, the imputation process consisted of stepwise imputation of animals to whole-genome sequence genotypes via references of GeneSeek Genomic Profiler, 50k, and HD genotypes. The whole-genome sequence reference consisted of 565 animals, comprised of 138 Holstein-Friesians, 99 Jerseys, 316 Holstein-Friesian × Jersey crossbreeds, and 12 from other breeds or crosses. Notably, the 6 F 1 sires included in our study were included in this wholegenome sequence reference and were sequenced with a target of 60× read-depth coverage. Phasing was undertaken using Beagle 4.0 (Browning and Browning, 2007), based on genotype probabilities, and variants were filtered to remove those where the allelic R 2 for missing genotypes was less than 0.95. Only variants located on Bos taurus autosomes were considered, resulting in a sequence reference comprising 19,659,361 segregating variants spanning all 29 autosomes. Imputation was carried out using Beagle 4.0 (Browning and Browning, 2007), ignoring pedigree information, and SNP with allelic R 2 < 0.7 were removed after each imputation step. The overall median imputation allelic R 2 for the wider set of 153,357 animals was 0.986, but was 0.992 for the 685 genotyped animals included in this study.

Genome-Wide Association Studies
Before conducting GWAS, adjusted fatty acid and protein phenotypes were generated for directly measured and FT-MIR predicted traits. The generation of the adjusted phenotypes was based on 1 or more samples measured on the same cow, which were fitted to a univariate pedigree-based repeated measures model in ASReml-R (Butler et al., 2009), as follows: where y is a vector of the measured or predicted trait, X, Z, and W are design matrices for the fixed, additive genetic, and permanent environment effects; b is the fixed effect of DIM (represented as 35-d windows from the start of lactation) within season (2003,2004) for the trait; u is a vector of random additive genetic effects with u ∼ N 0 A the length of the vector e, σ u 2 is the additive genetic variance, σ p 2 is the permanent environment variance, and σ e 2 is the residual variance. Adjusted phenotypes used in the GWAS were the average of y over all observations for a cow minus the relevant fixed effects.
For each directly measured fatty acid or protein trait and its corresponding FT-MIR predicted trait, a GWAS was conducted using Bolt-LMM software (Loh et al., 2015). Before conducting GWAS, a minor allele frequency threshold of 1% based on allele frequencies in the 685-animal study population was applied, resulting in 14,990,779 imputed sequence variants included in each GWAS. To assess the additive effect of each SNP, mixed model association statistics were evaluated under an infinitesimal model. To account for population structure, a genomic relationship matrix based on a subset of 42,374 SNPs was simultaneously fitted. That subset of SNP was derived by applying a minor allele frequency threshold of 1% to the 50k SNP-chip imputation reference (previously described). A leaveone-segment-out approach was used to avoid proximal contamination in the GWAS, whereby a 5-Mbp region flanking the sequence variant of interest was excluded from the set of SNPs used to estimate the genomic relationship matrix.
An adjusted Bonferroni threshold was adopted to determine variants with significant associations for each trait. Because a Bonferroni correction threshold based on all 14,990,779 variants is highly conservative, a modified threshold was evaluated based on the effective number of independent variants, as proposed by Duggal et al. (2008) and implemented in other studies (Zhu et al., 2017;. The effective number of independent variants were identified using a sliding window approach in Plink software (Purcell et al., 2007), with an R 2 threshold of 0.9, a window size of 100 kb and a step size of 5 variants. These criteria resulted in a set of 2,303,435 variants and enabled the calculation of an adjusted Bonferroni threshold which considered all tests across 2,303,435 variants as independent. Based on α = 0.05, this resulted in a nominal P-value of 4.3e-09 and a corresponding Bonferroni threshold of −log 10 (4.3e-09) = 8.36. Whole-genome sequence resolution genotypes within a 1Mbp window were annotated using SnpEff (version 4.3t;build 11-24-2017;Cingolani et al., 2012) and the Ensembl UMD3.1.86 gene annotations to assess the candidacy of QTL identified from the GWAS for each trait. We used a linkage disequilibrium (LD)-based approach to prioritize variants, similar to that described by Lopdell et al. (2017) because the association rankings of candidate variants are expected to be affected by phenotyping, genotyping, and imputation errors. Specifically, we identified QTL regions where the most highly associated variant was in high LD (R 2 > 0.7) with either a splice region variant, or a moderate or high impact coding variant, according to SnpEff classification.

Trait Prediction Models
Cross-validation prediction model accuracies R cv 2 ( ) were assessed for untreated spectra, as well as for spectra treated using SNV transformation, multiplicative scatter correction, or first-derivative treatments (Appendix Table A1). Window sizes of 15 data points (7 points either side) had consistently higher R cv 2 values, compared to other window sizes, so only these have been presented. Applying treatments to spectral data resulted in marginally higher R cv 2 values on average, compared to not treating spectra, and treating spectra with a SNV and first-derivative transformation prior to fitting PLS models resulted in the highest average R cv 2 value and was thus used in all further analysis. Descriptive statistics of fatty acid and protein traits, and goodness of fit measures of PLS calibration models (applied to SNV + first-derivative transformed spectra) for training and validation datasets are presented in Table  1.
For individual fatty acids, coefficient of determination values for the validation dataset R cv 2 ( ) were generally higher for short-chain fatty acids (C4 to C8), ranging from 0.54 to 0.62, compared with medium-chain fatty acids (C10 to C14), which ranged from 0.30 to 0.63, and long-chain fatty acids (C16 to C18), which ranged from 0.18 to 0.57. Concentrations of individual saturated fatty acids were typically higher and had higher average R cv 2 values, compared with individual unsaturated fatty acids. For grouped fatty acids, R cv 2 values were higher for UFA and SFA groups, compared to PUFA; additionally, for fatty acids grouped by carbon chain length, the highest R cv 2 value of 0.65 was observed for SCFA. It is notable that although we found an overall trend of higher R cv 2 values coinciding with lower RPE cv values, there were exceptions to this. For example, among individual fatty acids, C16:1 had a particularly low R cv 2 of 0.18, but an RPE cv of 0.13, which was comparable to other traits such as C10:0 and C12:0, which had R cv 2 values of ~0.60. This highlights the difference between R cv 2 and RPE cv as accuracy metrics, the former indicating how well the prediction model explains the variation in the directly measured trait, whereas the latter provides a comparison of how similar the predicted values are to the directly mea-sured trait values. In the present study, most comparisons of accuracy with other studies will be based on R cv 2 values because that is the accuracy metric that is most commonly reported; however, the example above shows that other metrics can be valuable for assessing FT-MIR prediction model accuracy.
The R cv 2 values we report are consistent with those from previous studies where fatty acids were expressed as a proportion of total fat content, with our values being similar to those reported by Soyeurt et al. (2006), but lower than those reported in other studies (Rutten et al., 2009;Lopez-Villalobos et al., 2014;Bonfatti et al., 2016). In the present study, for grouped SCFA, MCFA, and LCFA, R cv 2 values were lower than in other studies (Rutten et al., 2009;Lopez-Villalobos et al., 2014;Bonfatti et al., 2016). Accuracies for fatty acids predicted using FT-MIR spectra were variable in previous studies and were affected by factors such as the production system and the breed composition diversity present in calibration samples, the number of samples used to develop calibration equations, and the variability of fatty acid composition present in the calibration samples. Rutten et al. (2009)   R t 2 = coefficient of determination between actual and predicted trait values in the training dataset; RPE t = relative prediction error between actual and predicted trait values in the training dataset. 3 R cv 2 = coefficient of determination between actual and predicted trait values in the validation dataset; RPE cv = relative prediction error between actual and predicted trait values in the validation dataset. 4 SCFA = short-chain fatty acids, sum of C4:0, C6:0, and C8:0; MCFA = medium-chain fatty acids, sum of 10:0, 10:1, 12:0, 12:1, 14:0, and 14:1; LCFA = long-chain fatty acids, sum of C18 fatty acids. that prediction accuracy could be improved by increasing the sample size of their study, and by increasing the range of variation present in the fatty acids. Importantly, studies with the highest accuracies were those where the range of fatty acid values present in the validation samples were encompassed within the range of fatty acid values represented in calibration samples.
For individual milk proteins, R cv 2 values were generally lower than for fatty acids, ranging from 0.19 for β-CN to 0.68 for β-LG. Notably, although the R cv 2 values for β-CN and β-LG were very different, the RPE cv values for these 2 traits were similar (0.11 and 0.10, respectively). The R cv 2 values we report for individual milk proteins were typically higher than those reported in previous studies of individual milk proteins, with the exceptions of β-CN and lactoferrin (Lf), which were consistently lower here than in other studies (De Marchi et al., 2009;Lopez-Villalobos et al., 2009;Rutten et al., 2011;Soyeurt et al., 2012;Bonfatti et al., 2016;McDermott et al., 2016). Fuentes-Pila et al. (1996) suggest that a RPE of lower than 0.1 is an indicator of satisfactory prediction, a RPE between 0.1 to 0.2 is an indicator of relatively good or acceptable predictions, and a RPE greater than 0.2 is an indicator of unsatisfactory prediction. Based on these criteria, 21 of 23 individual and grouped fatty acids and all 6 protein fractions had good or satisfactory predictions in the validation datasets. Although the guidelines proposed by Fuentes-Pila et al. (1996) are useful as an indicator of prediction acceptability, they are potentially less meaningful when we are considering the value of incorporating FT-MIR predicted traits into animal breeding programs. This is because FT-MIR predictions can provide indicator traits across large numbers of animals at little or no cost, whereas it may be infeasible to directly measure these traits across even a small number of animals. Moreover, when we are considering the potential for incorporating an FT-MIR predicted trait into a breeding program, we are not only interested in the phenotypic correlation between the directly measured and FT-MIR predicted trait, but also the heritability of the FT-MIR predicted trait, and the genetic correlation between the directly measured and FT-MIR predicted trait.

Genetic Parameters of Directly Measured and FT-MIR Predicted Traits
Estimates of variance components for directly measured and FT-MIR predicted fatty acid and protein traits are shown in Table 2 and Appendix Table A2. Heritability estimates (h 2 ) for the majority of traits were moderate to high, with 17 of the directly measured traits and 20 of the FT-MIR predicted traits having an estimated heritability greater than 0.3. Because this is an F 2 study, genetic variances will include a segregation variance component that would typically inflate these values compared to what would be seen in a study of purebred animals. In general, lower heritability and repeatability estimates were observed for directly measured traits, compared to FT-MIR predicted traits. This was driven by higher total variation σ T 2 ( ) in the directly measured traits, coupled with a lower magnitude increase in the additive genetic variance component σ u 2 ( ) , compared to the FT-MIR predicted traits.
Despite this, the genetic correlations between measured and predicted traits remained high and were mostly greater than 0.75.
Fatty Acid Traits. In fatty acid traits, the lowest heritability estimates were observed for C18:0 and LCFA, with heritability estimates of 0.07 for the directly measured traits, and heritability estimates of 0.14 and 0.15 in the FT-MIR predicted traits, respectively. Although heritability estimates were typically higher in the FT-MIR predicted traits, there were exceptions to this. In particular, C14:1 had an estimated heritability for the measured trait that was substantially higher than that of the FT-MIR predicted trait (0.55 vs. 0.26). Genetic correlations between directly measured and FT-MIR predicted traits (r a ) were greater than 0.85 for 18 of 23 individual and grouped fatty acids, and for 11 of these traits, the genetic correlation was greater than 0.95. The lowest genetic correlations were observed for C18:0 (r a = 0.72) and C18:3 cis-3 (r a = 0.74). In general, we found a consistent trend for individual and grouped fatty acids, where lower genetic correlations coincided with low R cv 2 values.
Although several studies have reported genetic parameter estimates for directly measured or FT-MIR predicted fatty acid traits, or both, these studies vary in the specific individual fatty acids (if any) presented, and whether or not they present parameter estimates for grouped fatty acids. Many studies report genetic parameter estimates for FT-MIR predicted traits only (Soyeurt et al., 2007b;Lopez-Villalobos et al., 2014;Narayana et al., 2017;Fleming et al., 2018), with only 2 studies reporting genetic parameters for both directly measured and FT-MIR predicted traits (Rutten et al., 2010;Bonfatti et al., 2017b). These latter 2 studies also report genetic correlations between directly measured and FT-MIR predicted fatty acids, with Bonfatti et al. (2017b) presenting these for individual and grouped fatty acids, whereas Rutten et al. (2010) presented these for individual fatty acids only.
The heritability, repeatability, and genetic correlation estimates we report in the present study were broadly Tiplady et al.: BOVINE MILK FATTY ACIDS, PROTEINS: GENETIC CHARACTERISTICS consistent with those from previous studies. For directly measured fatty acids, the heritability estimates we report were typically higher than those reported by Bonfatti et al. (2017b), but lower than those reported by Rutten et al. (2010). For FT-MIR predicted fatty acids, the heritability and repeatability estimates we report for individual and grouped fatty acids were similar to those presented by Lopez-Villalobos et al. (2014), but lower than those presented by Narayana et al. (2017) and higher than those presented in other studies (Soyeurt et al., 2007b;Bonfatti et al., 2017b). Compared with other studies that report genetic correlations between directly measured and FT-MIR predicted fatty acids, the genetic correlations we report were similar, with standard errors of a similar magnitude (Rutten et al., 2010;Bonfatti et al., 2017b). The moderate to high heritability estimates we report, alongside high genetic correlations between directly measured and FT-MIR predicted fatty acid traits, indicate that there is genetic variation in the FT-MIR predicted traits that could potentially be exploited in animal breeding programs, and, in most cases, that selection for an FT-MIR predicted fatty acid trait would be expected to provide genetic gain in the actual fatty acid trait of interest.
Individual Milk Protein Traits. Heritability estimates were moderate to high for nearly all directly measured and FT-MIR predicted individual milk proteins ( Table 2). The exception to this was for directly measured β-CN, which had a heritability of 0.14. The highest heritability estimates were for β-LG, with h 2 = 0.63 and h 2 = 0.70 for directly measured and FT-MIR predicted β-LG, respectively. In general, heritability estimates for measured and FT-MIR predicted proteins were similar. An exception to this was β-CN, which had  heritability estimates for the directly measured and FT-MIR predicted trait of 0.14 and 0.38, respectively. Another exception was Lf, which had an estimated heritability for the measured trait that was substantially higher than that of the FT-MIR predicted trait (0.59 vs. 0.30). With the exceptions of α-LA and Lf, genetic correlations between directly measured and FT-MIR predicted individual milk proteins were greater than 0.8. In general, as we observed for fatty acid traits, we found a trend of low R cv 2 values, coinciding with low genetic correlations between directly measured and FT-MIR predicted traits.
There are few studies that report genetic parameters for directly measured or FT-MIR predicted milk proteins, or both, but those studies vary in the breed composition of the cows. Specifically, study populations include Dutch Holstein-Friesians (Schopen et al., 2009), Danish Holsteins and Jerseys (Buitenhuis et al., 2016), Italian Simmentals (Bonfatti et al., 2017b), or French Montbéliarde, Normande, and Holstein cows (Sanchez et al., 2017a). Studies also vary in that some report on individual proteins as a proportion of total protein or whey protein (Schopen et al., 2009;Buitenhuis et al., 2016), whereas other studies report on individual proteins as a proportion of total protein or as a proportion of total milk volume (Bonfatti et al., 2017b;Sanchez et al., 2017a). The heritability estimates we report for directly measured α-CN, β-CN, and κ-CN were lower than those previously reported by Bonfatti et al. (2017b), but the heritability estimates we report for directly measured α-LA and β-LG were substantially higher. In contrast, for FT-MIR predicted α-CN, β-CN, and κ-CN, the heritability estimates we report were consistently higher than those reported by Bonfatti et al. (2017b), but were similar to those reported by Sanchez et al. (2017a).
The only study to report genetic correlations between directly measured and FT-MIR predicted milk proteins that we are aware of is that of Bonfatti et al. (2017b). The genetic correlations that we report were higher than in that study. Specifically, for the protein fractions we studied, genetic correlations ranged from 0.76 for α-LA to 0.995 for β-LG, whereas in Bonfatti et al. (2017b), genetic correlations for these traits ranged from 0.24 for α-LA to 0.48 for β-LG. Interestingly, Bonfatti et al. (2017b) reported moderate heritability estimates for directly measured milk proteins (0.12 to 0.59), but much lower heritability estimates for FT-MIR predicted milk proteins (0.07 to 0.21). In contrast, the heritability estimates we observed for directly measured proteins (0.14 to 0.63) were similar to (and often lower than) the heritability estimates we observed for FT-MIR predicted proteins (0.30 to 0.70). These differ-ences in heritability were likely due to factors related to differences in the breed composition and population structure between the 2 studies (i.e., Italian Simmental cows from herds enrolled in the Italian national milk recording program versus Holstein-Friesian Jersey F 2 cows from a single research herd).
Moderate to high heritability estimates and high genetic correlations between directly measured and FT-MIR predicted milk proteins in our study indicate that indirect selection on FT-MIR predicted milk proteins could be used in animal breeding programs to achieve desired changes to milk protein composition. Moreover, high genetic correlations from pedigree-based models imply that directly measured and FT-MIR predicted traits may have a similar underlying genetic architecture and that genes contributing to the traits are likely to be co-inherited (Lynch and Walsh, 1998). To assess this directly, we conducted GWAS on directly measured traits and their corresponding FT-MIR predictions, and compared the QTL for each trait.

Sequence-Based Genome-Wide Association Analyses
Previously, there have been several GWAS that used a range of genotype densities for fatty acids in bovine milk samples determined by GC (Bouwman et al., 2011;Buitenhuis et al., 2014;Palombo et al., 2018) or fatty acids predicted using FT-MIR spectra (Olsen et al., 2017;Cruz et al., 2019;Iung et al., 2019;Sanchez et al., 2019;Freitas et al., 2020). Similarly, multiple GWAS have been conducted on protein fractions in milk samples determined by HPLC (Schopen et al., 2011;Buitenhuis et al., 2016;Pegolo et al., 2018;Zhou et al., 2019) or FT-MIR predicted protein fractions (Sanchez et al., 2017b. Each of those studies was conducted using either the directly measured trait (GC-based for fatty acids; HPLC-based for protein fractions) or the FT-MIR predicted trait, though none of these presented comparisons between the GWAS for directly measured and FT-MIR predicted traits. In the present study, we have sought to make these comparisons using imputed whole-genome sequence genotypes from an F 2 study population to enhance our ability to identify trait QTL and candidate causative mutations.
For each of 17 individual fatty acids, 6 grouped fatty acids, and 6 protein traits, GWAS were conducted using 14,990,779 imputed sequence variants. These analyses resulted in the identification of 40,946 variants with significant effects for directly measured traits, and 18,843 variants with significant association effects for the FT-MIR predicted traits. We found more than twice as many variants with significant effects for directly measured traits, compared with FT-MIR predicted traits, which was largely due to 20,949 variants with significant effects on BTA26 for directly measured traits compared with only 110 variants with significant effects on BTA26 for FT-MIR predicted traits. It was also notable that we detected 3,579 variants with significant effects on BTA22 for directly measured Lf, but no variants with significant effects on BTA22 for FT-MIR predicted traits. Manhattan plots showing the strength of association signals are presented in Figures 1-4 for individual fatty acids, Figure 5 for grouped fatty acids, and Figure 6 for individual protein traits. To assess the candidacy of QTL, relevant protein coding variants that were in high LD (R 2 > 0.7) with the most highly associated variant from each peak were identified. The most highly associated variant from each trait QTL and any relevant protein coding variants are shown in Table 3 for directly measured fatty acid and protein traits, and Table 4 for FT-MIR predicted fatty acid and protein traits. Effect sizes and minor allele frequency details for relevant variants and effects are provided in Appendix Table A3 for fatty acids and Appendix Table  A4 for protein traits.
Short-Chain Fatty Acids. Prominent peaks were observed on BTA17 for the short-chain fatty acids, C4:0, and C6:0 (Tables 3 and 4; Figure 1). For directly measured and FT-MIR predicted C4:0, these peaks were underpinned by the same QTL at chromosome (Chr) 17:53.03 Mbp (rs461037541). A peak of similar magnitude was also observed for FT-MIR predicted C6:0 at a nearby locus (rs207997694), with a less significant peak for directly measured C6:0 at that same locus. Other significant effects were also observed at this locus for directly measured SCFA (P-value = 1.2e-14) and FT-MIR predicted SCFA (P-value = 7.1e-22). The 2 implicated loci for the peaks on BTA17 were situated between the AACS and BRI3BP genes, and visual examination revealed several significant variants across both genes. The AACS gene codes for the enzyme acetoacetyl-CoA synthetase, which forms an important metabolic link between the ketone body acetoacetate on one hand, and the tricarboxylic acid cycle and fat synthesis on the other (Bergman, 1971). As this gene is highly expressed in both adipose and mammary tissue (NCBI Bioprojects PRJEB4337 and PRJEB2445), AACS makes a good candidate for the causal gene underlying fatty acid QTL in this region. Knutsen et al. (2018) also reported an effect for C4:0 fatty acids in this region and suggested that the QTL may be the result of a regulatory effect.
Medium-Chain Fatty Acids. Significant effects were observed on BTA11, BTA19, and BTA26 for medium-chain fatty acids (Tables 3 and 4; Figure 2). The peak on BTA11 was underpinned by a Chr11: 103 .30 locus (rs41255687) and was observed for FT-MIR predicted C12:1, but was absent for directly measured C12:1. This locus was in high LD (R 2 > 0.98) with 2 missense mutations in the PAEP gene, which encodes the major whey protein, β-LG. One of the missense mutations reported (rs109625649; V134A) is a variant that distinguishes the 'A' and 'B' haplotypes of β-LG (Caroli et al., 2009), where the 'A' haplotype is known to be associated with higher levels of β-LG expression. The PAEP locus has also been linked to FT-MIR wavenumbers characterized by carboxylic C=O bond stretching . This type of bond is found in both fats and proteins, strongly suggesting that the peak observed for the FT-MIR predicted phenotype is a false positive due to contamination of the signal by varying levels of β-LG expression.
Several QTL were identified for directly measured and FT-MIR predicted medium-chain fatty acids (C10:0, C12:0, C14:0) on BTA19 that were in high LD (R 2 > 0.97), with a missense mutation (rs41921160) in the CCDC57 gene (Tables 3 and 4; Figure 2). Significant effects were also observed in this region for FT-MIR predicted C8:0 (P-value = 8.9e-10; Figure 1), and directly measured (P-value = 1.4e-13) and FT-MIR predicted MCFA (P-value = 9.2e-13; Figure 5). The CCDC57 encodes a coiled-coil domain-containing protein that is expressed in the bovine mammary gland (Medrano et al., 2010). Previous studies have implicated the same or a nearby locus to the one reported here as having a significant association for fatty acids (Bouwman et al., 2014;Knutsen et al., 2018;Palombo et al., 2018) and fat composition (Tribout et al., 2020) in bovine milk. Significant effects have also been reported at a nearby locus for several FT-MIR wavenumbers, characterized by carboxylic C=O bond stretching . Bouwman et al. (2014) examined this region in depth using HD genotypes and identified 2 possible genes underlying an effect for C14:0, CCDC57 and FASN. The missense mutation we have highlighted (rs41921160) is located within the same region as the most highly associated variants in the study by Bouwman et al. (2014), and was in perfect LD with the set of 8 intronic HD variants with the most highly associated effects. On closer examination of the association effects for C10:0, C12:0, and C14:0 in our study, we determined that alongside the most highly associated variants in the QTL peaks, there were 47 other imputed whole-genome sequence variants between 51,306,219 and 51,330,072 bp that were in perfect LD with one another (including the missense variant rs41921160), with only marginally less significant P-values. A small cluster of association effects near to or in the FASN gene were also observed, with the most significant of these for directly measured C14:0 being at 51,380,689 bp, but the P-value for that predicted individual short-chain fatty acid traits. Dark and light blue data points represent association signals for GC-based traits and red data points represent association signals for FT-MIR predicted traits. Chromosomes and genomic position based on the UMD3.1 Bos taurus reference genome are represented on the x-axis. The strength of association signals is represented as the −log 10 (P-value) on the y-axis. The horizontal red line shows the Bonferroni significance threshold of −log 10 (4.3e-09). effect was not significant (P-value = 2.4e-07). To assess whether multiple QTL were present in this region, we repeated the GWAS, correcting for the rs136424304 locus by including it as a covariate in the Bolt-LMM model. This resulted in no significant effects remaining in a 1 Mbp region around the Chr19: 51 .32 locus, and the association effect near the FASN gene at 51,380,689 bp, dropping in significance to a P-value of 3.9e-02. Although our analysis provides evidence that the effect in this region may be underpinned by a missense mutation in the CCDC57 gene, the functional candidacy of FASN remains and such an effect would need to be confirmed by functional analysis.
Multiple QTL were identified for directly measured medium-chain fatty acids on BTA26 (Table 3; Figure  2). The most significant effects were observed at Chr26: 21 .15 Mbp for directly measured C10:1 (rs41255688; P-value=1.8e-48) and C14:1 (rs385285356; P-value = 6.1e-61). These loci were in high LD (R 2 = 0.92) with a splice region variant (rs41255693) in the SCD gene. The SCD gene was also identified as encompassing other effects with less significant P-values for directly measured C10:0, C14:0, SFA, and UFA (Table 3), and FT-MIR predicted UFA (Table 4). Stearoyl-CoA desaturase is an enzyme that plays an important role in biosynthesis of monounsaturated fatty acids (Bernard et al., 2006;Paton and Ntambi, 2009), and has previously been reported in other studies of fatty acids in bovine milk (Mele et al., 2007;Moioli et al., 2007;Schennink et al., 2008;Kgwatalala et al., 2009;Conte et al., 2010;Bouwman et al., 2011). The strong effect we see for directly measured C14:1 in the SCD gene is unsurprising, given that C14:0 in milk fat is predominantly derived from de novo synthesis in the mammary gland, meaning that almost all C14:1 cis-9 is likely to have been synthesized by stearoyl-CoA desaturase (Bernard et al., 2006). Interestingly, although we found a significant effect for FT-MIR predicted UFA at a nearby locus that was also in high LD with the rs41255693 splice region variant (R 2 = 0.91), no other effects were identified within the SCD gene for individual FT-MIR predicted fatty acids. A peak for FT-MIR predicted C14:1 was tagged by a nearby Chr26: 21 .17 Mbp locus (rs209445650 ; Table 4). However, the LD between the rs209445650 locus and the splice region variant identified for directly measured fatty acids (rs41255693) was moderately low (R 2 = 0 .32). Moreover, in a recent GWAS of individual FT-MIR wavenumbers, there was no evidence of an association effect linked to the SCD gene , indicating that changes in milk composition due to this gene may be difficult to detect using FT-MIR spectral data. However, we may also view the absence of FT-MIR predicted trait QTL in the SCD gene within the context of trait prediction accuracy. The largest QTL underpinned by SCD in directly measured fatty acids were for C10:1 (P-value = 1.8e-48) and C14:1 (P-value = 6.1e-61). The prediction accuracies for these traits were relatively poor: C10:1 R cv 2 ( = 0.30; RPE cv = 0.16) and C14:1 R cv 2 ( = 0.41; RPE cv = 0.23; Table 1). Also, it is notable that for C10:1 and C14:1, the heritability estimates of the FT-MIR predictions were lower than those for direct measurements of these traits. This contrasts with the typical pattern for nearly all other fatty acids where the heritability for the FT-MIR prediction was greater than the heritability for the directly measured trait. In particular, the heritability estimate for directly measured C14:1 was 0.55, whereas the heritability estimate for FT-MIR predicted C14:1 was 0.26 (Table 2). Low prediction accuracy and a substantially lower heritability estimate for FT-MIR predicted C14:1 may in part be explained by C14:1 being at relatively low concentrations in milk samples, particularly compared with saturated fatty acids. Specifically, C14:1 had a mean concentration of 0.75 g/100 g of total fat, compared to mean concentrations of 1.54 to 27.64 g/100 g of total fat for the individual saturated fatty acids included in this study (Table 1). Potentially, it may be possible to improve trait prediction accuracies, herita- bility estimates, and QTL identification for C14:1 by basing FT-MIR predictions on the ratio of C14:1 to C14:0, as in the study by Arnould et al. (2009a). In that study, they highlight that genetic variation and heritability estimates change throughout lactation for the ratio of C14:1 to C14:0, so it may also be valuable to examine other methods of accounting for stage of lactation such as using Legendre polynomials within random regression models.
One further QTL was observed for directly measured C10:1 and C12:1 on BTA26 at a Chr26: 26 .46 Mbp locus (rs445758306; Table 3; Figure 2). This locus was in high LD (R 2 = 0.76) with a missense mutation in the ITP-RIP gene (rs379463458). The ITPRIP gene modulates intracellular messaging by binding the inositol 1,4,5-triphosphate receptor ITPR. This gene has not previously been reported in GWAS of bovine milk composition, and the potential role it may play in the regulation of bovine milk fatty acids is unclear. An alternative potential candidate gene that the Chr26: 26 .46 Mbp locus maps close to is SORCS3, which encodes the sortilinrelated receptor SorCS3. Sortilins are involved in regulating glucose transport into cells in response to insulin (Huang et al., 2013). A potential mechanism by which this gene could influence milk fatty acid concentrations is via changing the supply of glucose available for the pentose phosphate pathway, which in turn provides the nicotinamide adenine dinucleotide phosphate necessary for fatty acid synthesis.
Long-Chain Fatty Acids. Two QTL were identified on BTA14 for directly measured individual long-chain fatty acids (Table 3; Figures 3 and 4). One of these was at a Chr14: 1 .80 Mbp (rs385135066) locus that had a significant effect for directly measured C16:0 (P-value = 1.2e-12). This locus was in high LD (R 2 = 0.74) with missense mutations in the SLC52A2 and DGAT1 genes. The other QTL was for directly measured C18:1 cis-9 at a Chr14: 1 .76 Mbp (rs208417762) locus, that was also in high LD (R 2 = 0.92) with missense mutations in the SLC52A2 and DGAT1 genes. Closer examination of association effects for FT-MIR predicted C16:0 revealed evidence of a peak at this locus, but the peak was marginally below the significance threshold. Notably, in the present study, the identified protein coding mutation in the SLC52A2 gene (rs134364612) was in perfect LD with the DGAT1 K232A polymorphism (rs109234250), which has been attributed to changes in bovine milk fat composition (Grisart et al., 2002;Fink et al., 2020) and  fatty acids (Bouwman et al., 2011;Buitenhuis et al., 2014;Li et al., 2014). The DGAT1 gene encodes diacylglycerol O-acyltransferase 1, an enzyme that catalyzes the final step in triglyceride production, thus making this a compelling candidate for the observed effects.
Two further QTL were identified for FT-MIR predicted C16:0 and C18:3 cis-3 at Chr27: 36 .20 Mbp loci that were not in high LD with a splice region variant, or a moderate or high impact coding variant (Table 4; Figures 3 and 4). However, the locus for C18:3 cis-3 (rs110950972) was in perfect LD with a 5′ untranslated region (rs208675276) in GPAT4, and the locus for C16:0 was also in high LD (R 2 = 0.997) with that same 5′ untranslated region. Interestingly, we found no evidence of QTL for the corresponding directly measured traits (Figures 3 and 4). The Chr27: 36 .20 Mbp loci are situated in the GPAT4 gene, which encodes the triglyceride synthesis enzyme glycerol-3-phosphate acyltransferase 4. As the milk fat percentage and other QTL at this locus have previously been shown to operate via a mechanism linked to gene expression (Littlejohn et al., 2014), it is not surprising that no significant coding mutations were identified in GPAT4.

Other Grouped Fatty Acid Effects.
Further significant effects were observed for directly measured SFA and UFA at a Chr19: 36 .19 Mbp locus (rs110980742), that was in high LD (R 2 > 0.81) with 2 missense mutations in the UTP18 gene (Table 3; Figure 5). This effect was not observed in any other individual or grouped fatty acid traits. The UTP18 gene is involved in the nucleolar processing of pre-18S ribosomal RNA, and has not previously been reported in GWAS of bovine milk composition. The signal at Chr19: 36 .19 is close to the locus of the KCNJ12 gene, which has a similar function to the KCNJ2 gene that has previously been shown to affect milk phenotypes , although a mechanism by which this gene could affect fatty acids is unclear.
Individual Milk Proteins. Significant effects were observed on BTA6, BTA11, BTA14 and BTA22 for individual milk proteins (Tables 3 and 4; Figure 6). Four QTL were identified on BTA6, 2 of which were for directly measured and FT-MIR predicted α-CN, and the other 2 for directly measured and FT-MIR predicted κ-CN, respectively. The effects for α-CN were observed at a Chr6: 87 .13 Mbp locus (rs109500363) that was in high LD (R 2 =0.92) with a missense mutation in the CSN1S1 gene (rs43703010). As the CSN1S1 gene codes for the α-CN protein (along with CSN1S2), it is not surprising that genetic signals affecting α-CN were enriched at this locus. Interestingly, FT-MIR predicted κ-CN also had a significant effect in the same region that was also in high LD (R 2 = 0.79) with rs43703010. The effect for directly measured κ-CN was observed at a Chr6: 87 .41 Mbp locus (rs110794953), which was in high LD (R 2 > 0.98) with 2 missense mutations in the CSN3 gene (rs43703015 and rs43703016). The CSN3 gene encodes κ-CN, an abundantly expressed milk protein. One of the missense mutations reported here (rs43703015) has previously been associated with milk composition traits and differential expression in mammary tissue (MacLeod et al., 2016). Significant effects have also been reported at this locus for a number of FT-MIR wavenumbers characterized by amide III and phosphate bands, C-H stretching vibrations of CH2 and -CH3, and N-H bending and C-N stretching in the amide II band .
Several QTL were identified for individual milk proteins on BTA11 that were in high LD (R 2 > 0.95) with missense mutations in the PAEP gene (rs110066229; rs109625649; Tables 3 and 4; Figure 6). Of these, the QTL with the most significant effects were observed for directly measured β-LG (P-value = 8.7e-117) and FT-MIR predicted β-LG (P-value = 5.4e-116). Smaller association effects were also observed for directly measured α-CN (P-value = 5.6e-10) and FT-MIR predicted β-CN (P-value = 8.3e-19). One of the implicated missense mutations in the PAEP gene was the V134A PAEP mutation (rs109625649) that distinguishes the 'A' and 'B' haplotypes of β-LG (previously described). This locus is likely driven by a regulatory effect, with a promoter variant reported to be in LD with the V134A mutation previously reported (Lum et al., 1997) to affect the binding of the Activator Protein-2 transcription factor. An expression QTL (eQTL) for PAEP was also reported in lactating bovine mammary tissue Davis et al., 2022).
One further QTL of interest was for directly measured Lf at a Chr22: 53 .54 Mbp locus (rs43765460; Table 3; Figure 6). The association effect at this locus had a P-value of 1.8e-41, but we found no relevant splice region variant, or moderate or high impact coding variant ascribed to this effect. However, the rs43765460 locus is a synonymous variant in the LTF gene. Using our previously published mammary RNA sequence dataset and eQTL mapping methodology (Lopdell et al., 2017;Tiplady et al., 2021), we identified the presence of a co-localized expression-based effect for LTF in this region. The rs43765460 locus we identified was in high LD with the top associated eQTL variant for Lf (R 2 = 0.88), and the Pearson correlation between the − log 10 (P-values) of the directly measured Lf QTL, and the −log 10 (P-values) of the Lf eQTL within a 1 Mbp region flanking the rs43765460 variant was 0.92. The LTF gene is a major iron-binding protein in milk that is linked to iron homeostasis and plays a key role in immune system response and cell growth. Previous studies have shown that the LTF gene is linked to changes in Lf concentrations in bovine milk (Bahar et al., 2011;Pawlik et al., 2014). Notably, there was no evidence of an association effect at or near this locus for FT-MIR predicted Lf (Table 4). Further, in a recent GWAS of individual FT-MIR wavenumbers, there was also no evidence of an association effect linked to the LTF gene , indicating that changes in milk composition due to this gene may not be easily detectable using FT-MIR spectral data. However, it is also important to note that prediction accuracies for Lf in the present study were relatively poor R cv 2 ( = 0.36; RPE cv = 0.19; Table 1), and the heritability estimate for FT-MIR predicted Lf was only 0.30, compared to the heritability estimate for directly measured Lf, which was 0.59 (Table 2). This pattern is similar to that which we observed for C14:1 and the SCD gene. That is, the component was in relatively low concentrations in the milk sample, model prediction accuracy  was relatively poor, the heritability for the measured trait was substantially higher than for the FT-MIR predicted trait, and a compelling peak was observed for the directly measured trait; however, no corresponding peak was observed for the FT-MIR predicted trait.

Perspectives on the Use of FT-MIR Trait Predictions in Dairy Cattle Selection
Utilizing FT-MIR predictions for fatty acids and proteins in milk can provide indicator traits across large numbers of animals at little or no marginal cost, because FT-MIR spectral data are already generated as part of routine milk testing to predict total fat and protein concentrations. The alternative to using FT-MIR trait predictions is to directly measure traits, which may be infeasible across even relatively small numbers of animals. Phenotypic correlations between directly measured and FT-MIR predicted traits provide a useful indication of the utility of FT-MIR trait predictions, particularly for herd management and milk processability traits. However, for breeding programs, we are also interested in the heritability of the FT-MIR predicted trait and the genetic correlation between the directly measured and FT-MIR predicted trait. This is because the heritability of the FT-MIR predicted trait defines the level of genetic variation present, whereas the genetic correlation between the directly measured and FT-MIR predicted trait defines the breeding progress we could expect in the directly measured trait if we were to select animals based on the FT-MIR predicted trait. Specifically, within the context of dairy cattle progeny test schemes, the genetic correlation will limit the relative amount of selection response that will result from using FT-MIR predictions instead of directly measured traits (Rutten et al., 2010). Based on this assumption, the genetic gain from selection using FT-MIR predictions for all traits we have studied would be greater than 70% of the gains achievable by direct selection on these traits; additionally, for 21 of the 29 traits, the genetic gains achievable would be greater than 85% of the gains achievable by direct selection.  However, it is important to note that this assumes that there is no true difference in heritability between the directly measured and FT-MIR predicted trait. For traits such as Lf and C14:1 where the estimated heritability of the direct measurement was lower than the heritability of the FT-MIR prediction, the genetic gain achievable would also be lower.
Although we observed high genetic correlations between directly measured and FT-MIR predicted traits in this study, the QTL underlying each trait were not always the same. An example of this includes where we observed a large association effect within the GPAT4 gene on BTA27 for FT-MIR predicted C18:3 cis-3, but no corresponding association effect was observed for directly measured C18:3 cis-3 (Figure 4). Similarly, a large association effect was observed for FT-MIR predicted β-CN within the PAEP gene on BTA11, but no corresponding association effect was observed in directly measured β-CN ( Figure 6). The presence of QTL with significant effects in an FT-MIR predicted trait only are not entirely surprising, given that FT-MIR predicted traits are a weighted linear function of absorbance values for individual wavenumbers, each of which may be underpinned by multiple genetic signals and QTL (Wang and Bovenhuis, 2018;Benedet et al., 2019;Zaalberg et al., 2020;Tiplady et al., 2021). This means that when a spectral wavenumber is included in a trait prediction equation, multiple genetic signals will also be present, some of which are specifically related to the trait of interest and some that are not. It is important that when FT-MIR predicted traits are used as proxies for other traits, we are mindful of this, particularly when using SNP-based approaches in our estimation of breeding values, whereby the impact will be determined by the relative proportion of genetic variation captured by each SNP and the interaction of additive effects between SNPs.
Instances also arose where a QTL was observed for a directly measured trait, but we found no corresponding QTL observed in the FT-MIR predicted trait. Examples of this include large association effects within the SCD gene for directly measured C10:1 and C14:1, but no corresponding association effects for individual FT-MIR predicted fatty acids (Figure 2). Similarly, a large association effect was observed within the LTF gene for directly measured Lf, but a corresponding association effect for FT-MIR predicted Lf was absent ( Figure 6).
In these examples, there was a consistent pattern where we have a component in relatively low concentrations in the milk sample, with relatively poor model prediction accuracies and lower heritability estimates for the FT-MIR predicted trait, compared with the directly measured trait (Tables 1 and 2). Although it might be argued that the failure to detect QTL in the SCD and LTF genes was because the calibration equations were inadequate for the task of quantifying the milk component targets (C10:1, C14:1, and Lf), it is also notable that in a previous GWAS we conducted on individual FT-MIR wavenumbers , no significant associations were identified between FT-MIR wavenumbers and variants within the SCD and LTF genes. Potentially, this means that changes in milk composition attributable to these 2 genes may be difficult to quantify directly using FT-MIR wavenumber spectra. For Lf to be detected using FT-MIR spectral data, it needs to provide a unique signal that distinguishes it from other whey proteins in solution that are at much higher concentrations. However, when the mean concentration of Lf is around 0.1 g/L and the major whey protein β-LG is at a 20-to 40-fold higher concentration, it is not surprising that a QTL is seen within the PAEP gene and not within the LTF gene.
With the growing interest in using FT-MIR spectral data to predict molecules at low concentrations in milk, it is important to understand that the predictive performance of these models may be limited, compared with models for predicting major milk components such as total fat and protein concentrations (Grelet et al., 2021). In the context of the present study, we have shown that for many fatty acids and protein traits, model prediction accuracies are moderate, but that genetic correlations between directly measured and FT-MIR predicted fatty acid and protein fractions are typically high. However, it is also clear that phenotypic variation between directly measured and FT-MIR predicted traits may be underpinned by differing genetic architecture. This may be due to several related factors including the trait of interest being at low concentrations in the milk sample, low prediction model accuracy, or that the trait is not easily detectable using FT-MIR spectroscopy. Improving calibration equations is central to optimizing our use of FT-MIR spectra to generate proxies for traits of interest to the industry such as fatty acids and protein fractions. Collaboration between research groups to generate data sets that include data from a range of herds that capture differences in climate, management systems, diet, and breed composition might improve calibration equations (Grelet et al., 2021). However, a key barrier to consolidating FT-MIR spectral data sets from different research groups is variation in spectral measurements between instruments and within instruments across time. Standardization of individual FT-MIR spectra wavenumbers using reference samples can effectively address these sources of variation (Grelet et al., 2015(Grelet et al., , 2017Tiplady et al., 2019); however, outside the European OptiMIR network, reference sample sharing and standardization is not common practice. Other approaches, such as those offered by Foss or Bentley, are appealing in that they are not reliant on perishable milk samples. However, as far as we are aware, the effectiveness of these procedures has not been independently evaluated. Validation of these within-instrument standardization procedures is important, because if the procedures work well, they could facilitate the consolidation of spectral data from different networks/ countries, and assist with the development of better prediction equations and improve trait prediction accuracies.

Study Limitations
In this study, we developed PLS prediction equations and compared the genetic characteristics of directly measured fatty acids and protein fractions to the same traits predicted from FT-MIR spectra. There are several areas of refinement that might improve prediction equations and the identification of QTL. First, before the development of prediction equations, we assessed several mathematical treatments of spectra, but we only assessed the prediction accuracies of those treatments using PLS models. Although PLS is a widely used method for developing calibration models from FT-MIR spectra, it may be possible to develop better prediction models for some traits by employing Bayesian or other machine learning approaches, as demonstrated in other studies of milk composition (Bonfatti et al., 2017a;El Jabri et al., 2019;Frizzarin et al., 2021), or animal health and feed intake traits (Dórea et al., 2018;Brand et al., 2021;Contla Hernández et al., 2021). Second, it is expected that increasing the number of samples in the study and including data from different herds would also improve trait prediction accuracies, particularly for fatty acids and protein fractions at low concentrations in milk samples. Extending the study to include data from different herds would also facilitate a more robust validation strategy. Although the cow-independent validation approach we have used is commonly practiced in studies of FT-MIR spectra trait prediction, it has been shown that record-or cow-independent validation can overinflate prediction accuracies, compared with herd-independent validation (Dórea et al., 2018;Lahart et al., 2019;Luke et al., 2019;Wang and Bovenhuis, 2019). Improving and validating the prediction equations we have developed in this study are important steps for future research to confirm their utility for prediction and use in future breeding programs.
Other potential refinements for the present study relate to genomic information and the strategy for identifying QTL. Specifically, we have used data sets mapped to the UMD3.1 genome; however, it is expected that the improved sequence continuity and per-base accu-racy of the ARS-UCD1.2 reference genome (Rosen et al., 2020) may yield a few additional QTL and reveal additional candidate mutations given improvements in accompanying transcript annotations. Also, the approach we used to identify QTL could be extended to account for nonadditive QTL in a similar manner to that outlined in Reynolds et al. (2021). Finally, the approach we used to identify causative genes and variants only considered protein-altering variants as candidates, which we acknowledge is relatively simple and crude, and that many of the identified signals could be underpinned by regulatory effects (e.g., gene expressionbased mechanisms). It is expected that integration of other functional data sets such as mammary eQTL and ChIP-seq data sets could map additional molecular QTL and enhance fine mapping and candidate variant identification .

CONCLUSIONS
We developed PLS calibration equations to predict bovine fatty acids and protein fractions in milk samples, and compared the genetic architecture underlying directly measured traits to that of corresponding FT-MIR predicted traits. Low to moderate prediction accuracies were observed, indicating that the potential application of using FT-MIR prediction equations for some traits may be limited. However, for most traits, heritability estimates were moderate to high, indicating that genetic variation exists that could potentially be exploited for the purposes of animal selection. Moreover, high genetic correlations between directly measured and FT-MIR predicted fatty acids and individual milk proteins indicated that selection based on FT-MIR predicted traits could provide high rates of genetic gain in the corresponding trait of interest. Trait QTL for fatty acids were identified with likely candidates in the DGAT1, CCDC57, SCD, and GPAT4 genes, but QTL underpinned by SCD were largely absent in FT-MIR predicted fatty acids. Similarly, likely candidates were identified for directly measured proteins in the CSN1S1, CSN3, PAEP, and LTF genes, but the QTL for CSN3 and LTF were absent in corresponding FT-MIR predicted traits. This highlighted that, in some instances, phenotypic variation for directly measured and FT-MIR predicted traits were underpinned by differing genetic architecture and segregation of alleles at QTL.
for the processing and analysis of milk samples, as well as the staff at Fonterra Research and Development Centre (Palmerston North, New Zealand) for milk analyses. Kathryn also thanks the wider LIC R&D team and fellow students for underlying technical support and thoughtful discussion, and Tracey Monehan (R&D Programme Manager, LIC) for overseeing the funding for this work. We also gratefully acknowledge the use of New Zealand eScience Infrastructure (NeSI) highperformance computing for this research. This research was funded through BoviQuest, a joint venture between LIC and ViaLactia Biosciences Ltd., a subsidiary (    1 Trait definitions and units as described in Table 1. 2 Untreated = untreated spectral data; first derivative = spectra pretreated with a first-order Savitzky-Golay derivative with a window of 7 data points either side; MSC = spectra pretreated with multiplicative scatter correction; MSC + first = spectra pretreated with MSC, followed by first-derivative transformation; SNV = spectra pretreated with a standard normal variate transformation; SNV + first = spectra pretreated with SNV followed by first-derivative transformation. 3 Cube-root transformation of lactoferrin (Lf).
Individual fatty acid (g/100 g of total fat) Trait definitions and units as described in Table 1. Standard errors shown in parentheses.