# Sequence-based GWAS and post-GWAS analyses reveal a key role of SLC37A1, ANKH, and regulatory regions on bovine milk mineral content

Marie-Pierre Sanchez, Dominique Rocha, Mathieu Charles, Mekki Boussaha, Chris Hozé, Mickaël Brochard, Agnes Delacroix-Buchet, Philippe Grosperrin, Didier Boichard

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marie-Pierre Sanchez, Dominique Rocha, Mathieu Charles, Mekki Boussaha, Chris Hozé, et al. Sequence-based GWAS and post-GWAS analyses reveal a key role of SLC37A1, ANKH, and regulatory regions on bovine milk mineral content. Scientific Reports, 2021, 11 (1), 10.1038/s41598-021-87078-1. hal-03191450

HAL Id: hal-03191450

https://hal.inrae.fr/hal-03191450

Submitted on 27 Feb 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Sequence-based GWAS and post-GWAS analyses reveal a key role of SLC37A1, ANKH, and regulatory regions on bovine milk mineral content

Marie-Pierre Sanchez, Dominique Rocha, Mathieu Charles, Mekki Boussaha, Chris Hozé, Mickaël Brochard, Agnes Delacroix-Buchet, Philippe Grosperrin, Didier Boichard

► To cite this version:

Marie-Pierre Sanchez, Dominique Rocha, Mathieu Charles, Mekki Boussaha, Chris Hozé, et al. Sequence-based GWAS and post-GWAS analyses reveal a key role of SLC37A1, ANKH, and regulatory regions on bovine milk mineral content. Scientific Reports, 2021, 11 (1), 10.1038/s41598-021-87078-1. hal-03191450

HAL Id: hal-03191450

https://hal.inrae.fr/hal-03191450v1

Submitted on 27 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

# Sequence-based GWAS and post-GWAS analyses reveal a key role of SLC37A1, ANKH, and regulatory regions on bovine milk mineral content

Marie-Pierre Sanchez<fnref n="1" />✉, Dominique Rocha<fnref n="1" />, Mathieu Charles<fnref n="1" />, Mekki Boussaha<fnref n="1" />, Chris Hozé<fnref n="1" />,², Mickaël Brochard³, Agnès Delacroix-Buchet<fnref n="1" />, Philippe Grosperrin⁴ & Didier Boichard<fnref n="1" />

The mineral composition of bovine milk plays an important role in determining its nutritional and cheese-making value. Concentrations of the main minerals predicted from mid-infrared spectra produced during milk recording, combined with cow genotypes, provide a unique opportunity to decipher the genetic determinism of these traits. The present study included 1 million test-day predictions of Ca, Mg, P, K, Na, and citrate content from 126,876 Montbéliarde cows, of which 19,586 had genotype data available. All investigated traits were highly heritable (0.50–0.58), with the exception of Na (0.32). A sequence-based genome-wide association study (GWAS) detected 50 QTL (18 affecting two to five traits) and positional candidate genes and variants, mostly located in non-coding sequences. In silico post-GWAS analyses highlighted 877 variants that could be regulatory SNPs altering transcription factor (TF) binding sites or located in non-coding RNA (mainly lncRNA). Furthermore, we found 47 positional candidate genes and 45 TFs highly expressed in mammary gland compared to 90 other bovine tissues. Among the mammary-specific genes, SLC37A1 and ANKH, encoding proteins involved in ion transport were located in the most significant QTL. This study therefore highlights a comprehensive set of functional candidate genes and variants that affect milk mineral content.

Bovine milk contains many essential nutrients, such as lactose (~ 48 g/L), fatty acids (~ 37 g/L), proteins (~ 34 g/L), and minerals (~ 9 g/L). Although less abundant than other solid components of milk, the major minerals—potassium (K), calcium (Ca), phosphorus (P), sodium (Na), and magnesium (Mg)—have an important effect on both human health and the cheese-making process. In humans, all of these minerals are necessary for many vital functions and therefore for the maintenance of good health. Dairy products can represent an important source of minerals in the human diet, especially of well-assimilated Ca<fnref n="1" />. In milk, minerals are found either in solution (soluble fraction) or in colloidal form (insoluble fraction). Some minerals are exclusively found in the soluble fraction (e.g., K and Na) while others exist in both fractions (e.g., Ca, P, and Mg). In the soluble fraction, Ca, P, and Mg exist in different forms, including ions and salts (phosphates and citrates), while in colloidal form, they are associated with casein molecules in the micelles and play a role in the structure and stability of these assemblages during the cheesemaking process². Higher mineral concentrations are therefore associated with improved coagulation properties of milk³,⁴ and could enhance the nutritional value for human consumers.

Despite the potential benefits to human nutrition and milk processing, little is known about the genetic factors that influence milk mineral composition, mainly because the determination of mineral content via reference analyses is costly and time-consuming. A number of studies have reported genetic variation in milk mineral composition⁵⁻¹¹, but to our knowledge, only two studies have conducted genome-wide association studies (GWAS) to investigate the genomic regions associated with these traits⁹,¹². Both of those prior studies used

<footnote n="1">¹Université Paris-Saclay, INRAE, AgroParisTech, GABI, 78350 Jouy-en-Josas, France. ²Allice, 75012 Paris, France. ³Umotest, 01250 Ceyzériat, France. ⁴Conseil Elevage 25-90, 25640 Roulans, France. ⁵email: marie-pierre.sanchez@inrae.fr</footnote>
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Sequence-based GWAS and post-GWAS analyses reveal a key role of SLC37A1, ANKH, and regulatory regions on bovine milk mineral content

Marie-Pierre Sanchez<fnref n="1" />✉, Dominique Rocha<fnref n="1" />, Mathieu Charles<fnref n="1" />, Mekki Boussaha<fnref n="1" />, Chris Hozé<fnref n="1" />,², Mickaël Brochard³, Agnès Delacroix-Buchet<fnref n="1" />, Philippe Grosperrin⁴ & Didier Boichard<fnref n="1" />

The mineral composition of bovine milk plays an important role in determining its nutritional and cheese-making value. Concentrations of the main minerals predicted from mid-infrared spectra produced during milk recording, combined with cow genotypes, provide a unique opportunity to decipher the genetic determinism of these traits. The present study included 1 million test-day predictions of Ca, Mg, P, K, Na, and citrate content from 126,876 Montbéliarde cows, of which 19,586 had genotype data available. All investigated traits were highly heritable (0.50–0.58), with the exception of Na (0.32). A sequence-based genome-wide association study (GWAS) detected 50 QTL (18 affecting two to five traits) and positional candidate genes and variants, mostly located in non-coding sequences. In silico post-GWAS analyses highlighted 877 variants that could be regulatory SNPs altering transcription factor (TF) binding sites or located in non-coding RNA (mainly lncRNA). Furthermore, we found 47 positional candidate genes and 45 TFs highly expressed in mammary gland compared to 90 other bovine tissues. Among the mammary-specific genes, SLC37A1 and ANKH, encoding proteins involved in ion transport were located in the most significant QTL. This study therefore highlights a comprehensive set of functional candidate genes and variants that affect milk mineral content.

Bovine milk contains many essential nutrients, such as lactose (~ 48 g/L), fatty acids (~ 37 g/L), proteins (~ 34 g/L), and minerals (~ 9 g/L). Although less abundant than other solid components of milk, the major minerals—potassium (K), calcium (Ca), phosphorus (P), sodium (Na), and magnesium (Mg)—have an important effect on both human health and the cheese-making process. In humans, all of these minerals are necessary for many vital functions and therefore for the maintenance of good health. Dairy products can represent an important source of minerals in the human diet, especially of well-assimilated Ca<fnref n="1" />. In milk, minerals are found either in solution (soluble fraction) or in colloidal form (insoluble fraction). Some minerals are exclusively found in the soluble fraction (e.g., K and Na) while others exist in both fractions (e.g., Ca, P, and Mg). In the soluble fraction, Ca, P, and Mg exist in different forms, including ions and salts (phosphates and citrates), while in colloidal form, they are associated with casein molecules in the micelles and play a role in the structure and stability of these assemblages during the cheesemaking process². Higher mineral concentrations are therefore associated with improved coagulation properties of milk³,⁴ and could enhance the nutritional value for human consumers.

Despite the potential benefits to human nutrition and milk processing, little is known about the genetic factors that influence milk mineral composition, mainly because the determination of mineral content via reference analyses is costly and time-consuming. A number of studies have reported genetic variation in milk mineral composition⁵⁻¹¹, but to our knowledge, only two studies have conducted genome-wide association studies (GWAS) to investigate the genomic regions associated with these traits⁹,¹². Both of those prior studies used

<footnote n="1">¹Université Paris-Saclay, INRAE, AgroParisTech, GABI, 78350 Jouy-en-Josas, France. ²Allice, 75012 Paris, France. ³Umotest, 01250 Ceyzériat, France. ⁴Conseil Elevage 25-90, 25640 Roulans, France. ⁵email: marie-pierre.sanchez@inrae.fr</footnote>
<table><caption>**Table 1.** Mean, standard deviation (SD), and accuracy (R²), estimated by cross-validation, of mid-infrared (MIR) predictions for concentrations of minerals and citrate in milk from Montbéliarde cows.</caption><thead><tr><th>Trait</th><th>Abbrev</th><th>Mean</th><th>SD</th><th>R<sup>2</sup> MIR prediction</th></tr></thead><tbody><tr><td>Calcium, in mg kg<sup>-1</sup> of milk</td><td>Ca</td><td>1161.4</td><td>92.6</td><td>0.82</td></tr><tr><td>Phosphorus, in mg kg<sup>-1</sup> of milk</td><td>P</td><td>1007.0</td><td>77.5</td><td>0.75</td></tr><tr><td>Magnesium, in mg kg<sup>-1</sup> of milk</td><td>Mg</td><td>1473.8</td><td>104.5</td><td>0.77</td></tr><tr><td>Potassium, in mg kg<sup>-1</sup> of milk</td><td>K</td><td>100.5</td><td>7.2</td><td>0.68</td></tr><tr><td>Sodium, in mg kg<sup>-1</sup> of milk</td><td>Na</td><td>341.7</td><td>44.5</td><td>0.44</td></tr><tr><td>Citrate, in g kg<sup>-1</sup> of milk</td><td>Citrate</td><td>8.27</td><td>1.49</td><td>0.90</td></tr></tbody></table>

<table><caption>**Table 2.** Estimates of heritability (in bold, diagonal) and genetic correlation (above the diagonal) for the concentrations of minerals and citrate in milk (SE < 0.01).</caption><thead><tr><th></th><th>Ca</th><th>P</th><th>Mg</th><th>K</th><th>Na</th><th>Citrate</th></tr></thead><tbody><tr><td>Ca</td><td><strong>0.50</strong></td><td>0.34</td><td>0.60</td><td>-0.22</td><td>-0.20</td><td>0.57</td></tr><tr><td>P</td><td></td><td><strong>0.56</strong></td><td>0.58</td><td>0.39</td><td>-0.23</td><td>-0.16</td></tr><tr><td>Mg</td><td></td><td></td><td><strong>0.52</strong></td><td>0.11</td><td>-0.02</td><td>0.59</td></tr><tr><td>K</td><td></td><td></td><td></td><td><strong>0.53</strong></td><td>-0.05</td><td>0.01</td></tr><tr><td>Na</td><td></td><td></td><td></td><td></td><td><strong>0.32</strong></td><td>-0.15</td></tr><tr><td>Citrate</td><td></td><td></td><td></td><td></td><td></td><td><strong>0.48</strong></td></tr></tbody></table>

777 k SNPs and a relatively small sample of cows. As an alternative to reference analyses, mid-infrared (MIR) spectrometry can predict various milk components, including mineral fractions¹³⁻¹⁵, quickly and cheaply. Because of these advantages, milk MIR spectra are routinely recorded and stored. The combination of this technology with i) the widespread genotyping of cows for genomic selection, ii) the availability of whole-genome sequence (WGS) data from the 1000 Bull Genomes Project¹⁶, and iii) ever-increasing knowledge of the bovine genome¹⁷,¹⁸ creates the possibility of large-scale, high-resolution analyses for identification of the genes and variants that affect the mineral content of milk. We have previously applied this approach—whole-genome sequence-based GWAS combined with MIR predictions—to investigate milk protein composition¹⁹ and cheese-making traits²⁰, and in both cases we succeeded in highlighting functional candidate genes. In particular, the genes *SLC37A1* and *ANKH* were strongly linked with milk quality; both encode transmembrane proteins involved in ion transport and are therefore likely to have an effect on milk mineral composition. However, in these genes, as well as in other genes we identified, the variants with the most significant effects were mostly found in non-coding regions with limited annotation, which made it difficult to distinguish the causal variant. This pattern is quite general and many studies have reported the major role of regulatory variants in the architecture of complex traits²¹. To address this challenge, further investigation of non-coding regions is needed, particularly with respect to binding sites of transcription factors and non-coding RNA which could regulate gene expression.

The main objective of this study was to identify the best candidate genes and variants that might affect the content of Ca, P, Mg, K, Na, and citrate in milk, as predicted from MIR spectra. For this, we first conducted a GWAS on imputed WGS data of 19,586 Montbéliarde cows, and then performed post-GWAS analyses using different sources of annotation data to further refine our results.

## Results

We analyzed six traits predicted from MIR spectra in Montbéliarde cows, representing the mineral (Ca, P, Mg, K, and Na) and citrate content of milk. MIR prediction equations originated from the Optimir project¹⁴,¹⁵. The accuracies of these MIR predictions, as assessed by the coefficient of determination (R²) in a validation population (Table 1), ranged from 0.68 to 0.90, with the exception of Na (0.44).

**Heritability and genetic correlation estimates for mineral and citrate content.** Genetic parameters, i.e. heritabilities (h²) and genetic correlations (r_g), were estimated for milk mineral and citrate content, as predicted from more than 1 million test-day records from 126,873 cows (Table 2). At the test-day level, heritability estimates were moderate for Na content (h²=0.32) but higher for other minerals (h²=0.50 to 0.56) and citrate (h²=0.48). Na was negatively and poorly correlated with other minerals (-0.23 ≤ r_g ≤ -0.02) and citrate (r_g = -0.15), while the levels of most other minerals were generally positively correlated (0.11 ≤ r_g ≤ 0.60); the exception was the relationship between Ca and K (r_g = -0.22). Values of the genetic correlation between citrate and mineral levels in milk ranged quite broadly, depending on the mineral: null with K (r_g = 0.01 with K), slightly negative with Na (r_g = -0.15) and P (r_g = -0.16), and highly positive with Ca (r_g = 0.57) and Mg (r_g = 0.59).
<table><caption>**Table 1.** Mean, standard deviation (SD), and accuracy (R²), estimated by cross-validation, of mid-infrared (MIR) predictions for concentrations of minerals and citrate in milk from Montbéliarde cows.</caption><thead><tr><th>Trait</th><th>Abbrev</th><th>Mean</th><th>SD</th><th>R<sup>2</sup> MIR prediction</th></tr></thead><tbody><tr><td>Calcium, in mg kg<sup>-1</sup> of milk</td><td>Ca</td><td>1161.4</td><td>92.6</td><td>0.82</td></tr><tr><td>Phosphorus, in mg kg<sup>-1</sup> of milk</td><td>P</td><td>1007.0</td><td>77.5</td><td>0.75</td></tr><tr><td>Magnesium, in mg kg<sup>-1</sup> of milk</td><td>Mg</td><td>1473.8</td><td>104.5</td><td>0.77</td></tr><tr><td>Potassium, in mg kg<sup>-1</sup> of milk</td><td>K</td><td>100.5</td><td>7.2</td><td>0.68</td></tr><tr><td>Sodium, in mg kg<sup>-1</sup> of milk</td><td>Na</td><td>341.7</td><td>44.5</td><td>0.44</td></tr><tr><td>Citrate, in g kg<sup>-1</sup> of milk</td><td>Citrate</td><td>8.27</td><td>1.49</td><td>0.90</td></tr></tbody></table>

<table><caption>**Table 2.** Estimates of heritability (in bold, diagonal) and genetic correlation (above the diagonal) for the concentrations of minerals and citrate in milk (SE < 0.01).</caption><thead><tr><th></th><th>Ca</th><th>P</th><th>Mg</th><th>K</th><th>Na</th><th>Citrate</th></tr></thead><tbody><tr><td>Ca</td><td><strong>0.50</strong></td><td>0.34</td><td>0.60</td><td>-0.22</td><td>-0.20</td><td>0.57</td></tr><tr><td>P</td><td></td><td><strong>0.56</strong></td><td>0.58</td><td>0.39</td><td>-0.23</td><td>-0.16</td></tr><tr><td>Mg</td><td></td><td></td><td><strong>0.52</strong></td><td>0.11</td><td>-0.02</td><td>0.59</td></tr><tr><td>K</td><td></td><td></td><td></td><td><strong>0.53</strong></td><td>-0.05</td><td>0.01</td></tr><tr><td>Na</td><td></td><td></td><td></td><td></td><td><strong>0.32</strong></td><td>-0.15</td></tr><tr><td>Citrate</td><td></td><td></td><td></td><td></td><td></td><td><strong>0.48</strong></td></tr></tbody></table>

777 k SNPs and a relatively small sample of cows. As an alternative to reference analyses, mid-infrared (MIR) spectrometry can predict various milk components, including mineral fractions¹³⁻¹⁵, quickly and cheaply. Because of these advantages, milk MIR spectra are routinely recorded and stored. The combination of this technology with i) the widespread genotyping of cows for genomic selection, ii) the availability of whole-genome sequence (WGS) data from the 1000 Bull Genomes Project¹⁶, and iii) ever-increasing knowledge of the bovine genome¹⁷,¹⁸ creates the possibility of large-scale, high-resolution analyses for identification of the genes and variants that affect the mineral content of milk. We have previously applied this approach—whole-genome sequence-based GWAS combined with MIR predictions—to investigate milk protein composition¹⁹ and cheese-making traits²⁰, and in both cases we succeeded in highlighting functional candidate genes. In particular, the genes SLC37A1 and ANKH were strongly linked with milk quality; both encode transmembrane proteins involved in ion transport and are therefore likely to have an effect on milk mineral composition. However, in these genes, as well as in other genes we identified, the variants with the most significant effects were mostly found in non-coding regions with limited annotation, which made it difficult to distinguish the causal variant. This pattern is quite general and many studies have reported the major role of regulatory variants in the architecture of complex traits²¹. To address this challenge, further investigation of non-coding regions is needed, particularly with respect to binding sites of transcription factors and non-coding RNA which could regulate gene expression.

The main objective of this study was to identify the best candidate genes and variants that might affect the content of Ca, P, Mg, K, Na, and citrate in milk, as predicted from MIR spectra. For this, we first conducted a GWAS on imputed WGS data of 19,586 Montbéliarde cows, and then performed post-GWAS analyses using different sources of annotation data to further refine our results.

## Results

We analyzed six traits predicted from MIR spectra in Montbéliarde cows, representing the mineral (Ca, P, Mg, K, and Na) and citrate content of milk. MIR prediction equations originated from the Optimir project¹⁴,¹⁵. The accuracies of these MIR predictions, as assessed by the coefficient of determination (R²) in a validation population (Table 1), ranged from 0.68 to 0.90, with the exception of Na (0.44).

**Heritability and genetic correlation estimates for mineral and citrate content.** Genetic parameters, i.e. heritabilities (h²) and genetic correlations (r_g), were estimated for milk mineral and citrate content, as predicted from more than 1 million test-day records from 126,873 cows (Table 2). At the test-day level, heritability estimates were moderate for Na content (h²=0.32) but higher for other minerals (h²=0.50 to 0.56) and citrate (h²=0.48). Na was negatively and poorly correlated with other minerals (-0.23 ≤ r_g ≤ -0.02) and citrate (r_g = -0.15), while the levels of most other minerals were generally positively correlated (0.11 ≤ r_g ≤ 0.60); the exception was the relationship between Ca and K (r_g = -0.22). Values of the genetic correlation between citrate and mineral levels in milk ranged quite broadly, depending on the mineral: null with K (r_g = 0.01 with K), slightly negative with Na (r_g = -0.15) and P (r_g = -0.16), and highly positive with Ca (r_g = 0.57) and Mg (r_g = 0.59).
<figure><img src="image_4.png" /><figcaption>**Figure 1.** –log₁₀(*P*) value of the effect of variants on milk mineral (Ca, P, Mg, K, and Na) or citrate content plotted against their position on *Bos taurus* autosomes.</figcaption></figure>

**QTL identified in GWAS.** We conducted single-trait GWASs on imputed WGSs from 19,586 cows (12,907,802 variants) and primarily identified 96 trait x region combinations with significant effects (-log₁₀(P) ≥ 8.4) on milk mineral and citrate composition by applying the procedure described in the Methods section (Fig. 1). In two genomic regions (a 12Mbp-region on BTA1 and a 27Mbp-region on BTA20), the single marker analyses identified multiple trait x region combinations (10 and 15, respectively). Therefore, conditional analyses including the most significant variant identified in each region was applied to decipher if significant
<figure><img src="image_3.png" /><figcaption>**Figure 1.** –log₁₀(*P*) value of the effect of variants on milk mineral (Ca, P, Mg, K, and Na) or citrate content plotted against their position on *Bos taurus* autosomes.</figcaption></figure>

**QTL identified in GWAS.** We conducted single-trait GWASs on imputed WGSs from 19,586 cows (12,907,802 variants) and primarily identified 96 trait x region combinations with significant effects (-log₁₀(P) ≥ 8.4) on milk mineral and citrate composition by applying the procedure described in the Methods section (Fig. 1). In two genomic regions (a 12Mbp-region on BTA1 and a 27Mbp-region on BTA20), the single marker analyses identified multiple trait x region combinations (10 and 15, respectively). Therefore, conditional analyses including the most significant variant identified in each region was applied to decipher if significant
<table><caption>**Table 3.** QTL identified for milk mineral and citrate content on *Bos taurus* (BTA) autosomes 1 to 10 (when a QTL region affected multiple traits, the most significantly affected trait is indicated in bold).</caption><thead><tr><th rowspan="2">Region</th><th rowspan="2">QTL</th><th rowspan="2">Trait</th><th rowspan="2">BTA</th><th colspan="2">Confidence Interval of the QTL</th><th rowspan="2"># variants</th><th rowspan="2"># genes</th><th rowspan="2">Position (bp)</th><th rowspan="2">ID</th><th rowspan="2">p</th><th rowspan="2">R<sup>2</sup></th><th rowspan="2">-log<sub>10</sub>P</th><th rowspan="2">b</th><th rowspan="2">SE</th><th rowspan="2">Name of gene(s)</th><th rowspan="2">Impact</th></tr><tr><th>From (bp)</th><th>To (bp)</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>P</td><td>1</td><td>135,384,680</td><td>136,262,601</td><td>8</td><td>3</td><td>136,262,601</td><td>rs210630730</td><td>0.199</td><td>0.69</td><td>11.5</td><td>8.57</td><td>1.23</td><td>TMEM108</td><td>Intron</td></tr><tr><td>2</td><td>4</td><td>Mg</td><td>1</td><td>142,826,156</td><td>142,837,365</td><td>9</td><td>1</td><td>142,834,737</td><td>rs109459130</td><td>0.576</td><td>0.85</td><td>64.1</td><td>1.48</td><td>0.09</td><td>SLC37A1</td><td>Intron</td></tr><tr><td>2</td><td>5</td><td>P</td><td>1</td><td><strong>142,826,156</strong></td><td><strong>142,835,551</strong></td><td>7</td><td><strong>1</strong></td><td><strong>142,834,737</strong></td><td><strong>rs109459130</strong></td><td><strong>0.576</strong></td><td><strong>0.85</strong></td><td><strong>219.3</strong></td><td><strong>30.2</strong></td><td><strong>0.95</strong></td><td><strong>SLC37A1</strong></td><td><strong>Intron</strong></td></tr><tr><td>2</td><td>6</td><td>K</td><td>1</td><td>142,826,156</td><td>142,838,477</td><td>31</td><td>1</td><td>142,835,551</td><td>rs109717634</td><td>0.613</td><td>1.00</td><td>138.1</td><td>22.4</td><td>0.89</td><td>SLC37A1</td><td>Intron</td></tr><tr><td>2</td><td>7</td><td>Na</td><td>1</td><td>142,826,156</td><td>142,835,551</td><td>7</td><td>1</td><td>142,835,551</td><td>rs109717634</td><td>0.613</td><td>1.00</td><td>19.5</td><td>-3.77</td><td>0.41</td><td>SLC37A1</td><td>Intron</td></tr><tr><td>3</td><td>9</td><td>Mg</td><td>1</td><td>145,114,056</td><td>146,804,023</td><td>64</td><td>3</td><td>146,775,338</td><td>rs43048006</td><td>0.630</td><td>0.33</td><td>10.3</td><td>-0.70</td><td>0.11</td><td>-</td><td>-</td></tr><tr><td>4</td><td>11</td><td>Na</td><td>2</td><td>5,738,340</td><td>7,483,560</td><td>104</td><td>7</td><td>5,739,439</td><td>rs209271012</td><td>0.111</td><td>0.86</td><td>10.4</td><td>3.85</td><td>0.58</td><td>LOC104971101</td><td>IncRNA</td></tr><tr><td>5</td><td>12</td><td>Na</td><td>2</td><td>131,111,763</td><td>131,742,987</td><td>206</td><td>9</td><td>131,197,793</td><td>rs133533873</td><td>0.672</td><td>0.99</td><td>13.8</td><td>2.84</td><td>0.37</td><td>ALPL</td><td>Intron</td></tr><tr><td>6</td><td>13</td><td>Na</td><td>3</td><td>15,459,107</td><td>15,540,579</td><td>7</td><td>6</td><td>15,470,670</td><td>rs110073735</td><td>0.026</td><td>0.92</td><td>8.9</td><td>-5.20</td><td>0.86</td><td>EFNA1 /<br>LOC107132270</td><td>Intron /<br>IncRNA</td></tr><tr><td>7</td><td>14</td><td>Ca</td><td>4</td><td>48,711,128</td><td>48,835,496</td><td>22</td><td>4</td><td>48,711,294</td><td>rs381111559</td><td>0.291</td><td>0.59</td><td>13.0</td><td>-10.5</td><td>1.42</td><td>SLC26A4</td><td>Intron</td></tr><tr><td>8</td><td>15</td><td>Na</td><td>4</td><td>75,153,659</td><td>77,188,370</td><td>14</td><td>7</td><td>77,147,308</td><td>rs386097640</td><td>0.024</td><td>0.87</td><td>9.5</td><td>6.65</td><td>1.06</td><td>GCK</td><td>Intron</td></tr><tr><td>8</td><td>16</td><td>Ca</td><td>4</td><td>77,147,308</td><td>77,257,775</td><td>9</td><td>7</td><td>77,173,366</td><td>rs133824313</td><td>0.025</td><td>0.98</td><td>9.1</td><td>14.8</td><td>2.40</td><td>GCK</td><td>Intron</td></tr><tr><td>8</td><td>17</td><td>Mg</td><td>4</td><td><strong>77,037,103</strong></td><td><strong>77,188,370</strong></td><td><strong>8</strong></td><td><strong>4</strong></td><td><strong>77,173,366</strong></td><td><strong>rs133824313</strong></td><td><strong>0.025</strong></td><td><strong>0.98</strong></td><td><strong>11.1</strong></td><td><strong>1.34</strong></td><td><strong>0.20</strong></td><td><strong>GCK</strong></td><td><strong>Intron</strong></td></tr><tr><td>9</td><td>18</td><td>Ca</td><td>4</td><td>91,807,789</td><td>92,192,288</td><td>182</td><td>6</td><td>92,192,288</td><td>rs42451422</td><td>0.028</td><td>0.95</td><td>9.5</td><td>13.4</td><td>2.13</td><td>LRRC4 / SND1</td><td>Intron /<br>Intron</td></tr><tr><td>10</td><td>19</td><td>Ca</td><td>5</td><td>108,256,181</td><td>109,217,271</td><td>87</td><td>11</td><td>109,141,258</td><td>rs383141568</td><td>0.570</td><td>0.82</td><td>8.4</td><td>5.73</td><td>0.98</td><td>BID</td><td>Intron</td></tr><tr><td>11</td><td>20</td><td>Ca</td><td>5</td><td>111,637,039</td><td>111,848,859</td><td>47</td><td>7</td><td>111,745,513</td><td>rs384814366</td><td>0.231</td><td>0.94</td><td>16.6</td><td>-9.38</td><td>1.11</td><td>MRTFA = MKL1</td><td>Intron</td></tr><tr><td>12</td><td>21</td><td>Ca</td><td>5</td><td><strong>116,403,262</strong></td><td><strong>116,463,104</strong></td><td><strong>48</strong></td><td><strong>8</strong></td><td><strong>116,429,840</strong></td><td><strong>rs382268272</strong></td><td><strong>0.049</strong></td><td><strong>0.58</strong></td><td><strong>60.5</strong></td><td><strong>-34.8</strong></td><td><strong>2.11</strong></td><td><strong>LOC101903383</strong></td><td><strong>IncRNA</strong></td></tr><tr><td>12</td><td>22</td><td>Na</td><td>5</td><td>116,403,262</td><td>116,907,862</td><td>72</td><td>17</td><td>116,429,922</td><td>rs110716509</td><td>0.050</td><td>0.67</td><td>9.3</td><td>5.05</td><td>0.81</td><td>LOC101903383</td><td>IncRNA</td></tr><tr><td>12</td><td>23</td><td>Citrate</td><td>5</td><td>116,403,262</td><td>116,468,630</td><td>50</td><td>8</td><td>116,439,760</td><td>rs110048000</td><td>0.042</td><td>0.56</td><td>9.5</td><td>-0.25</td><td>0.04</td><td>PPARA</td><td>Intron</td></tr><tr><td>12</td><td>24</td><td>Mg</td><td>5</td><td>116,403,806</td><td>116,463,104</td><td>36</td><td>6</td><td>116,448,403</td><td>rs110659119</td><td>0.048</td><td>0.97</td><td>24.2</td><td>-1.38</td><td>0.13</td><td>PPARA</td><td>Intron</td></tr><tr><td>12</td><td>25</td><td>P</td><td>5</td><td>116,420,172</td><td>117,323,874</td><td>28</td><td>5</td><td>116,448,403</td><td>rs110659119</td><td>0.048</td><td>0.97</td><td>19.5</td><td>-13.5</td><td>1.46</td><td>PPARA</td><td>Intron</td></tr><tr><td>13</td><td>26</td><td>Ca</td><td>5</td><td><strong>119,066,644</strong></td><td><strong>119,069,138</strong></td><td><strong>25</strong></td><td><strong>1</strong></td><td><strong>119,068,337</strong></td><td><strong>rs134583661</strong></td><td><strong>0.046</strong></td><td><strong>0.95</strong></td><td><strong>48.4</strong></td><td><strong>-24.8</strong></td><td><strong>1.69</strong></td><td>-</td><td><strong>Intergenic</strong></td></tr><tr><td>13</td><td>27</td><td>Mg</td><td>5</td><td>119,066,644</td><td>119,069,138</td><td>25</td><td>1</td><td>119,068,337</td><td>rs134583661</td><td>0.046</td><td>0.95</td><td>20.1</td><td>-1.29</td><td>0.14</td><td>-</td><td>Intergenic</td></tr><tr><td>13</td><td>28</td><td>P</td><td>5</td><td>119,066,644</td><td>119,069,138</td><td>25</td><td>1</td><td>119,068,337</td><td>rs134583661</td><td>0.046</td><td>0.95</td><td>17.9</td><td>-13.2</td><td>1.50</td><td>-</td><td>Intergenic</td></tr><tr><td>14</td><td>29</td><td>K</td><td>6</td><td>40,012,262</td><td>40,518,049</td><td>18</td><td>2</td><td>40,024,540</td><td>rs43610452</td><td>0.064</td><td>0.99</td><td>11.7</td><td>-11.6</td><td>1.66</td><td>SLIT2</td><td>Intron</td></tr><tr><td>15</td><td>30</td><td>P</td><td>6</td><td>42,855,236</td><td>43,032,274</td><td>17</td><td>3</td><td>42,870,371</td><td>rs42610629</td><td>0.523</td><td>0.97</td><td>9.2</td><td>5.59</td><td>0.91</td><td>-</td><td>Intergenic</td></tr><tr><td>16</td><td>31</td><td>K</td><td>6</td><td>43,297,568</td><td>43,322,676</td><td>12</td><td>1</td><td>43,322,676</td><td>rs43463132</td><td>0.778</td><td>0.97</td><td>40.7</td><td>13.2</td><td>0.98</td><td>-</td><td>Intergenic</td></tr><tr><td>17</td><td>32</td><td>P</td><td>6</td><td>45,018,531</td><td>45,401,570</td><td>78</td><td>4</td><td>45,056,055</td><td>rs378008740</td><td>0.713</td><td>0.76</td><td>30.6</td><td>12.7</td><td>1.09</td><td>-</td><td>Intergenic</td></tr><tr><td>17</td><td>33</td><td>K</td><td>6</td><td><strong>45,018,882</strong></td><td><strong>45,401,570</strong></td><td><strong>66</strong></td><td><strong>4</strong></td><td><strong>45,328,581</strong></td><td><strong>rs135848932</strong></td><td><strong>0.174</strong></td><td><strong>1.00</strong></td><td><strong>120.0</strong></td><td><strong>-27.4</strong></td><td><strong>1.17</strong></td><td><strong>SEL1L3</strong></td><td><strong>Intron</strong></td></tr><tr><td>17</td><td>34</td><td>Ca</td><td>6</td><td>45,328,392</td><td>45,411,116</td><td>17</td><td>4</td><td>45,393,335</td><td>rs134399380</td><td>0.183</td><td>0.91</td><td>8.6</td><td>7.24</td><td>1.22</td><td>-</td><td>Intergenic</td></tr><tr><td>17</td><td>35</td><td>Na</td><td>6</td><td>45,215,035</td><td>45,401,570</td><td>17</td><td>4</td><td>45,401,485</td><td>rs108972810</td><td>0.177</td><td>0.99</td><td>10.2</td><td>-3.39</td><td>0.52</td><td>SMIM20 /<br>LOC112447060</td><td>Upstream<br>BS DOF5.3 /<br>IncRNA</td></tr><tr><td>18</td><td>36</td><td>K</td><td>6</td><td>47,687,737</td><td>49,094,600</td><td>79</td><td>1</td><td>48,145,964</td><td>rs43461552</td><td>0.336</td><td>0.92</td><td>26.9</td><td>-10.6</td><td>0.98</td><td>-</td><td>Intergenic</td></tr><tr><td>19</td><td>37</td><td>K</td><td>6</td><td>50,277,420</td><td>51,315,936</td><td>82</td><td>2</td><td>51,311,430</td><td>rs110769289</td><td>0.555</td><td>0.53</td><td>9.1</td><td>-5.72</td><td>0.93</td><td>-</td><td>Intergenic</td></tr><tr><td>20</td><td>38</td><td>K</td><td>6</td><td>56,960,135</td><td>56,970,541</td><td>6</td><td>1</td><td>56,968,022</td><td>rs134533189</td><td>0.091</td><td>0.41</td><td>9.4</td><td>-8.37</td><td>1.34</td><td>RELL1</td><td>Intron</td></tr><tr><td>21</td><td>39</td><td>Na</td><td>6</td><td>85,409,619</td><td>86,104,900</td><td>312</td><td>15</td><td>85,512,651</td><td>rs110772556</td><td>0.660</td><td>0.98</td><td>8.6</td><td>2.79</td><td>0.47</td><td>-</td><td>Intergenic</td></tr><tr><td>21</td><td>40</td><td>Mg</td><td>6</td><td><strong>85,587,260</strong></td><td><strong>85,599,456</strong></td><td><strong>41</strong></td><td><strong>2</strong></td><td><strong>85,594,649</strong></td><td><strong>rs209134503</strong></td><td><strong>0.481</strong></td><td><strong>0.85</strong></td><td><strong>29.8</strong></td><td><strong>-1.09</strong></td><td><strong>0.10</strong></td><td><strong>ODAM</strong></td><td><strong>Synonymous</strong></td></tr><tr><td>21</td><td>41</td><td>Ca</td><td>6</td><td>85,406,891</td><td>85,599,456</td><td>44</td><td>5</td><td>85,594,717</td><td>rs135113545</td><td>0.485</td><td>0.79</td><td>17.0</td><td>-10.6</td><td>1.24</td><td>ODAM</td><td>Intron</td></tr><tr><td>21</td><td>42</td><td>P</td><td>6</td><td>85,564,407</td><td>87,276,896</td><td>1004</td><td>6</td><td>87,033,142</td><td>rs110121930</td><td>0.613</td><td>1.00</td><td>12.3</td><td>-7.88</td><td>1.09</td><td>-</td><td>Intergenic</td></tr><tr><td>22</td><td>43</td><td>Mg</td><td>6</td><td>88,380,434</td><td>88,436,951</td><td>4</td><td>1</td><td>88,390,406</td><td>rs133413027</td><td>0.252</td><td>0.99</td><td>9.8</td><td>-0.56</td><td>0.09</td><td>-</td><td>Intergenic</td></tr><tr><td>23</td><td>44</td><td>Citrate</td><td>6</td><td>117,331,679</td><td>117,456,476</td><td>30</td><td>7</td><td>117,407,851</td><td>rs43489995</td><td>0.306</td><td>0.30</td><td>10.4</td><td>0.22</td><td>0.03</td><td>TMEM175 / DGKQ</td><td>Downstream /<br>Synonymous</td></tr><tr><td>23</td><td>45</td><td>Ca</td><td>6</td><td><strong>117,331,679</strong></td><td><strong>117,457,454</strong></td><td><strong>42</strong></td><td><strong>6</strong></td><td><strong>117,442,167</strong></td><td>-</td><td><strong>0.407</strong></td><td><strong>0.41</strong></td><td><strong>19.8</strong></td><td><strong>12.7</strong></td><td><strong>1.37</strong></td><td><strong>GAK</strong></td><td><strong>Intron</strong></td></tr><tr><td>24</td><td>46</td><td>Na</td><td>7</td><td>1,055,458</td><td>1,365,056</td><td>387</td><td>3</td><td>1,057,563</td><td>rs210648426</td><td>0.346</td><td>0.91</td><td>14.3</td><td>-2.85</td><td>0.36</td><td>RASGEF1C</td><td>Intron</td></tr><tr><td>25</td><td>47</td><td>Ca</td><td>7</td><td>37,353,946</td><td>39,223,460</td><td>24</td><td>16</td><td>39,194,875</td><td>rs41567828</td><td>0.068</td><td>0.99</td><td>9.1</td><td>12.8</td><td>2.08</td><td>LOC107132625</td><td>IncRNA</td></tr><tr><td>25</td><td>48</td><td>P</td><td>7</td><td><strong>38,917,914</strong></td><td><strong>40,541,631</strong></td><td><strong>10</strong></td><td><strong>5</strong></td><td><strong>40,541,609</strong></td><td><strong>rs1116258782</strong></td><td><strong>0.009</strong></td><td><strong>0.22</strong></td><td><strong>16.9</strong></td><td><strong>31.7</strong></td><td><strong>3.71</strong></td><td>-</td><td><strong>Intergenic</strong></td></tr><tr><td>25</td><td>49</td><td>Ca</td><td>7</td><td>39,528,343</td><td>42,832,365</td><td>179</td><td>16</td><td>41,321,781</td><td>rs209051255</td><td>0.059</td><td>0.76</td><td>11.6</td><td>18.4</td><td>2.63</td><td>ENS-<br>BTAG00000053872</td><td>Upstream BS<br>FOXD3</td></tr><tr><td>26</td><td>50</td><td>Ca</td><td>7</td><td>43,090,012</td><td>44,848,268</td><td>250</td><td>19</td><td>44,836,642</td><td>rs210107035</td><td>0.068</td><td>0.98</td><td>17.1</td><td>21.6</td><td>2.51</td><td>FSTL4 /<br>LOC101902643</td><td>Intron /<br>IncRNA</td></tr></tbody></table>

variants were due to linkage disequilibrium (LD) with the same causal mutation or to the presence of multiple causal mutations. These analyses led to the exclusion of 4 and 9 trait x region combinations on BTA1 and BTA20, respectively, resulting in 83 independent QTL. These QTL corresponded to 50 genomic regions located on *Bos taurus* (BTA) autosomes 1, 2, 3, 4, 5, 6, and 7 (Table 3) and 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 22, 24, 25, 27, 28, and 29 (Table 4). Of these, 18 regions had effects on several traits (2 to 5), while 32 affected only one trait. The
<table><caption>**Table 3.** QTL identified for milk mineral and citrate content on *Bos taurus* (BTA) autosomes 1 to 10 (when a QTL region affected multiple traits, the most significantly affected trait is indicated in bold).</caption><thead><tr><th rowspan="2">Region</th><th rowspan="2">QTL</th><th rowspan="2">Trait</th><th rowspan="2">BTA</th><th colspan="2">Confidence Interval of the QTL</th><th rowspan="2"># variants</th><th rowspan="2"># genes</th><th rowspan="2">Position (bp)</th><th rowspan="2">ID</th><th rowspan="2">p</th><th rowspan="2">R<sup>2</sup></th><th rowspan="2">-log<sub>10</sub>P</th><th rowspan="2">b</th><th rowspan="2">SE</th><th rowspan="2">Name of gene(s)</th><th rowspan="2">Impact</th></tr><tr><th>From (bp)</th><th>To (bp)</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>P</td><td>1</td><td>135,384,680</td><td>136,262,601</td><td>8</td><td>3</td><td>136,262,601</td><td>rs210630730</td><td>0.199</td><td>0.69</td><td>11.5</td><td>8.57</td><td>1.23</td><td>TMEM108</td><td>Intron</td></tr><tr><td>2</td><td>4</td><td>Mg</td><td>1</td><td>142,826,156</td><td>142,837,365</td><td>9</td><td>1</td><td>142,834,737</td><td>rs109459130</td><td>0.576</td><td>0.85</td><td>64.1</td><td>1.48</td><td>0.09</td><td>SLC37A1</td><td>Intron</td></tr><tr><td>2</td><td>5</td><td>P</td><td>1</td><td><strong>142,826,156</strong></td><td><strong>142,835,551</strong></td><td>7</td><td><strong>1</strong></td><td><strong>142,834,737</strong></td><td><strong>rs109459130</strong></td><td><strong>0.576</strong></td><td><strong>0.85</strong></td><td><strong>219.3</strong></td><td><strong>30.2</strong></td><td><strong>0.95</strong></td><td><strong>SLC37A1</strong></td><td><strong>Intron</strong></td></tr><tr><td>2</td><td>6</td><td>K</td><td>1</td><td>142,826,156</td><td>142,838,477</td><td>31</td><td>1</td><td>142,835,551</td><td>rs109717634</td><td>0.613</td><td>1.00</td><td>138.1</td><td>22.4</td><td>0.89</td><td>SLC37A1</td><td>Intron</td></tr><tr><td>2</td><td>7</td><td>Na</td><td>1</td><td>142,826,156</td><td>142,835,551</td><td>7</td><td>1</td><td>142,835,551</td><td>rs109717634</td><td>0.613</td><td>1.00</td><td>19.5</td><td>-3.77</td><td>0.41</td><td>SLC37A1</td><td>Intron</td></tr><tr><td>3</td><td>9</td><td>Mg</td><td>1</td><td>145,114,056</td><td>146,804,023</td><td>64</td><td>3</td><td>146,775,338</td><td>rs43048006</td><td>0.630</td><td>0.33</td><td>10.3</td><td>-0.70</td><td>0.11</td><td>-</td><td>-</td></tr><tr><td>4</td><td>11</td><td>Na</td><td>2</td><td>5,738,340</td><td>7,483,560</td><td>104</td><td>7</td><td>5,739,439</td><td>rs209271012</td><td>0.111</td><td>0.86</td><td>10.4</td><td>3.85</td><td>0.58</td><td>LOC104971101</td><td>IncRNA</td></tr><tr><td>5</td><td>12</td><td>Na</td><td>2</td><td>131,111,763</td><td>131,742,987</td><td>206</td><td>9</td><td>131,197,793</td><td>rs133533873</td><td>0.672</td><td>0.99</td><td>13.8</td><td>2.84</td><td>0.37</td><td>ALPL</td><td>Intron</td></tr><tr><td>6</td><td>13</td><td>Na</td><td>3</td><td>15,459,107</td><td>15,540,579</td><td>7</td><td>6</td><td>15,470,670</td><td>rs110073735</td><td>0.026</td><td>0.92</td><td>8.9</td><td>-5.20</td><td>0.86</td><td>EFNA1 /<br>LOC107132270</td><td>Intron /<br>IncRNA</td></tr><tr><td>7</td><td>14</td><td>Ca</td><td>4</td><td>48,711,128</td><td>48,835,496</td><td>22</td><td>4</td><td>48,711,294</td><td>rs381111559</td><td>0.291</td><td>0.59</td><td>13.0</td><td>-10.5</td><td>1.42</td><td>SLC26A4</td><td>Intron</td></tr><tr><td>8</td><td>15</td><td>Na</td><td>4</td><td>75,153,659</td><td>77,188,370</td><td>14</td><td>7</td><td>77,147,308</td><td>rs386097640</td><td>0.024</td><td>0.87</td><td>9.5</td><td>6.65</td><td>1.06</td><td>GCK</td><td>Intron</td></tr><tr><td>8</td><td>16</td><td>Ca</td><td>4</td><td>77,147,308</td><td>77,257,775</td><td>9</td><td>7</td><td>77,173,366</td><td>rs133824313</td><td>0.025</td><td>0.98</td><td>9.1</td><td>14.8</td><td>2.40</td><td>GCK</td><td>Intron</td></tr><tr><td>8</td><td>17</td><td>Mg</td><td>4</td><td><strong>77,037,103</strong></td><td><strong>77,188,370</strong></td><td><strong>8</strong></td><td><strong>4</strong></td><td><strong>77,173,366</strong></td><td><strong>rs133824313</strong></td><td><strong>0.025</strong></td><td><strong>0.98</strong></td><td><strong>11.1</strong></td><td><strong>1.34</strong></td><td><strong>0.20</strong></td><td><strong>GCK</strong></td><td><strong>Intron</strong></td></tr><tr><td>9</td><td>18</td><td>Ca</td><td>4</td><td>91,807,789</td><td>92,192,288</td><td>182</td><td>6</td><td>92,192,288</td><td>rs42451422</td><td>0.028</td><td>0.95</td><td>9.5</td><td>13.4</td><td>2.13</td><td>LRRC4 / SND1</td><td>Intron /<br>Intron</td></tr><tr><td>10</td><td>19</td><td>Ca</td><td>5</td><td>108,256,181</td><td>109,217,271</td><td>87</td><td>11</td><td>109,141,258</td><td>rs383141568</td><td>0.570</td><td>0.82</td><td>8.4</td><td>5.73</td><td>0.98</td><td>BID</td><td>Intron</td></tr><tr><td>11</td><td>20</td><td>Ca</td><td>5</td><td>111,637,039</td><td>111,848,859</td><td>47</td><td>7</td><td>111,745,513</td><td>rs384814366</td><td>0.231</td><td>0.94</td><td>16.6</td><td>-9.38</td><td>1.11</td><td>MRTFA = MKL1</td><td>Intron</td></tr><tr><td>12</td><td>21</td><td>Ca</td><td>5</td><td><strong>116,403,262</strong></td><td><strong>116,463,104</strong></td><td><strong>48</strong></td><td><strong>8</strong></td><td><strong>116,429,840</strong></td><td><strong>rs382268272</strong></td><td><strong>0.049</strong></td><td><strong>0.58</strong></td><td><strong>60.5</strong></td><td><strong>-34.8</strong></td><td><strong>2.11</strong></td><td><strong>LOC101903383</strong></td><td><strong>IncRNA</strong></td></tr><tr><td>12</td><td>22</td><td>Na</td><td>5</td><td>116,403,262</td><td>116,907,862</td><td>72</td><td>17</td><td>116,429,922</td><td>rs110716509</td><td>0.050</td><td>0.67</td><td>9.3</td><td>5.05</td><td>0.81</td><td>LOC101903383</td><td>IncRNA</td></tr><tr><td>12</td><td>23</td><td>Citrate</td><td>5</td><td>116,403,262</td><td>116,468,630</td><td>50</td><td>8</td><td>116,439,760</td><td>rs110048000</td><td>0.042</td><td>0.56</td><td>9.5</td><td>-0.25</td><td>0.04</td><td>PPARA</td><td>Intron</td></tr><tr><td>12</td><td>24</td><td>Mg</td><td>5</td><td>116,403,806</td><td>116,463,104</td><td>36</td><td>6</td><td>116,448,403</td><td>rs110659119</td><td>0.048</td><td>0.97</td><td>24.2</td><td>-1.38</td><td>0.13</td><td>PPARA</td><td>Intron</td></tr><tr><td>12</td><td>25</td><td>P</td><td>5</td><td>116,420,172</td><td>117,323,874</td><td>28</td><td>5</td><td>116,448,403</td><td>rs110659119</td><td>0.048</td><td>0.97</td><td>19.5</td><td>-13.5</td><td>1.46</td><td>PPARA</td><td>Intron</td></tr><tr><td>13</td><td>26</td><td>Ca</td><td>5</td><td><strong>119,066,644</strong></td><td><strong>119,069,138</strong></td><td><strong>25</strong></td><td><strong>1</strong></td><td><strong>119,068,337</strong></td><td><strong>rs134583661</strong></td><td><strong>0.046</strong></td><td><strong>0.95</strong></td><td><strong>48.4</strong></td><td><strong>-24.8</strong></td><td><strong>1.69</strong></td><td><strong>-</strong></td><td><strong>Intergenic</strong></td></tr><tr><td>13</td><td>27</td><td>Mg</td><td>5</td><td>119,066,644</td><td>119,069,138</td><td>25</td><td>1</td><td>119,068,337</td><td>rs134583661</td><td>0.046</td><td>0.95</td><td>20.1</td><td>-1.29</td><td>0.14</td><td>-</td><td>Intergenic</td></tr><tr><td>13</td><td>28</td><td>P</td><td>5</td><td>119,066,644</td><td>119,069,138</td><td>25</td><td>1</td><td>119,068,337</td><td>rs134583661</td><td>0.046</td><td>0.95</td><td>17.9</td><td>-13.2</td><td>1.50</td><td>-</td><td>Intergenic</td></tr><tr><td>14</td><td>29</td><td>K</td><td>6</td><td>40,012,262</td><td>40,518,049</td><td>18</td><td>2</td><td>40,024,540</td><td>rs43610452</td><td>0.064</td><td>0.99</td><td>11.7</td><td>-11.6</td><td>1.66</td><td>SLIT2</td><td>Intron</td></tr><tr><td>15</td><td>30</td><td>P</td><td>6</td><td>42,855,236</td><td>43,032,274</td><td>17</td><td>3</td><td>42,870,371</td><td>rs42610629</td><td>0.523</td><td>0.97</td><td>9.2</td><td>5.59</td><td>0.91</td><td>-</td><td>Intergenic</td></tr><tr><td>16</td><td>31</td><td>K</td><td>6</td><td>43,297,568</td><td>43,322,676</td><td>12</td><td>1</td><td>43,322,676</td><td>rs43463132</td><td>0.778</td><td>0.97</td><td>40.7</td><td>13.2</td><td>0.98</td><td>-</td><td>Intergenic</td></tr><tr><td>17</td><td>32</td><td>P</td><td>6</td><td>45,018,531</td><td>45,401,570</td><td>78</td><td>4</td><td>45,056,055</td><td>rs378008740</td><td>0.713</td><td>0.76</td><td>30.6</td><td>12.7</td><td>1.09</td><td>-</td><td>Intergenic</td></tr><tr><td>17</td><td>33</td><td>K</td><td>6</td><td><strong>45,018,882</strong></td><td><strong>45,401,570</strong></td><td><strong>66</strong></td><td><strong>4</strong></td><td><strong>45,328,581</strong></td><td><strong>rs135848932</strong></td><td><strong>0.174</strong></td><td><strong>1.00</strong></td><td><strong>120.0</strong></td><td><strong>-27.4</strong></td><td><strong>1.17</strong></td><td><strong>SEL1L3</strong></td><td><strong>Intron</strong></td></tr><tr><td>17</td><td>34</td><td>Ca</td><td>6</td><td>45,328,392</td><td>45,411,116</td><td>17</td><td>4</td><td>45,393,335</td><td>rs134399380</td><td>0.183</td><td>0.91</td><td>8.6</td><td>7.24</td><td>1.22</td><td>-</td><td>Intergenic</td></tr><tr><td>17</td><td>35</td><td>Na</td><td>6</td><td>45,215,035</td><td>45,401,570</td><td>17</td><td>4</td><td>45,401,485</td><td>rs108972810</td><td>0.177</td><td>0.99</td><td>10.2</td><td>-3.39</td><td>0.52</td><td>SMIM20 /<br>LOC112447060</td><td>Upstream<br>BS DOF5.3 /<br>IncRNA</td></tr><tr><td>18</td><td>36</td><td>K</td><td>6</td><td>47,687,737</td><td>49,094,600</td><td>79</td><td>1</td><td>48,145,964</td><td>rs43461552</td><td>0.336</td><td>0.92</td><td>26.9</td><td>-10.6</td><td>0.98</td><td>-</td><td>Intergenic</td></tr><tr><td>19</td><td>37</td><td>K</td><td>6</td><td>50,277,420</td><td>51,315,936</td><td>82</td><td>2</td><td>51,311,430</td><td>rs110769289</td><td>0.555</td><td>0.53</td><td>9.1</td><td>-5.72</td><td>0.93</td><td>-</td><td>Intergenic</td></tr><tr><td>20</td><td>38</td><td>K</td><td>6</td><td>56,960,135</td><td>56,970,541</td><td>6</td><td>1</td><td>56,968,022</td><td>rs134533189</td><td>0.091</td><td>0.41</td><td>9.4</td><td>-8.37</td><td>1.34</td><td>RELL1</td><td>Intron</td></tr><tr><td>21</td><td>39</td><td>Na</td><td>6</td><td>85,409,619</td><td>86,104,900</td><td>312</td><td>15</td><td>85,512,651</td><td>rs110772556</td><td>0.660</td><td>0.98</td><td>8.6</td><td>2.79</td><td>0.47</td><td>-</td><td>Intergenic</td></tr><tr><td>21</td><td>40</td><td>Mg</td><td>6</td><td><strong>85,587,260</strong></td><td><strong>85,599,456</strong></td><td><strong>41</strong></td><td><strong>2</strong></td><td><strong>85,594,649</strong></td><td><strong>rs209134503</strong></td><td><strong>0.481</strong></td><td><strong>0.85</strong></td><td><strong>29.8</strong></td><td><strong>-1.09</strong></td><td><strong>0.10</strong></td><td><strong>ODAM</strong></td><td><strong>Synonymous</strong></td></tr><tr><td>21</td><td>41</td><td>Ca</td><td>6</td><td>85,406,891</td><td>85,599,456</td><td>44</td><td>5</td><td>85,594,717</td><td>rs135113545</td><td>0.485</td><td>0.79</td><td>17.0</td><td>-10.6</td><td>1.24</td><td>ODAM</td><td>Intron</td></tr><tr><td>21</td><td>42</td><td>P</td><td>6</td><td>85,564,407</td><td>87,276,896</td><td>1004</td><td>6</td><td>87,033,142</td><td>rs110121930</td><td>0.613</td><td>1.00</td><td>12.3</td><td>-7.88</td><td>1.09</td><td>-</td><td>Intergenic</td></tr><tr><td>22</td><td>43</td><td>Mg</td><td>6</td><td>88,380,434</td><td>88,436,951</td><td>4</td><td>1</td><td>88,390,406</td><td>rs133413027</td><td>0.252</td><td>0.99</td><td>9.8</td><td>-0.56</td><td>0.09</td><td>-</td><td>Intergenic</td></tr><tr><td>23</td><td>44</td><td>Citrate</td><td>6</td><td>117,331,679</td><td>117,456,476</td><td>30</td><td>7</td><td>117,407,851</td><td>rs43489995</td><td>0.306</td><td>0.30</td><td>10.4</td><td>0.22</td><td>0.03</td><td>TMEM175 / DGKQ</td><td>Downstream /<br>Synonymous</td></tr><tr><td>23</td><td>45</td><td>Ca</td><td>6</td><td><strong>117,331,679</strong></td><td><strong>117,457,454</strong></td><td><strong>42</strong></td><td><strong>6</strong></td><td><strong>117,442,167</strong></td><td><strong>-</strong></td><td><strong>0.407</strong></td><td><strong>0.41</strong></td><td><strong>19.8</strong></td><td><strong>12.7</strong></td><td><strong>1.37</strong></td><td><strong>GAK</strong></td><td><strong>Intron</strong></td></tr><tr><td>24</td><td>46</td><td>Na</td><td>7</td><td>1,055,458</td><td>1,365,056</td><td>387</td><td>3</td><td>1,057,563</td><td>rs210648426</td><td>0.346</td><td>0.91</td><td>14.3</td><td>-2.85</td><td>0.36</td><td>RASGEF1C</td><td>Intron</td></tr><tr><td>25</td><td>47</td><td>Ca</td><td>7</td><td>37,353,946</td><td>39,223,460</td><td>24</td><td>16</td><td>39,194,875</td><td>rs41567828</td><td>0.068</td><td>0.99</td><td>9.1</td><td>12.8</td><td>2.08</td><td>LOC107132625</td><td>IncRNA</td></tr><tr><td>25</td><td>48</td><td>P</td><td>7</td><td><strong>38,917,914</strong></td><td><strong>40,541,631</strong></td><td><strong>10</strong></td><td><strong>5</strong></td><td><strong>40,541,609</strong></td><td><strong>rs1116258782</strong></td><td><strong>0.009</strong></td><td><strong>0.22</strong></td><td><strong>16.9</strong></td><td><strong>31.7</strong></td><td><strong>3.71</strong></td><td><strong>-</strong></td><td><strong>Intergenic</strong></td></tr><tr><td>25</td><td>49</td><td>Ca</td><td>7</td><td>39,528,343</td><td>42,832,365</td><td>179</td><td>16</td><td>41,321,781</td><td>rs209051255</td><td>0.059</td><td>0.76</td><td>11.6</td><td>18.4</td><td>2.63</td><td>ENS-<br>BTAG00000053872</td><td>Upstream BS<br>FOXD3</td></tr><tr><td>26</td><td>50</td><td>Ca</td><td>7</td><td>43,090,012</td><td>44,848,268</td><td>250</td><td>19</td><td>44,836,642</td><td>rs210107035</td><td>0.068</td><td>0.98</td><td>17.1</td><td>21.6</td><td>2.51</td><td>FSTL4 /<br>LOC101902643</td><td>Intron /<br>IncRNA</td></tr></tbody></table>

variants were due to linkage disequilibrium (LD) with the same causal mutation or to the presence of multiple causal mutations. These analyses led to the exclusion of 4 and 9 trait x region combinations on BTA1 and BTA20, respectively, resulting in 83 independent QTL. These QTL corresponded to 50 genomic regions located on *Bos taurus* (BTA) autosomes 1, 2, 3, 4, 5, 6, and 7 (Table 3) and 11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 22, 24, 25, 27, 28, and 29 (Table 4). Of these, 18 regions had effects on several traits (2 to 5), while 32 affected only one trait. The
<table><caption>**Table 4.** QTL identified for milk mineral and citrate content on *Bos taurus* (BTA) autosomes 11 to 29 (when a QTL affected multiple traits, the most significantly affected trait is indicated in bold).</caption><thead><tr><th rowspan="2">Region</th><th rowspan="2">QTL</th><th rowspan="2">Trait</th><th rowspan="2">BTA</th><th colspan="3">Confidence Interval of the QTL</th><th colspan="9">Variant with the most significant effect (frequency p and effect b of the ALT allele)</th><th rowspan="2">Name of gene(s)</th><th rowspan="2">Impact</th></tr><tr><th>From (bp)</th><th>To (bp)</th><th># variants</th><th># genes</th><th>Position (bp)</th><th>ID</th><th>p</th><th>R<sup>2</sup></th><th>-log<sub>10</sub>P</th><th>b</th><th>SE</th></tr></thead><tbody><tr><td>27</td><td>51</td><td>K</td><td>11</td><td>13,798,363</td><td>15,368,263</td><td>563</td><td>7</td><td>14,405,285</td><td>rs110177775</td><td>0.102</td><td>0.96</td><td>10.3</td><td>8.04</td><td>1.22</td><td>-</td><td>Intergenic</td></tr><tr><td>28</td><td>52</td><td>Citrate</td><td>11</td><td>103,239,082</td><td>103,277,147</td><td>279</td><td>5</td><td>103,244,306</td><td>rs435710868</td><td>0.514</td><td>0.85</td><td>23.9</td><td>0.21</td><td>0.02</td><td>ENS-<br/>BTAG00000048091</td><td>Down-<br/>stream</td></tr><tr><td>28</td><td>53</td><td>P</td><td>11</td><td>103,239,082</td><td>103,277,147</td><td>278</td><td>6</td><td>103,250,796</td><td>rs379692402</td><td>0.473</td><td>0.72</td><td>25.0</td><td>-11.8</td><td>1.13</td><td>-</td><td>Intergenic</td></tr><tr><td>29</td><td>54</td><td>Na</td><td>12</td><td>68,052,340</td><td>69,715,887</td><td>87</td><td>7</td><td>68,475,505</td><td>rs133520200</td><td>0.257</td><td>1.00</td><td>9.2</td><td>-2.91</td><td>0.47</td><td>-</td><td>Intergenic</td></tr><tr><td>30</td><td>55</td><td>Na</td><td>12</td><td>73,666,737</td><td>73,725,910</td><td>47</td><td>1</td><td>73,674,445</td><td>rs136814628</td><td>0.680</td><td>0.99</td><td>8.5</td><td>2.53</td><td>0.43</td><td>HS6ST3</td><td>Intron</td></tr><tr><td>31</td><td>56</td><td>K</td><td>13</td><td>54,405,085</td><td>54,637,805</td><td>45</td><td>8</td><td>54,431,192</td><td>rs133004626</td><td>0.623</td><td>0.95</td><td>10.2</td><td>-6.48</td><td>0.99</td><td>-</td><td>Intergenic</td></tr><tr><td>32</td><td>57</td><td>Citrate</td><td>14</td><td>238,942</td><td>613,906</td><td>35</td><td>15</td><td>453,437</td><td>rs384162250</td><td>0.512</td><td>1.00</td><td>11.0</td><td>0.10</td><td>0.01</td><td>PPP1R16A</td><td>Upstream</td></tr><tr><td>32</td><td>58</td><td>K</td><td>14</td><td>232,311</td><td>639,312</td><td>79</td><td>20</td><td>515,185</td><td>rs460225555</td><td>0.501</td><td>0.96</td><td>14.8</td><td>5.60</td><td>0.70</td><td>VPS28</td><td>Down-<br/>stream</td></tr><tr><td>32</td><td>59</td><td>P</td><td>14</td><td>271,581</td><td>758,897</td><td>64</td><td>20</td><td>577,322</td><td>rs133788084</td><td>0.212</td><td>0.21</td><td>13.8</td><td>12.6</td><td>1.64</td><td>FBXL6</td><td>Upstream</td></tr><tr><td>32</td><td>60</td><td>Ca</td><td>14</td><td>284,571</td><td>700,497</td><td>27</td><td>7</td><td>666,982</td><td>rs133196323</td><td>0.025</td><td>0.34</td><td>13.0</td><td>32.2</td><td>4.31</td><td>MROH1</td><td>Intron</td></tr><tr><td>33</td><td>61</td><td>Na</td><td>15</td><td>52,059,631</td><td>52,851,538</td><td>466</td><td>6</td><td>52,461,335</td><td>rs210668186</td><td>0.218</td><td>0.60</td><td>9.2</td><td>-2.65</td><td>0.43</td><td>FCHSD2</td><td>Intron</td></tr><tr><td>34</td><td>62</td><td>K</td><td>17</td><td>50,508,807</td><td>50,988,913</td><td>227</td><td>6</td><td>50,825,796</td><td>rs474259058</td><td>0.055</td><td>0.95</td><td>21.3</td><td>15.1</td><td>1.56</td><td>BRI3BP</td><td>Intron</td></tr><tr><td>35</td><td>63</td><td>Ca</td><td>18</td><td>10,563,963</td><td>10,833,324</td><td>11</td><td>4</td><td>10,833,324</td><td>rs137549452</td><td>0.045</td><td>0.91</td><td>11.3</td><td>-12.9</td><td>1.86</td><td>KLHL36 / COTL1</td><td>Upstream /<br/>Upstream</td></tr><tr><td>36</td><td>64</td><td>Ca</td><td>19</td><td>56,448,441</td><td>56,609,243</td><td>340</td><td>8</td><td>56,524,233</td><td>rs207749796</td><td>0.335</td><td>0.87</td><td>16.1</td><td>7.39</td><td>0.89</td><td>FADS6</td><td>Intron</td></tr><tr><td>37</td><td>65</td><td>Na</td><td>19</td><td>60,550,831</td><td>60,870,548</td><td>95</td><td>2</td><td>60,555,853</td><td>rs435138644</td><td>0.783</td><td>0.61</td><td>23.6</td><td>4.40</td><td>0.43</td><td>-</td><td>Intergenic</td></tr><tr><td>37</td><td>66</td><td>Ca</td><td>19</td><td>60,009,085</td><td>60,870,548</td><td>41</td><td>1</td><td>60,561,566</td><td>rs41923848</td><td>0.843</td><td>0.73</td><td>11.2</td><td>-7.19</td><td>1.05</td><td>-</td><td>Intergenic</td></tr><tr><td>38</td><td>67</td><td>Ca</td><td>20</td><td>31,888,449</td><td>31,888,449</td><td>1</td><td>1</td><td>31,888,449</td><td>rs385640152</td><td>0.009</td><td>0.91</td><td>17.0</td><td>-25.1</td><td>2.92</td><td>GHR /<br/>LOC112443004</td><td>Missense /<br/>lncRNA</td></tr><tr><td>39</td><td>68</td><td>Mg</td><td>20</td><td>31,888,449</td><td>31,888,449</td><td>1</td><td>2</td><td>31,888,449</td><td>rs385640152</td><td>0.009</td><td>0.91</td><td>13.9</td><td>-1.85</td><td>0.24</td><td>GHR /<br/>LOC112443004</td><td>Missense /<br/>lncRNA</td></tr><tr><td>39</td><td>69</td><td>P</td><td>20</td><td>31,888,449</td><td>31,888,449</td><td>1</td><td>1</td><td>31,888,449</td><td>rs385640152</td><td>0.009</td><td>0.91</td><td>10.3</td><td>-17.0</td><td>2.59</td><td>GHR /<br/>LOC112443004</td><td>Missense /<br/>lncRNA</td></tr><tr><td>40</td><td>75</td><td>Ca</td><td>20</td><td>58,185,895</td><td>58,207,145</td><td>21</td><td>1</td><td>58,189,663</td><td>rs110048176</td><td>0.170</td><td>1.00</td><td>25.3</td><td>13.1</td><td>1.24</td><td>ENS-<br/>BTAG00000048498</td><td>lncRNA</td></tr><tr><td>40</td><td>76</td><td>Mg</td><td>20</td><td>58,185,895</td><td>58,388,462</td><td>65</td><td>3</td><td>58,189,663</td><td>rs110048176</td><td>0.170</td><td>1.00</td><td>117.2</td><td>2.34</td><td>0.10</td><td>ENS-<br/>BTAG00000048498</td><td>lncRNA</td></tr><tr><td>40</td><td>77</td><td>Citrate</td><td>20</td><td>58,185,895</td><td>58,388,462</td><td>81</td><td>3</td><td>58,204,929</td><td>rs109956167</td><td>0.166</td><td>0.90</td><td>187.2</td><td>0.73</td><td>0.02</td><td>ENS-<br/>BTAG00000048498</td><td>lncRNA</td></tr><tr><td>40</td><td>78</td><td>K</td><td>20</td><td>58,185,895</td><td>58,388,849</td><td>69</td><td>2</td><td>58,386,888</td><td>rs134021638</td><td>0.054</td><td>0.92</td><td>26.3</td><td>17.4</td><td>1.61</td><td>ANKH</td><td>Intron</td></tr><tr><td>41</td><td>81</td><td>Citrate</td><td>20</td><td>63,223,010</td><td>64,759,735</td><td>87</td><td>2</td><td>63,223,010</td><td>rs109390768</td><td>0.595</td><td>0.98</td><td>14.9</td><td>0.14</td><td>0.02</td><td>-</td><td>Intergenic</td></tr><tr><td>41</td><td>82</td><td>Mg</td><td>20</td><td>63,223,010</td><td>64,753,386</td><td>13</td><td>1</td><td>63,223,010</td><td>rs109390768</td><td>0.595</td><td>0.98</td><td>9.0</td><td>0.47</td><td>0.08</td><td>-</td><td>Intergenic</td></tr><tr><td>42</td><td>85</td><td>Mg</td><td>21</td><td>39,683,083</td><td>41,254,513</td><td>144</td><td>2</td><td>39,933,138</td><td>rs134950504</td><td>0.070</td><td>0.99</td><td>9.3</td><td>-0.85</td><td>0.14</td><td>-</td><td>Intergenic</td></tr><tr><td>43</td><td>86</td><td>P</td><td>22</td><td>28,516,349</td><td>30,130,198</td><td>35</td><td>3</td><td>30,130,198</td><td>rs108953480</td><td>0.966</td><td>0.90</td><td>8.6</td><td>11.4</td><td>1.92</td><td>FOXP1 / ENS-<br/>BTAG00000052330</td><td>Upstream /<br/>lncRNA</td></tr><tr><td>44</td><td>87</td><td>P</td><td>22</td><td>32,774,174</td><td>32,787,867</td><td>10</td><td>2</td><td>32,786,684</td><td>rs208411747</td><td>0.040</td><td>0.92</td><td>20.6</td><td>-14.5</td><td>1.53</td><td>FAM19A4</td><td>Intron</td></tr><tr><td>44</td><td>88</td><td>Mg</td><td>22</td><td>32,772,562</td><td>33,329,180</td><td>45</td><td>4</td><td>32,786,854</td><td>rs209198296</td><td>0.046</td><td>0.82</td><td>13.3</td><td>-1.03</td><td>0.14</td><td>FAM19A4</td><td>Intron</td></tr><tr><td>45</td><td>89</td><td>K</td><td>24</td><td>49,837,381</td><td>50,315,415</td><td>274</td><td>8</td><td>50,074,695</td><td>rs448650804</td><td>0.324</td><td>0.61</td><td>11.1</td><td>5.18</td><td>0.76</td><td>-</td><td>Intergenic</td></tr><tr><td>46</td><td>90</td><td>P</td><td>24</td><td>58,272,345</td><td>58,396,158</td><td>102</td><td>4</td><td>58,306,855</td><td>rs380879212</td><td>0.362</td><td>0.56</td><td>12.8</td><td>5.26</td><td>0.71</td><td>LMAN1</td><td>Intron</td></tr><tr><td>47</td><td>91</td><td>Ca</td><td>25</td><td>20,480,835</td><td>20,482,817</td><td>2</td><td>1</td><td>20,480,835</td><td>rs109069510</td><td>0.888</td><td>0.99</td><td>8.4</td><td>-8.01</td><td>1.36</td><td>HS3ST2</td><td>Intron</td></tr><tr><td>48</td><td>92</td><td>Mg</td><td>27</td><td>36,511,563</td><td>36,592,556</td><td>65</td><td>6</td><td>36,522,002</td><td>rs378026790</td><td>0.488</td><td>0.96</td><td>15.2</td><td>0.63</td><td>0.08</td><td>GPAT4</td><td>Upstream</td></tr><tr><td>49</td><td>93</td><td>K</td><td>28</td><td>6,518,121</td><td>6,534,988</td><td>2</td><td>2</td><td>6,518,121</td><td>rs42033936</td><td>0.224</td><td>0.98</td><td>13.0</td><td>-6.44</td><td>0.86</td><td>KCNK1</td><td>Intron</td></tr><tr><td>49</td><td>94</td><td>Na</td><td>28</td><td>5,987,406</td><td>6,889,024</td><td>58</td><td>3</td><td>6,518,121</td><td>rs42033936</td><td>0.224</td><td>0.98</td><td>10.8</td><td>-2.70</td><td>0.40</td><td>KCNK1</td><td>Intron</td></tr><tr><td>50</td><td>95</td><td>Ca</td><td>29</td><td>9,193,304</td><td>9,907,275</td><td>162</td><td>7</td><td>9,517,882</td><td>rs380735416</td><td>0.163</td><td>0.41</td><td>13.1</td><td>8.67</td><td>1.16</td><td>PICALM</td><td>Upstream</td></tr><tr><td>50</td><td>96</td><td>Na</td><td>29</td><td>8,731,208</td><td>10,144,903</td><td>486</td><td>11</td><td>9,923,696</td><td>rs385315120</td><td>0.166</td><td>0.87</td><td>14.8</td><td>-4.28</td><td>0.54</td><td>TMEM126B</td><td>Down-<br/>stream</td></tr></tbody></table>

four regions with the most significant effects (-log₁₀(P) ≥ 50) each affected four to five traits; these were located on BTA 1 (~142.8 Mbp, for P, K, Mg, and Na; Fig. 2), 20 (~58.2 Mbp, for citrate, Mg, K, and Ca; Fig. 2), 6 (~45.3 Mbp, for K, Ca, P, and Na), and 5 (~116.4 Mbp, for Ca, Mg, P, citrate, and Na). Ten other regions with smaller but still highly significant effects (20 ≤ -log₁₀(P) < 50) were located on BTA 5 (~119.0 Mbp), 6 (~43.3, 48.1, and 85.6 Mbp), 11 (~103.2 Mbp), 17 (~50.8 Mbp), 19 (~60.6 Mbp), 20 (~55.9 and 60.8 Mbp), and 22 (~32.8 Mbp); these also generally affected multiple milk components. All remaining significant QTL regions (8.4 ≤ -log₁₀(P) < 20) were spread across 22 different autosomes. Overall, 6 to 21 different QTL were identified for each trait, and their cumulative effects explained from 19.8% to 45.5% of the genetic variance of a given trait (Table 5). The most significant QTL, located on BTA1 and BTA20, explained more than 23% of the genetic variance of P and citrate, respectively.
<table><caption>**Table 4.** QTL identified for milk mineral and citrate content on *Bos taurus* (BTA) autosomes 11 to 29 (when a QTL affected multiple traits, the most significantly affected trait is indicated in bold).</caption><thead><tr><th rowspan="2">Region</th><th rowspan="2">QTL</th><th rowspan="2">Trait</th><th rowspan="2">BTA</th><th colspan="3">Confidence Interval of the QTL</th><th colspan="10">Variant with the most significant effect (frequency p and effect b of the ALT allele)</th><th rowspan="2">Name of gene(s)</th><th rowspan="2">Impact</th></tr><tr><th>From (bp)</th><th>To (bp)</th><th># variants</th><th># genes</th><th>Position (bp)</th><th>ID</th><th>p</th><th>R<sup>2</sup></th><th>-log<sub>10</sub>P</th><th>b</th><th>SE</th></tr></thead><tbody><tr><td>27</td><td>51</td><td>K</td><td>11</td><td>13,798,363</td><td>15,368,263</td><td>563</td><td>7</td><td>14,405,285</td><td>rs110177775</td><td>0.102</td><td>0.96</td><td>10.3</td><td>8.04</td><td>1.22</td><td>-</td><td>Intergenic</td></tr><tr><td>28</td><td>52</td><td>Citrate</td><td>11</td><td>103,239,082</td><td>103,277,147</td><td>279</td><td>5</td><td>103,244,306</td><td>rs435710868</td><td>0.514</td><td>0.85</td><td>23.9</td><td>0.21</td><td>0.02</td><td>ENS-<br/>BTAG00000048091</td><td>Down-<br/>stream</td></tr><tr><td>28</td><td>53</td><td>P</td><td>11</td><td>103,239,082</td><td>103,277,147</td><td>278</td><td>6</td><td>103,250,796</td><td>rs379692402</td><td>0.473</td><td>0.72</td><td>25.0</td><td>-11.8</td><td>1.13</td><td>-</td><td>Intergenic</td></tr><tr><td>29</td><td>54</td><td>Na</td><td>12</td><td>68,052,340</td><td>69,715,887</td><td>87</td><td>7</td><td>68,475,505</td><td>rs133520200</td><td>0.257</td><td>1.00</td><td>9.2</td><td>-2.91</td><td>0.47</td><td>-</td><td>Intergenic</td></tr><tr><td>30</td><td>55</td><td>Na</td><td>12</td><td>73,666,737</td><td>73,725,910</td><td>47</td><td>1</td><td>73,674,445</td><td>rs136814628</td><td>0.680</td><td>0.99</td><td>8.5</td><td>2.53</td><td>0.43</td><td>HS6ST3</td><td>Intron</td></tr><tr><td>31</td><td>56</td><td>K</td><td>13</td><td>54,405,085</td><td>54,637,805</td><td>45</td><td>8</td><td>54,431,192</td><td>rs133004626</td><td>0.623</td><td>0.95</td><td>10.2</td><td>-6.48</td><td>0.99</td><td>-</td><td>Intergenic</td></tr><tr><td>32</td><td>57</td><td>Citrate</td><td>14</td><td>238,942</td><td>613,906</td><td>35</td><td>15</td><td>453,437</td><td>rs384162250</td><td>0.512</td><td>1.00</td><td>11.0</td><td>0.10</td><td>0.01</td><td>PPP1R16A</td><td>Upstream</td></tr><tr><td>32</td><td>58</td><td>K</td><td>14</td><td>232,311</td><td>639,312</td><td>79</td><td>20</td><td>515,185</td><td>rs460225555</td><td>0.501</td><td>0.96</td><td>14.8</td><td>5.60</td><td>0.70</td><td>VPS28</td><td>Down-<br/>stream</td></tr><tr><td>32</td><td>59</td><td>P</td><td>14</td><td>271,581</td><td>758,897</td><td>64</td><td>20</td><td>577,322</td><td>rs133788084</td><td>0.212</td><td>0.21</td><td>13.8</td><td>12.6</td><td>1.64</td><td>FBXL6</td><td>Upstream</td></tr><tr><td>32</td><td>60</td><td>Ca</td><td>14</td><td>284,571</td><td>700,497</td><td>27</td><td>7</td><td>666,982</td><td>rs133196323</td><td>0.025</td><td>0.34</td><td>13.0</td><td>32.2</td><td>4.31</td><td>MROH1</td><td>Intron</td></tr><tr><td>33</td><td>61</td><td>Na</td><td>15</td><td>52,059,631</td><td>52,851,538</td><td>466</td><td>6</td><td>52,461,335</td><td>rs210668186</td><td>0.218</td><td>0.60</td><td>9.2</td><td>-2.65</td><td>0.43</td><td>FCHSD2</td><td>Intron</td></tr><tr><td>34</td><td>62</td><td>K</td><td>17</td><td>50,508,807</td><td>50,988,913</td><td>227</td><td>6</td><td>50,825,796</td><td>rs474259058</td><td>0.055</td><td>0.95</td><td>21.3</td><td>15.1</td><td>1.56</td><td>BRI3BP</td><td>Intron</td></tr><tr><td>35</td><td>63</td><td>Ca</td><td>18</td><td>10,563,963</td><td>10,833,324</td><td>11</td><td>4</td><td>10,833,324</td><td>rs137549452</td><td>0.045</td><td>0.91</td><td>11.3</td><td>-12.9</td><td>1.86</td><td>KLHL36 / COTL1</td><td>Upstream /<br/>Upstream</td></tr><tr><td>36</td><td>64</td><td>Ca</td><td>19</td><td>56,448,441</td><td>56,609,243</td><td>340</td><td>8</td><td>56,524,233</td><td>rs207749796</td><td>0.335</td><td>0.87</td><td>16.1</td><td>7.39</td><td>0.89</td><td>FADS6</td><td>Intron</td></tr><tr><td>37</td><td>65</td><td>Na</td><td>19</td><td>60,550,831</td><td>60,870,548</td><td>95</td><td>2</td><td>60,555,853</td><td>rs435138644</td><td>0.783</td><td>0.61</td><td>23.6</td><td>4.40</td><td>0.43</td><td>-</td><td>Intergenic</td></tr><tr><td>37</td><td>66</td><td>Ca</td><td>19</td><td>60,009,085</td><td>60,870,548</td><td>41</td><td>1</td><td>60,561,566</td><td>rs41923848</td><td>0.843</td><td>0.73</td><td>11.2</td><td>-7.19</td><td>1.05</td><td>-</td><td>Intergenic</td></tr><tr><td>38</td><td>67</td><td>Ca</td><td>20</td><td>31,888,449</td><td>31,888,449</td><td>1</td><td>1</td><td>31,888,449</td><td>rs385640152</td><td>0.009</td><td>0.91</td><td>17.0</td><td>-25.1</td><td>2.92</td><td>GHR /<br/>LOC112443004</td><td>Missense /<br/>lncRNA</td></tr><tr><td>39</td><td>68</td><td>Mg</td><td>20</td><td>31,888,449</td><td>31,888,449</td><td>1</td><td>2</td><td>31,888,449</td><td>rs385640152</td><td>0.009</td><td>0.91</td><td>13.9</td><td>-1.85</td><td>0.24</td><td>GHR /<br/>LOC112443004</td><td>Missense /<br/>lncRNA</td></tr><tr><td>39</td><td>69</td><td>P</td><td>20</td><td>31,888,449</td><td>31,888,449</td><td>1</td><td>1</td><td>31,888,449</td><td>rs385640152</td><td>0.009</td><td>0.91</td><td>10.3</td><td>-17.0</td><td>2.59</td><td>GHR /<br/>LOC112443004</td><td>Missense /<br/>lncRNA</td></tr><tr><td>40</td><td>75</td><td>Ca</td><td>20</td><td>58,185,895</td><td>58,207,145</td><td>21</td><td>1</td><td>58,189,663</td><td>rs110048176</td><td>0.170</td><td>1.00</td><td>25.3</td><td>13.1</td><td>1.24</td><td>ENS-<br/>BTAG00000048498</td><td>lncRNA</td></tr><tr><td>40</td><td>76</td><td>Mg</td><td>20</td><td>58,185,895</td><td>58,388,462</td><td>65</td><td>3</td><td>58,189,663</td><td>rs110048176</td><td>0.170</td><td>1.00</td><td>117.2</td><td>2.34</td><td>0.10</td><td>ENS-<br/>BTAG00000048498</td><td>lncRNA</td></tr><tr><td>40</td><td>77</td><td>Citrate</td><td>20</td><td>58,185,895</td><td>58,388,462</td><td>81</td><td>3</td><td>58,204,929</td><td>rs109956167</td><td>0.166</td><td>0.90</td><td>187.2</td><td>0.73</td><td>0.02</td><td>ENS-<br/>BTAG00000048498</td><td>lncRNA</td></tr><tr><td>40</td><td>78</td><td>K</td><td>20</td><td>58,185,895</td><td>58,388,849</td><td>69</td><td>2</td><td>58,386,888</td><td>rs134021638</td><td>0.054</td><td>0.92</td><td>26.3</td><td>17.4</td><td>1.61</td><td>ANKH</td><td>Intron</td></tr><tr><td>41</td><td>81</td><td>Citrate</td><td>20</td><td>63,223,010</td><td>64,759,735</td><td>87</td><td>2</td><td>63,223,010</td><td>rs109390768</td><td>0.595</td><td>0.98</td><td>14.9</td><td>0.14</td><td>0.02</td><td>-</td><td>Intergenic</td></tr><tr><td>41</td><td>82</td><td>Mg</td><td>20</td><td>63,223,010</td><td>64,753,386</td><td>13</td><td>1</td><td>63,223,010</td><td>rs109390768</td><td>0.595</td><td>0.98</td><td>9.0</td><td>0.47</td><td>0.08</td><td>-</td><td>Intergenic</td></tr><tr><td>42</td><td>85</td><td>Mg</td><td>21</td><td>39,683,083</td><td>41,254,513</td><td>144</td><td>2</td><td>39,933,138</td><td>rs134950504</td><td>0.070</td><td>0.99</td><td>9.3</td><td>-0.85</td><td>0.14</td><td>-</td><td>Intergenic</td></tr><tr><td>43</td><td>86</td><td>P</td><td>22</td><td>28,516,349</td><td>30,130,198</td><td>35</td><td>3</td><td>30,130,198</td><td>rs108953480</td><td>0.966</td><td>0.90</td><td>8.6</td><td>11.4</td><td>1.92</td><td>FOXP1 / ENS-<br/>BTAG00000052330</td><td>Upstream /<br/>lncRNA</td></tr><tr><td>44</td><td>87</td><td>P</td><td>22</td><td>32,774,174</td><td>32,787,867</td><td>10</td><td>2</td><td>32,786,684</td><td>rs208411747</td><td>0.040</td><td>0.92</td><td>20.6</td><td>-14.5</td><td>1.53</td><td>FAM19A4</td><td>Intron</td></tr><tr><td>44</td><td>88</td><td>Mg</td><td>22</td><td>32,772,562</td><td>33,329,180</td><td>45</td><td>4</td><td>32,786,854</td><td>rs209198296</td><td>0.046</td><td>0.82</td><td>13.3</td><td>-1.03</td><td>0.14</td><td>FAM19A4</td><td>Intron</td></tr><tr><td>45</td><td>89</td><td>K</td><td>24</td><td>49,837,381</td><td>50,315,415</td><td>274</td><td>8</td><td>50,074,695</td><td>rs448650804</td><td>0.324</td><td>0.61</td><td>11.1</td><td>5.18</td><td>0.76</td><td>-</td><td>Intergenic</td></tr><tr><td>46</td><td>90</td><td>P</td><td>24</td><td>58,272,345</td><td>58,396,158</td><td>102</td><td>4</td><td>58,306,855</td><td>rs380879212</td><td>0.362</td><td>0.56</td><td>12.8</td><td>5.26</td><td>0.71</td><td>LMAN1</td><td>Intron</td></tr><tr><td>47</td><td>91</td><td>Ca</td><td>25</td><td>20,480,835</td><td>20,482,817</td><td>2</td><td>1</td><td>20,480,835</td><td>rs109069510</td><td>0.888</td><td>0.99</td><td>8.4</td><td>-8.01</td><td>1.36</td><td>HS3ST2</td><td>Intron</td></tr><tr><td>48</td><td>92</td><td>Mg</td><td>27</td><td>36,511,563</td><td>36,592,556</td><td>65</td><td>6</td><td>36,522,002</td><td>rs378026790</td><td>0.488</td><td>0.96</td><td>15.2</td><td>0.63</td><td>0.08</td><td>GPAT4</td><td>Upstream</td></tr><tr><td>49</td><td>93</td><td>K</td><td>28</td><td>6,518,121</td><td>6,534,988</td><td>2</td><td>2</td><td>6,518,121</td><td>rs42033936</td><td>0.224</td><td>0.98</td><td>13.0</td><td>-6.44</td><td>0.86</td><td>KCNK1</td><td>Intron</td></tr><tr><td>49</td><td>94</td><td>Na</td><td>28</td><td>5,987,406</td><td>6,889,024</td><td>58</td><td>3</td><td>6,518,121</td><td>rs42033936</td><td>0.224</td><td>0.98</td><td>10.8</td><td>-2.70</td><td>0.40</td><td>KCNK1</td><td>Intron</td></tr><tr><td>50</td><td>95</td><td>Ca</td><td>29</td><td>9,193,304</td><td>9,907,275</td><td>162</td><td>7</td><td>9,517,882</td><td>rs380735416</td><td>0.163</td><td>0.41</td><td>13.1</td><td>8.67</td><td>1.16</td><td>PICALM</td><td>Upstream</td></tr><tr><td>50</td><td>96</td><td>Na</td><td>29</td><td>8,731,208</td><td>10,144,903</td><td>486</td><td>11</td><td>9,923,696</td><td>rs385315120</td><td>0.166</td><td>0.87</td><td>14.8</td><td>-4.28</td><td>0.54</td><td>TMEM126B</td><td>Down-<br/>stream</td></tr></tbody></table>

four regions with the most significant effects (-log₁₀(P) ≥ 50) each affected four to five traits; these were located on BTA 1 (~142.8 Mbp, for P, K, Mg, and Na; Fig. 2), 20 (~58.2 Mbp, for citrate, Mg, K, and Ca; Fig. 2), 6 (~45.3 Mbp, for K, Ca, P, and Na), and 5 (~116.4 Mbp, for Ca, Mg, P, citrate, and Na). Ten other regions with smaller but still highly significant effects (20 ≤ -log₁₀(P) < 50) were located on BTA 5 (~119.0 Mbp), 6 (~43.3, 48.1, and 85.6 Mbp), 11 (~103.2 Mbp), 17 (~50.8 Mbp), 19 (~60.6 Mbp), 20 (~55.9 and 60.8 Mbp), and 22 (~32.8 Mbp); these also generally affected multiple milk components. All remaining significant QTL regions (8.4 ≤ -log₁₀(P) < 20) were spread across 22 different autosomes. Overall, 6 to 21 different QTL were identified for each trait, and their cumulative effects explained from 19.8% to 45.5% of the genetic variance of a given trait (Table 5). The most significant QTL, located on BTA1 and BTA20, explained more than 23% of the genetic variance of P and citrate, respectively.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** $-\log_{10}(P)$ values of the effects of variants on milk mineral (Ca, P, Mg, K, and Na) or citrate content plotted against their position on *Bos taurus* autosomes 1 and 20.</figcaption></figure>

<table><caption>**Table 5.** Number of QTL and the percentage of genetic variance explained for each trait.</caption><thead><tr><th rowspan="2">Trait</th><th rowspan="2">Number of QTL</th><th colspan="3">% of genetic variance explained by the QTL</th></tr><tr><th>Total</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td>Ca</td><td>21</td><td>33.1</td><td>0.43</td><td>5.0</td></tr><tr><td>P</td><td>14</td><td>41.7</td><td>0.27</td><td>23.3</td></tr><tr><td>Mg</td><td>13</td><td>31.0</td><td>0.42</td><td>10.6</td></tr><tr><td>K</td><td>14</td><td>42.2</td><td>0.66</td><td>13.6</td></tr><tr><td>Na</td><td>15</td><td>19.8</td><td>0.52</td><td>2.6</td></tr><tr><td>Citrate</td><td>6</td><td>32.9</td><td>0.76</td><td>23.1</td></tr></tbody></table>

**Annotation of GWAS peaks.** The sizes of the confidence intervals (CI) of QTL ranged from 1 bp to 3.3 Mbp; each CI contained between 1 and 1004 variants with significant effects (99.3 on average). While only ca. 6% of the tested variants were included on SNP chips, these chip variants (mainly HD) were disproportionately represented among the variants ranked in the top 10 in QTL peaks, and accounted for nearly 43% of the variants ranked first in the peaks, i.e. “top 1” variants (Fig. 3). For less than one-third of the QTL detected (24/83), the variant with the most significant effect was not located in a gene. Depending on the set of significant variants under consideration (all, top 100, top 50, top 10, or top 1), between 25.1% (top 1 variants) and 28.8% (all vari-
<figure><img src="image_4.png" /><figcaption>**Figure 2.** –log₁₀(*P*) values of the effects of variants on milk mineral (Ca, P, Mg, K, and Na) or citrate content plotted against their position on *Bos taurus* autosomes 1 and 20.</figcaption></figure>

<table><caption>**Table 5.** Number of QTL and the percentage of genetic variance explained for each trait.</caption><thead><tr><th rowspan="2">Trait</th><th rowspan="2">Number of QTL</th><th colspan="3">% of genetic variance explained by the QTL</th></tr><tr><th>Total</th><th>Min</th><th>Max</th></tr></thead><tbody><tr><td>Ca</td><td>21</td><td>33.1</td><td>0.43</td><td>5.0</td></tr><tr><td>P</td><td>14</td><td>41.7</td><td>0.27</td><td>23.3</td></tr><tr><td>Mg</td><td>13</td><td>31.0</td><td>0.42</td><td>10.6</td></tr><tr><td>K</td><td>14</td><td>42.2</td><td>0.66</td><td>13.6</td></tr><tr><td>Na</td><td>15</td><td>19.8</td><td>0.52</td><td>2.6</td></tr><tr><td>Citrate</td><td>6</td><td>32.9</td><td>0.76</td><td>23.1</td></tr></tbody></table>

**Annotation of GWAS peaks.** The sizes of the confidence intervals (CI) of QTL ranged from 1 bp to 3.3 Mbp; each CI contained between 1 and 1004 variants with significant effects (99.3 on average). While only ca. 6% of the tested variants were included on SNP chips, these chip variants (mainly HD) were disproportionately represented among the variants ranked in the top 10 in QTL peaks, and accounted for nearly 43% of the variants ranked first in the peaks, i.e. “top 1” variants (Fig. 3). For less than one-third of the QTL detected (24/83), the variant with the most significant effect was not located in a gene. Depending on the set of significant variants under consideration (all, top 100, top 50, top 10, or top 1), between 25.1% (top 1 variants) and 28.8% (all vari-
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Percentages of chip variants (50 K, HD, or “research” SNPs, in green) and imputed variants (in blue) found among the top 50 variants of QTL peaks. Dashed lines represent the total percentage of chip variants (6%, in green) and imputed variants (94%, in blue) of all variants that were tested in the GWAS.</figcaption></figure>

<table><caption>**Table 6.** Annotation of variants located within confidence intervals of the QTL (a single variant could have different annotations).</caption><thead><tr><th rowspan="2">Annotation</th><th colspan="2">ALL</th><th colspan="2">Top 100</th><th colspan="2">Top 50</th><th colspan="2">Top 10</th><th colspan="2">Top 1</th></tr><tr><th>N</th><th>%</th><th>N</th><th>%</th><th>N</th><th>%</th><th>N</th><th>%</th><th>N</th><th>%</th></tr></thead><tbody><tr><td>Intergenic</td><td>3705</td><td>28.8</td><td>1785</td><td>27.3</td><td>1062</td><td>26.0</td><td>269</td><td>25.1</td><td>29</td><td>25.9</td></tr><tr><td>Total protein coding</td><td>5575</td><td>63.3</td><td>2745</td><td>61.8</td><td>1792</td><td>62.9</td><td>494</td><td>65.7</td><td>49</td><td>60.5</td></tr><tr><td>Intron</td><td>3223</td><td>36.6</td><td>1427</td><td>32.1</td><td>980</td><td>34.4</td><td>267</td><td>35.5</td><td>33</td><td>40.7</td></tr><tr><td>Downstream</td><td>1203</td><td>11.7</td><td>635</td><td>12.1</td><td>368</td><td>11.1</td><td>71</td><td>8.0</td><td>4</td><td>4.2</td></tr><tr><td>Upstream</td><td>948</td><td>9.2</td><td>562</td><td>10.7</td><td>378</td><td>11.4</td><td>136</td><td>15.3</td><td>9</td><td>9.4</td></tr><tr><td>3' UTR</td><td>102</td><td>1.0</td><td>40</td><td>0.76</td><td>23</td><td>0.69</td><td>5</td><td>0.56</td><td>0</td><td>0</td></tr><tr><td>Missense</td><td>84</td><td>0.81</td><td>55</td><td>1.1</td><td>24</td><td>0.78</td><td>4</td><td>0.68</td><td>1</td><td>1.0</td></tr><tr><td>Synonymous</td><td>76</td><td>0.74</td><td>38</td><td>0.73</td><td>28</td><td>0.84</td><td>10</td><td>1.1</td><td>3</td><td>3.1</td></tr><tr><td>5' UTR</td><td>61</td><td>0.59</td><td>44</td><td>0.84</td><td>27</td><td>0.81</td><td>13</td><td>1.5</td><td>0</td><td>0</td></tr><tr><td>Splice region</td><td>11</td><td>0.11</td><td>5</td><td>0.10</td><td>1</td><td>0.03</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Splice acceptor</td><td>1</td><td>0.01</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Frameshift</td><td>1</td><td>0.01</td><td>0</td><td>0.00</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Total non-coding RNA</td><td>693</td><td>7.9</td><td>487</td><td>11.0</td><td>315</td><td>11.1</td><td>69</td><td>9.2</td><td>11</td><td>13.6</td></tr><tr><td>IncRNA</td><td>552</td><td>8.8</td><td>367</td><td>11.4</td><td>237</td><td>11.2</td><td>68</td><td>12.1</td><td>11</td><td>18.3</td></tr><tr><td>miRNA</td><td>69</td><td>1.1</td><td>50</td><td>1.5</td><td>42</td><td>2.0</td><td>1</td><td>0.18</td><td>0</td><td>0.00</td></tr><tr><td>scaRNA</td><td>4</td><td>0.06</td><td>4</td><td>0.12</td><td>4</td><td>0.19</td><td>0</td><td>0.00</td><td>0</td><td>0.00</td></tr><tr><td>snRNA</td><td>67</td><td>1.1</td><td>65</td><td>2.0</td><td>31</td><td>1.5</td><td>0</td><td>0.00</td><td>0</td><td>0.00</td></tr><tr><td>tRNA</td><td>1</td><td>0.02</td><td>1</td><td>0.03</td><td>1</td><td>0.05</td><td>0</td><td>0.00</td><td>0</td><td>0.00</td></tr><tr><td>Total</td><td>10,311</td><td>100</td><td>5234</td><td>100</td><td>3327</td><td>100</td><td>888</td><td>100</td><td>96</td><td>100</td></tr></tbody></table>

ants in the CI of the QTL) of variants were intergenic, while 60.5% (top 1) to 65.7% (top 10) of variants were located in genes, including upstream and downstream regions (Table 6). In all QTL, we therefore identified at least one positional candidate gene (1 to 20 per QTL, 5.3 on average). In total, 271 different genes were found within the CIs of the 83 QTL. In the four QTL regions with the most significant effects, the top-ranked variants were located in the genes *SLC37A1* (BTA1), *ANKH* (BTA20), *SEL1L3* (BTA6), and *PPARA* (BTA5); and the top-ranked variants in the remaining QTL were located in dozens of other genes (Tables 3 and 4). The majority of the variants located in genes were intronic, and were only rarely found in coding regions: only 1%, approximately, of variants were non-synonymous, i.e. 84 variants within the CIs of the QTL, and 55, 24, 4, and 1 among the top 100, top 50, top 10 and top 1 variants, respectively (Table 6). Out of these non-synonymous variants located in coding regions of genes, 4 are serious candidates because they were found among the 10 most significant variants of the QTL peaks: G446S in *SLC26A4* (ranked 2nd for the QTL found on BTA4 at 48.7Mbp for Ca), V298M in *ENSBTAG00000050954* (ranked 10th for the QTL found on BTA7 at 40.3Mbp for P), F267V in *ARHGAP39* (ranked 5th for the QTL found on BTA14 at 0.4Mbp for K and Citrate) and F257Y in *GHR* (ranked 1st for the QTL found on BTA20 at 31.9Mbp for P, Ca and Mg). We also found a relatively high proportion of variants
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Percentages of chip variants (50 K, HD, or “research” SNPs, in green) and imputed variants (in blue) found among the top 50 variants of QTL peaks. Dashed lines represent the total percentage of chip variants (6%, in green) and imputed variants (94%, in blue) of all variants that were tested in the GWAS.</figcaption></figure>

<table><caption>**Table 6.** Annotation of variants located within confidence intervals of the QTL (a single variant could have different annotations).</caption><thead><tr><th rowspan="2">Annotation</th><th colspan="2">ALL</th><th colspan="2">Top 100</th><th colspan="2">Top 50</th><th colspan="2">Top 10</th><th colspan="2">Top 1</th></tr><tr><th>N</th><th>%</th><th>N</th><th>%</th><th>N</th><th>%</th><th>N</th><th>%</th><th>N</th><th>%</th></tr></thead><tbody><tr><td>Intergenic</td><td>3705</td><td>28.8</td><td>1785</td><td>27.3</td><td>1062</td><td>26.0</td><td>269</td><td>25.1</td><td>29</td><td>25.9</td></tr><tr><td>Total protein coding</td><td>5575</td><td>63.3</td><td>2745</td><td>61.8</td><td>1792</td><td>62.9</td><td>494</td><td>65.7</td><td>49</td><td>60.5</td></tr><tr><td>Intron</td><td>3223</td><td>36.6</td><td>1427</td><td>32.1</td><td>980</td><td>34.4</td><td>267</td><td>35.5</td><td>33</td><td>40.7</td></tr><tr><td>Downstream</td><td>1203</td><td>11.7</td><td>635</td><td>12.1</td><td>368</td><td>11.1</td><td>71</td><td>8.0</td><td>4</td><td>4.2</td></tr><tr><td>Upstream</td><td>948</td><td>9.2</td><td>562</td><td>10.7</td><td>378</td><td>11.4</td><td>136</td><td>15.3</td><td>9</td><td>9.4</td></tr><tr><td>3' UTR</td><td>102</td><td>1.0</td><td>40</td><td>0.76</td><td>23</td><td>0.69</td><td>5</td><td>0.56</td><td>0</td><td>0</td></tr><tr><td>Missense</td><td>84</td><td>0.81</td><td>55</td><td>1.1</td><td>24</td><td>0.78</td><td>4</td><td>0.68</td><td>1</td><td>1.0</td></tr><tr><td>Synonymous</td><td>76</td><td>0.74</td><td>38</td><td>0.73</td><td>28</td><td>0.84</td><td>10</td><td>1.1</td><td>3</td><td>3.1</td></tr><tr><td>5' UTR</td><td>61</td><td>0.59</td><td>44</td><td>0.84</td><td>27</td><td>0.81</td><td>13</td><td>1.5</td><td>0</td><td>0</td></tr><tr><td>Splice region</td><td>11</td><td>0.11</td><td>5</td><td>0.10</td><td>1</td><td>0.03</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Splice acceptor</td><td>1</td><td>0.01</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Frameshift</td><td>1</td><td>0.01</td><td>0</td><td>0.00</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td>Total non-coding RNA</td><td>693</td><td>7.9</td><td>487</td><td>11.0</td><td>315</td><td>11.1</td><td>69</td><td>9.2</td><td>11</td><td>13.6</td></tr><tr><td>IncRNA</td><td>552</td><td>8.8</td><td>367</td><td>11.4</td><td>237</td><td>11.2</td><td>68</td><td>12.1</td><td>11</td><td>18.3</td></tr><tr><td>miRNA</td><td>69</td><td>1.1</td><td>50</td><td>1.5</td><td>42</td><td>2.0</td><td>1</td><td>0.18</td><td>0</td><td>0.00</td></tr><tr><td>scaRNA</td><td>4</td><td>0.06</td><td>4</td><td>0.12</td><td>4</td><td>0.19</td><td>0</td><td>0.00</td><td>0</td><td>0.00</td></tr><tr><td>snRNA</td><td>67</td><td>1.1</td><td>65</td><td>2.0</td><td>31</td><td>1.5</td><td>0</td><td>0.00</td><td>0</td><td>0.00</td></tr><tr><td>tRNA</td><td>1</td><td>0.02</td><td>1</td><td>0.03</td><td>1</td><td>0.05</td><td>0</td><td>0.00</td><td>0</td><td>0.00</td></tr><tr><td>Total</td><td>10,311</td><td>100</td><td>5234</td><td>100</td><td>3327</td><td>100</td><td>888</td><td>100</td><td>96</td><td>100</td></tr></tbody></table>

ants in the CI of the QTL) of variants were intergenic, while 60.5% (top 1) to 65.7% (top 10) of variants were located in genes, including upstream and downstream regions (Table 6). In all QTL, we therefore identified at least one positional candidate gene (1 to 20 per QTL, 5.3 on average). In total, 271 different genes were found within the CIs of the 83 QTL. In the four QTL regions with the most significant effects, the top-ranked variants were located in the genes *SLC37A1* (BTA1), *ANKH* (BTA20), *SEL1L3* (BTA6), and *PPARA* (BTA5); and the top-ranked variants in the remaining QTL were located in dozens of other genes (Tables 3 and 4). The majority of the variants located in genes were intronic, and were only rarely found in coding regions: only 1%, approximately, of variants were non-synonymous, i.e. 84 variants within the CIs of the QTL, and 55, 24, 4, and 1 among the top 100, top 50, top 10 and top 1 variants, respectively (Table 6). Out of these non-synonymous variants located in coding regions of genes, 4 are serious candidates because they were found among the 10 most significant variants of the QTL peaks: G446S in *SLC26A4* (ranked 2nd for the QTL found on BTA4 at 48.7Mbp for Ca), V298M in *ENSBTAG00000050954* (ranked 10th for the QTL found on BTA7 at 40.3Mbp for P), F267V in *ARHGAP39* (ranked 5th for the QTL found on BTA14 at 0.4Mbp for K and Citrate) and F257Y in *GHR* (ranked 1st for the QTL found on BTA20 at 31.9Mbp for P, Ca and Mg). We also found a relatively high proportion of variants
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Counts of significantly expressed genes (*P* < 10⁻⁴) in 91 tissues, showing (a) positional candidate genes and (b) transcription factors identified as putative regulators of these genes.</figcaption></figure>

located in non-coding RNA (7.9% of all variants in the CIs of the QTL and 13.6% of top 1 variants). As presented in Table 6, these variants were mostly located in long non-coding RNA (lncRNA), less frequently in micro RNA (miRNA) and small nuclear RNA (snRNA), and very rarely in small nucleolar RNA (small Cajal body-specific RNA, scaRNA), and transfer RNA (tRNA).

We then completed annotation of the QTL (1) identifying putative regulatory SNPs (rSNPs) located within transcription factor binding sites (TFBSs) according to the procedure described in the Methods section, and (2) using transcriptomic data available on the Cattle Gene Atlas website, http://cattlegeneatlas.roslin.ed.ac.uk/.

In the CIs of 40 of the 83 QTL, we identified 184 variants that overlapped with the TFBSs of 438 transcription factors (TFs), i.e. putative rSNPs (see Supplementary Table S1). These 438 TFs had 67 unique target genes among the 271 positional candidates identified in the QTL. Each TF targeted from 1 to 5 positional candidate genes, and 110 TFs had binding sites (BSs) at multiple loci. One TF (ASCL2) targeted 5 genes (*BHLHE23*, *DGKQ*, *RMND5B*, *bta-let-7a-3*, and *bta-mir-2443*), while 2, 21, and 85 TFs targeted 4, 3, and 2 genes, respectively. Two top 1 variants, rs108972810 (~ 45.4 Mbp in the upstream region of *SMIM20* on BTA6, for Na) and rs209051255 (~ 41.3 Mbp in the upstream region of *ENSBTAG00000053872* on BTA7, for Ca) were in the TFBSs for the TFs DOF and FOXD3, respectively. In total, for 22 QTL, top 10 variants were found in the TFBSs of 19 different target genes (including *CSN2*, *PAEP*, *GPAT4*, *BRI3BP*, and *PPARA*), which are potentially regulated by 47 TFs.

We then assessed the degree to which these positional candidate genes and TFs demonstrated tissue-specific expression. Of the 271 genes and 438 TFs associated with milk-mineral QTL, 203 and 160, respectively, were present in the Cattle Gene Atlas, which contains expression data from 91 different tissues or cell lines. Using the model described in the Methods section, we assessed the overexpression or tissue specificity of these genes by evaluating the *t*-statistic and the *P* value, *Pt*, associated with the effect of the tissue under consideration (see Supplementary Tables S2 and S3). Among the 203 positional candidate genes and 160 TFs with a putative regulatory effect on these genes, we calculated the number that were overexpressed in each tissue type (*Pt* < 10⁻⁴). The profiles obtained for both categories of genes, presented in Fig. 4, were quite similar. The tissues or cell types in which we found the highest number of genes associated with milk mineral composition were white blood cells (51 candidate genes and 59 TFs) and mammary gland (47 candidate genes and 45 TFs). In each tissue, we identified the most-specific genes by retaining the top 10% with respect to *t*-statistic values, i.e. 20 candidate genes and 16 TFs. Among the candidate genes that were most specific to mammary gland, we found the genes encoding the main milk proteins (*CSN1S1*, *CSN1S2*, *CSN2*, *CSN3*, and *PAEP*), and *SLC37A1* and *ANKH*, which are located within the QTL that had the most significant effects in the present study. The ASCL2 TF, which potentially regulates five genes located in the CIs of the QTL, was one of the most mammary-gland-specific TFs. Among genes specific to white blood cells, we found *COTL1*, *MKL1*, and *FOXP1*. *SLC37A1* was not among the 20 most-specific genes but it was ranked 22nd. In the upstream or intronic regions of these four genes, we identified the top-ranked variant of four different QTL, located on BTA18, BTA5, BTA22, and BTA1, respectively. Furthermore,
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Counts of significantly expressed genes (*P* < 10⁻⁴) in 91 tissues, showing (a) positional candidate genes and (b) transcription factors identified as putative regulators of these genes.</figcaption></figure>

located in non-coding RNA (7.9% of all variants in the CIs of the QTL and 13.6% of top 1 variants). As presented in Table 6, these variants were mostly located in long non-coding RNA (lncRNA), less frequently in micro RNA (miRNA) and small nuclear RNA (snRNA), and very rarely in small nucleolar RNA (small Cajal body-specific RNA, scaRNA), and transfer RNA (tRNA).

We then completed annotation of the QTL (1) identifying putative regulatory SNPs (rSNPs) located within transcription factor binding sites (TFBSs) according to the procedure described in the Methods section, and (2) using transcriptomic data available on the Cattle Gene Atlas website, http://cattlegeneatlas.roslin.ed.ac.uk/.

In the CIs of 40 of the 83 QTL, we identified 184 variants that overlapped with the TFBSs of 438 transcription factors (TFs), i.e. putative rSNPs (see Supplementary Table S1). These 438 TFs had 67 unique target genes among the 271 positional candidates identified in the QTL. Each TF targeted from 1 to 5 positional candidate genes, and 110 TFs had binding sites (BSs) at multiple loci. One TF (ASCL2) targeted 5 genes (*BHLHE23*, *DGKQ*, *RMND5B*, *bta-let-7a-3*, and *bta-mir-2443*), while 2, 21, and 85 TFs targeted 4, 3, and 2 genes, respectively. Two top 1 variants, rs108972810 (~ 45.4 Mbp in the upstream region of *SMIM20* on BTA6, for Na) and rs209051255 (~ 41.3 Mbp in the upstream region of *ENSBTAG00000053872* on BTA7, for Ca) were in the TFBSs for the TFs DOF and FOXD3, respectively. In total, for 22 QTL, top 10 variants were found in the TFBSs of 19 different target genes (including *CSN2*, *PAEP*, *GPAT4*, *BRI3BP*, and *PPARA*), which are potentially regulated by 47 TFs.

We then assessed the degree to which these positional candidate genes and TFs demonstrated tissue-specific expression. Of the 271 genes and 438 TFs associated with milk-mineral QTL, 203 and 160, respectively, were present in the Cattle Gene Atlas, which contains expression data from 91 different tissues or cell lines. Using the model described in the Methods section, we assessed the overexpression or tissue specificity of these genes by evaluating the *t*-statistic and the *P* value, *Pt*, associated with the effect of the tissue under consideration (see Supplementary Tables S2 and S3). Among the 203 positional candidate genes and 160 TFs with a putative regulatory effect on these genes, we calculated the number that were overexpressed in each tissue type (*Pt* < 10⁻⁴). The profiles obtained for both categories of genes, presented in Fig. 4, were quite similar. The tissues or cell types in which we found the highest number of genes associated with milk mineral composition were white blood cells (51 candidate genes and 59 TFs) and mammary gland (47 candidate genes and 45 TFs). In each tissue, we identified the most-specific genes by retaining the top 10% with respect to *t*-statistic values, i.e. 20 candidate genes and 16 TFs. Among the candidate genes that were most specific to mammary gland, we found the genes encoding the main milk proteins (*CSN1S1*, *CSN1S2*, *CSN2*, *CSN3*, and *PAEP*), and *SLC37A1* and *ANKH*, which are located within the QTL that had the most significant effects in the present study. The ASCL2 TF, which potentially regulates five genes located in the CIs of the QTL, was one of the most mammary-gland-specific TFs. Among genes specific to white blood cells, we found *COTL1*, *MKL1*, and *FOXP1*. *SLC37A1* was not among the 20 most-specific genes but it was ranked 22nd. In the upstream or intronic regions of these four genes, we identified the top-ranked variant of four different QTL, located on BTA18, BTA5, BTA22, and BTA1, respectively. Furthermore,
<figure><img src="image_8.png" /><figcaption>**Figure 5.** SLC37A1 and ANKH genotypes and genotype interaction effects on milk mineral (Ca, P, Mg, K, and Na) or citrate content.</figcaption></figure>

we identified 27 rSNPs for which both TF and target gene were highly specific to the mammary gland; these combinations were composed of 26 unique TFs and 15 unique target genes (see Supplementary Table S4).

**Focus on the SLC37A1 and ANKH gene regions.** As mentioned above, two QTL—located on BTA1 (142.8 Mbp) and BTA20 (58.2 Mbp) in the vicinity of the *SLC37A1* and *ANKH* genes, respectively—had very strong effects on milk mineral and citrate content. The *SLC37A1* region affected P, K, Mg, and Na levels while the *ANKH* region had effects on citrate, Mg, K, and Ca (Fig. 2). These two regions alone affected all six milk components analyzed in this study. They explained a high proportion of the genetic variance in P (23.3%), citrate (23.1%), Mg (17.9%, i.e. 7.3% for *SLC37A1* and 10.6% for *ANKH*), and K (15.4%, i.e. 13.6% for *SLC37A1* and 1.8% for *ANKH*), but had much more moderate effects on Na (2.6%) and Ca (2.2%).

In the *SLC37A1* region, the variant with the most significant effect was located in an intron; it was an imputed variant ($R^2=0.85$) for P and Mg (rs109459130 at 142,834,737 bp) and a chip variant (HD SNP) for K and Na (rs109717634 at 142,835,551 bp). For both variants, which were located 814 bp apart, the most frequent allele (0.58 and 0.61, respectively) increased the amount of all minerals affected by this region.

In the *ANKH* region, we identified three different top 1 variants, depending on the trait. For K, the top variant was rs134021638 (at 58,386,888 bp, imputed with $R^2=0.92$), located in an intronic region of *ANKH*; for Ca, Mg and citrate, the top variants were located in an *lncRNA ENSBTAG00000048498* (rs110048176 at 58,189,663 bp, HD SNP for Ca and Mg and rs109956167 at 58,204,929 bp, imputed variant with $R^2=0.90$ for citrate). *ENSBTAG00000048498* spans from 58,186,301 to 58,224,434 bp and is located around 83 kbp upstream of the *ANKH* gene (58,307,527–58,477,497 bp). In all cases, the alleles responsible for an increase in citrate and mineral content were not the most common, with a MAF of 0.05 for the *ANKH* variant (rs134021638) and 0.17 for the two other variants.

In both regions, we examined the top-ranked variant for the most-affected traits, i.e. rs109459130 and rs109956167, and tested the interaction effects of their genotypes on milk mineral and citrate content. The results, presented in Fig. 5, revealed no significant interaction for any of the traits analyzed ($P \ge 0.02$), suggesting that the effects of the two regions on milk composition were additive.

**Discussion**

To the best of our knowledge, this study is unique because it is the first GWAS of imputed whole-genome sequences based on the most-recent bovine reference genome with such a large population of cattle (19,586 cows); and it is the first attempt to investigate milk mineral content with a sequence-based GWAS, which here assessed a very large panel of genomic polymorphisms (12.9 million).

Using this approach, we identified 83 QTL that explained a substantial part of the genetic variance (up to 42.2%) for mineral and citrate content in cows' milk; in each QTL region, we then identified functional candidate genes and variants. Our results build on those of two previous studies that used GWAS to investigate the mineral composition of milk—based on mineral content measured with reference methods and HD 777 k SNP genotypes—in small populations of Holstein (371⁹ and 444¹²) or Jersey (321⁹) cows. The study of Buitenhuis et al.⁹ was dedicated to mineral composition while Kemper et al.¹² attempted to identify QTL that overlapped between milk production and composition traits. Only two QTL were shared between the two studies: a region located on BTA1 that affected P, for which both studies identified *SLC37A1* as the best functional candidate gene, and a region located near the *DGAT1* gene on BTA14 that affected Ca and P. In our study, we also identified these two regions, together with dozens of other regions located throughout the bovine genome.

By considering a very large dataset (more than 1 million test-day records of 126,873 cows), we estimated heritability values for mineral and citrate composition to be moderate to high, i.e. higher⁵,⁷,⁸,¹⁰ than or similar⁶,⁹ to those previously reported in the literature. Furthermore, we found strong genetic correlations among Ca, P, and Mg and among Ca, Mg, and citrate, which was consistent with the studies of Toffanin et al.⁵ and Denholm et al.⁸. Among minerals, Na had both the lowest heritability value (0.32 vs 0.48–0.56 for the other minerals)
<figure><img src="image_7.png" /><figcaption>**Figure 5.** *SLC37A1* and *ANKH* genotypes and genotype interaction effects on milk mineral (Ca, P, Mg, K, and Na) or citrate content.</figcaption></figure>

we identified 27 rSNPs for which both TF and target gene were highly specific to the mammary gland; these combinations were composed of 26 unique TFs and 15 unique target genes (see Supplementary Table S4).

**Focus on the SLC37A1 and ANKH gene regions.** As mentioned above, two QTL—located on BTA1 (142.8 Mbp) and BTA20 (58.2 Mbp) in the vicinity of the *SLC37A1* and *ANKH* genes, respectively—had very strong effects on milk mineral and citrate content. The *SLC37A1* region affected P, K, Mg, and Na levels while the *ANKH* region had effects on citrate, Mg, K, and Ca (Fig. 2). These two regions alone affected all six milk components analyzed in this study. They explained a high proportion of the genetic variance in P (23.3%), citrate (23.1%), Mg (17.9%, i.e. 7.3% for *SLC37A1* and 10.6% for *ANKH*), and K (15.4%, i.e. 13.6% for *SLC37A1* and 1.8% for *ANKH*), but had much more moderate effects on Na (2.6%) and Ca (2.2%).

In the *SLC37A1* region, the variant with the most significant effect was located in an intron; it was an imputed variant ($R^2=0.85$) for P and Mg (rs109459130 at 142,834,737 bp) and a chip variant (HD SNP) for K and Na (rs109717634 at 142,835,551 bp). For both variants, which were located 814 bp apart, the most frequent allele (0.58 and 0.61, respectively) increased the amount of all minerals affected by this region.

In the *ANKH* region, we identified three different top 1 variants, depending on the trait. For K, the top variant was rs134021638 (at 58,386,888 bp, imputed with $R^2=0.92$), located in an intronic region of *ANKH*; for Ca, Mg and citrate, the top variants were located in an *lncRNA ENSBTAG00000048498* (rs110048176 at 58,189,663 bp, HD SNP for Ca and Mg and rs109956167 at 58,204,929 bp, imputed variant with $R^2=0.90$ for citrate). *ENSBTAG00000048498* spans from 58,186,301 to 58,224,434 bp and is located around 83 kbp upstream of the *ANKH* gene (58,307,527–58,477,497 bp). In all cases, the alleles responsible for an increase in citrate and mineral content were not the most common, with a MAF of 0.05 for the *ANKH* variant (rs134021638) and 0.17 for the two other variants.

In both regions, we examined the top-ranked variant for the most-affected traits, i.e. rs109459130 and rs109956167, and tested the interaction effects of their genotypes on milk mineral and citrate content. The results, presented in Fig. 5, revealed no significant interaction for any of the traits analyzed ($P \ge 0.02$), suggesting that the effects of the two regions on milk composition were additive.

Discussion

To the best of our knowledge, this study is unique because it is the first GWAS of imputed whole-genome sequences based on the most-recent bovine reference genome with such a large population of cattle (19,586 cows); and it is the first attempt to investigate milk mineral content with a sequence-based GWAS, which here assessed a very large panel of genomic polymorphisms (12.9 million).

Using this approach, we identified 83 QTL that explained a substantial part of the genetic variance (up to 42.2%) for mineral and citrate content in cows' milk; in each QTL region, we then identified functional candidate genes and variants. Our results build on those of two previous studies that used GWAS to investigate the mineral composition of milk—based on mineral content measured with reference methods and HD 777 k SNP genotypes—in small populations of Holstein (371⁹ and 444¹²) or Jersey (321⁹) cows. The study of Buitenhuis et al.⁹ was dedicated to mineral composition while Kemper et al.¹² attempted to identify QTL that overlapped between milk production and composition traits. Only two QTL were shared between the two studies: a region located on BTA1 that affected P, for which both studies identified *SLC37A1* as the best functional candidate gene, and a region located near the *DGAT1* gene on BTA14 that affected Ca and P. In our study, we also identified these two regions, together with dozens of other regions located throughout the bovine genome.

By considering a very large dataset (more than 1 million test-day records of 126,873 cows), we estimated heritability values for mineral and citrate composition to be moderate to high, i.e. higher⁵,⁷,⁸,¹⁰ than or similar⁶,⁹ to those previously reported in the literature. Furthermore, we found strong genetic correlations among Ca, P, and Mg and among Ca, Mg, and citrate, which was consistent with the studies of Toffanin et al.⁵ and Denholm et al.⁸. Among minerals, Na had both the lowest heritability value (0.32 vs 0.48–0.56 for the other minerals)
and the smallest percentage of genetic variance explained by its associated QTL (19.8% vs 33.1–42.2% for the other minerals). Ca, P, and Mg shared more QTL than any other group of minerals: 7, 6, and 6 QTL were shared between Ca and P, Ca and Mg, and P and Mg, respectively, while 4 QTL had effects on all three. In contrast, despite the high degree of genetic correlation between Ca, Mg, and citrate, only 2 QTL were shared among these three traits; however, one of these was a QTL located on BTA20 (~58.2 Mbp) that explained 2.2%, 10.6%, and 23.1% of the genetic variance of Ca, Mg, and citrate, respectively. Overall, estimates of genetic correlations were consistent with QTL results, which probably reflects the common biological pathway of these minerals in milk. In colloidal form, Ca, P, and Mg are associated with caseins in the micelles while in the soluble fraction, Ca and Mg are associated with citrate. Further studies, using random regression models fitted across lactation, could investigate the pattern of genetic relationships between the different milk minerals during lactation and thus provide a better understanding of the underlying biological mechanisms.

The resolution of our study was high enough to identify a single or a few positional candidate genes in most of the 83 QTL we identified. In each of these regions, though, the variant with the most significant effect was not necessarily located in a gene. As an example, we detected the overrepresentation of chip SNPs at the top of the QTL peaks (42.7% of the top 1 variants), and the majority of these SNPs were located in intergenic regions. Chip SNPs are directly genotyped or have a higher imputation accuracy than the surrounding variants, which probably enabled more precise estimation of their effects. To identify the best candidate variant(s) for each QTL, we therefore considered not only the variant with the most significant effect, but the set of variants ranked at the top of the peak that were located in or close to genes. In most cases, we found that these variants were located in non-coding regions of genes; for example, 15.3% and 12.1% of all top 10 variants were located in upstream regions of genes or in lncRNA, respectively, i.e. in putative regulatory regions. To further refine our results, we consulted in silico annotations of rSNPs in upstream regions of genes as well as existing knowledge regarding genomic regions that are transcribed into non-coding RNA. Moreover, to support the putative roles of these variants, we also considered the tissue specificity of candidate genes and transcription factors. Unfortunately, the expression dataset, which was based on Ensembl release 94, did not contain all of the genes located in QTL regions. Here we present some examples of QTL in which the putative causal mutation was located in a regulatory region, with particular attention to functional candidate genes previously associated with milk composition and the QTL with the most significant effects in this study.

GPAT4 (*glycerol-3-phosphate acyltransferase 4*) was the best candidate gene for the QTL on BTA 27 with effects on Mg. Although it was not the most notable QTL in terms of its effect on minerals, this region, and this gene in particular, have been highlighted by previous studies as affecting milk composition (fat, protein, and lactose content)¹⁹,²²⁻²⁴. In these studies, five variants in linkage disequilibrium—four in the upstream region and one in the 5' UTR region of GPAT4—were highlighted as the best candidate causal variants. In our study, these variants were ranked in the top 5 in the GWAS peak; the variants ranked 3rd (rs209479876) and 5th (rs209855549) alter the WRKY48 and TWI transcription factor binding sites (TFBSs), respectively, and are more likely to be the causal variants. Daetwyler et al.²³ also highlighted rs209479876 as the best causal variant in this region because of the high probability that it overlaps a TFBS. We also confirmed that GPAT4 was overexpressed in the mammary glands (P = 3.10⁻¹²) compared to 90 other tissues or cell types, but expression data were not available for the TFs associated with this gene.

For the QTL identified on BTA11 at ~103.2 Mbp with effects on P and citrate, the top two variants in the peak for P (the most-affected trait) were intergenic. Instead, the variants ranked 3rd to 7th were located in the upstream region of the *PAEP* (*progestogen-associated endometrial protein*) gene; the 5th-ranked variant, rs110710904, probably alters a binding site of the Macho-1 transcription factor (P = 0.02). *PAEP*, which is one of the most mammary-gland-specific genes (Pt = 1.8.10⁻³⁷), encodes β-lactoglobulin, the most abundant whey protein in cow milk. Two non-synonymous variants in this gene, rs109625649 and rs110066229, were previously highlighted as the causal mutations underlying β-lactoglobulin concentration in milk²⁵. In our study, these were respectively ranked 157th and 171st in the peak, i.e. far below the top-ranked variants upstream. These results corroborate previous reports that these two putatively causative missense mutations did not explain all the effects of the region on milk composition¹⁹,²⁰ and highlight an rSNP as a likely causative variant.

BRI3BP (*BRI3 binding protein*), located on BTA 17 at ~50.8 Mbp, has been previously associated with de novo short chain fatty acid synthesis in bovine milk²⁶,²⁷. In our study, this region appeared to have effects on the K content of milk, and we found that *BRI3BP* was significantly overexpressed in mammary gland compared to the 90 other tissues investigated (P = 3.9.10⁻⁸). The first seven variants in the peak were located in an intronic region of *BRI3BP*. The 8th (rs477456528) and 9th (rs440703666) variants, instead, were located in the upstream region of the gene and rs440703666 probably alter the binding site of the transcription factor SPL8. This variant thus represents a very attractive functional candidate in this region.

In the four regions with the most significant effects on milk minerals, we identified *SLC37A1* (BTA 1), *ANKH* (BTA 20), *SEL1L3* / *SMIM20* (BTA 6), and *PPARA* (BTA 5) as the best candidate genes.

Five of the six traits analyzed in this study were affected by a QTL region at ~116.4 Mbp on BTA5. Here, a variant located at 116,438,773 bp (1000G_80994138) ranked 2nd for four traits. This variant was located in the upstream region of *PPARA* (*peroxisome proliferator activated receptor alpha*) and overlapped the binding sites for the transcription factors E2F4, RSC3, RSC30, TDA9, and TFDP1. In this QTL region, top 1 variants were located either in an intronic region of *PPARA* (citrate, Mg, and P) or in an lncRNA, *LOC101903383* (Ca and Na). *PPARA* encodes the PPAR-α transcription factor, a key regulator of lipid metabolism belonging to the superfamily of PPAR hormone receptors²⁸; expression of this TF in mammary gland was previously found to be associated with milk fatty-acid composition in dairy cows²⁹,³⁰. Here, we did not detect any mammary-gland specificity of *PPARA* but we did observe that this gene was expressed in this tissue. It is therefore possible that polymorphisms in the binding sites of this TF or in the lncRNA *LOC101903383*, located 1.3 kbp upstream of *PPARA*, could regulate the expression of *PPARA* and be responsible for the strong effects of this region on milk composition.
and the smallest percentage of genetic variance explained by its associated QTL (19.8% vs 33.1–42.2% for the other minerals). Ca, P, and Mg shared more QTL than any other group of minerals: 7, 6, and 6 QTL were shared between Ca and P, Ca and Mg, and P and Mg, respectively, while 4 QTL had effects on all three. In contrast, despite the high degree of genetic correlation between Ca, Mg, and citrate, only 2 QTL were shared among these three traits; however, one of these was a QTL located on BTA20 (~58.2 Mbp) that explained 2.2%, 10.6%, and 23.1% of the genetic variance of Ca, Mg, and citrate, respectively. Overall, estimates of genetic correlations were consistent with QTL results, which probably reflects the common biological pathway of these minerals in milk. In colloidal form, Ca, P, and Mg are associated with caseins in the micelles while in the soluble fraction, Ca and Mg are associated with citrate. Further studies, using random regression models fitted across lactation, could investigate the pattern of genetic relationships between the different milk minerals during lactation and thus provide a better understanding of the underlying biological mechanisms.

The resolution of our study was high enough to identify a single or a few positional candidate genes in most of the 83 QTL we identified. In each of these regions, though, the variant with the most significant effect was not necessarily located in a gene. As an example, we detected the overrepresentation of chip SNPs at the top of the QTL peaks (42.7% of the top 1 variants), and the majority of these SNPs were located in intergenic regions. Chip SNPs are directly genotyped or have a higher imputation accuracy than the surrounding variants, which probably enabled more precise estimation of their effects. To identify the best candidate variant(s) for each QTL, we therefore considered not only the variant with the most significant effect, but the set of variants ranked at the top of the peak that were located in or close to genes. In most cases, we found that these variants were located in non-coding regions of genes; for example, 15.3% and 12.1% of all top 10 variants were located in upstream regions of genes or in lncRNA, respectively, i.e. in putative regulatory regions. To further refine our results, we consulted in silico annotations of rSNPs in upstream regions of genes as well as existing knowledge regarding genomic regions that are transcribed into non-coding RNA. Moreover, to support the putative roles of these variants, we also considered the tissue specificity of candidate genes and transcription factors. Unfortunately, the expression dataset, which was based on Ensembl release 94, did not contain all of the genes located in QTL regions. Here we present some examples of QTL in which the putative causal mutation was located in a regulatory region, with particular attention to functional candidate genes previously associated with milk composition and the QTL with the most significant effects in this study.

GPAT4 (glycerol-3-phosphate acyltransferase 4) was the best candidate gene for the QTL on BTA 27 with effects on Mg. Although it was not the most notable QTL in terms of its effect on minerals, this region, and this gene in particular, have been highlighted by previous studies as affecting milk composition (fat, protein, and lactose content)¹⁹,²²⁻²⁴. In these studies, five variants in linkage disequilibrium—four in the upstream region and one in the 5' UTR region of GPAT4—were highlighted as the best candidate causal variants. In our study, these variants were ranked in the top 5 in the GWAS peak; the variants ranked 3rd (rs209479876) and 5th (rs209855549) alter the WRKY48 and TWI transcription factor binding sites (TFBSs), respectively, and are more likely to be the causal variants. Daetwyler et al.²³ also highlighted rs209479876 as the best causal variant in this region because of the high probability that it overlaps a TFBS. We also confirmed that GPAT4 was overexpressed in the mammary glands (P = 3.10⁻¹²) compared to 90 other tissues or cell types, but expression data were not available for the TFs associated with this gene.

For the QTL identified on BTA11 at ~103.2 Mbp with effects on P and citrate, the top two variants in the peak for P (the most-affected trait) were intergenic. Instead, the variants ranked 3rd to 7th were located in the upstream region of the *PAEP* (*progestogen-associated endometrial protein*) gene; the 5th-ranked variant, rs110710904, probably alters a binding site of the Macho-1 transcription factor (P = 0.02). *PAEP*, which is one of the most mammary-gland-specific genes (*Pt* = 1.8.10⁻³⁷), encodes β-lactoglobulin, the most abundant whey protein in cow milk. Two non-synonymous variants in this gene, rs109625649 and rs110066229, were previously highlighted as the causal mutations underlying β-lactoglobulin concentration in milk²⁵. In our study, these were respectively ranked 157th and 171st in the peak, i.e. far below the top-ranked variants upstream. These results corroborate previous reports that these two putatively causative missense mutations did not explain all the effects of the region on milk composition¹⁹,²⁰ and highlight an rSNP as a likely causative variant.

BRI3BP (*BRI3 binding protein*), located on BTA 17 at ~50.8 Mbp, has been previously associated with de novo short chain fatty acid synthesis in bovine milk²⁶,²⁷. In our study, this region appeared to have effects on the K content of milk, and we found that *BRI3BP* was significantly overexpressed in mammary gland compared to the 90 other tissues investigated (P = 3.9.10⁻⁸). The first seven variants in the peak were located in an intronic region of *BRI3BP*. The 8th (rs477456528) and 9th (rs440703666) variants, instead, were located in the upstream region of the gene and rs440703666 probably alter the binding site of the transcription factor SPL8. This variant thus represents a very attractive functional candidate in this region.

In the four regions with the most significant effects on milk minerals, we identified *SLC37A1* (BTA 1), *ANKH* (BTA 20), *SEL1L3* / *SMIM20* (BTA 6), and *PPARA* (BTA 5) as the best candidate genes.

Five of the six traits analyzed in this study were affected by a QTL region at ~116.4 Mbp on BTA5. Here, a variant located at 116,438,773 bp (1000G_80994138) ranked 2nd for four traits. This variant was located in the upstream region of *PPARA* (*peroxisome proliferator activated receptor alpha*) and overlapped the binding sites for the transcription factors E2F4, RSC3, RSC30, TDA9, and TFDP1. In this QTL region, top 1 variants were located either in an intronic region of *PPARA* (citrate, Mg, and P) or in an lncRNA, *LOC101903383* (Ca and Na). *PPARA* encodes the PPAR-α transcription factor, a key regulator of lipid metabolism belonging to the superfamily of PPAR hormone receptors²⁸; expression of this TF in mammary gland was previously found to be associated with milk fatty-acid composition in dairy cows²⁹,³⁰. Here, we did not detect any mammary-gland specificity of *PPARA* but we did observe that this gene was expressed in this tissue. It is therefore possible that polymorphisms in the binding sites of this TF or in the lncRNA *LOC101903383*, located 1.3 kbp upstream of *PPARA*, could regulate the expression of *PPARA* and be responsible for the strong effects of this region on milk composition.
On BTA6 at ~45.3 Mbp, we identified a QTL with very strong significant effects on K and to a lesser extent on P, Na, and Ca. Within the confidence interval of this QTL, we found 80 different potential variants, of which the majority (69) were intergenic, 9 were located in an intronic region of *SEL1L3* (*SEL1L family member 3*), and 2 were in the upstream region of *SMIM20* (*small integral membrane protein 20*). However, when we examined only the top 10 variants for the different traits, the number of distinct variants shrank to 18 (9 intergenic, 7 in *SEL1L3*, and 2 in *SMIM20*). Of these, two (rs108972810 at 45,401,485 bp and rs136498639 at 45,401,570 bp) were located in the upstream region of *SMIM20* and overlapped a TFBS. *SEL1L3* was previously identified as a candidate gene in a QTL region with effects on bovine milk protein composition¹⁹. Instead, to the best of our knowledge, this study is the first to propose *SMIM20* as a candidate gene for milk composition. However, the functional link between these genes and milk composition has yet to be established, as we found both genes to be endometrium-specific and underexpressed in mammary glands compared to the 90 other tissues investigated.

For milk mineral composition, the two best functional candidate genes highlighted by our analysis were *SLC37A1* (*solute carrier family 37, member A1*) and *ANKH* (*inorganic pyrophosphate transport regulator*). These two genes were previously found to be overexpressed in mammary glands relative to 17 other types of tissue³¹, and here we found both among the top 10% of mammary-gland-specific genes (Pt = 3.4.10⁻²⁶ and 6.6.10⁻¹⁶, respectively). This suggests that the main function of these genes occurs in epithelial cells of this tissue. *SLC37A1* and *ANKH* both encode transmembrane proteins involved in ion transport and have been found to play a role in inorganic anion transport²⁰. Variants located in both genes explained a large degree of the genetic variation in milk mineral content. Earlier studies also linked *SLC37A1* with mineral content, in particular that of P⁹,¹², while both genes have been implicated in determining milk protein composition in Holstein, Montbéliarde, and Normande cows¹⁹ and cheese-making traits in the same Montbéliarde cows we analyzed here²⁰. These results are consistent with what is known about the association of minerals with casein molecules in micelles; milk mineral composition is strongly related to milk protein composition² and therefore to cheese-making traits³,⁴.

In *SLC37A1*, 7, 31, 9, and 7 variants were located in the confidence intervals of QTL for P, K, Mg, and Na, respectively. These variants (32 of which were unique) were overwhelmingly located in introns of the gene, in a 12.3-kbp-region from 142,826,156 to 142,838,477 bp. Among these 32 variants, we were not able to distinguish the best candidate based on annotation, but two variants, 814 bp apart, were particularly highly ranked for all traits. Specifically, rs109459130, at 142,834,737 bp, was ranked 1st, 2nd, 1st, and 3rd for P, K, Mg, and Na, respectively, while rs109717634, at 142,835,551 bp, was ranked 2nd, 1st, 3rd, and 1st for P, K, Mg, and Na, respectively. Of these, the first (rs109459130) appeared to be the better candidate because it was top-ranked for the most affected trait, P, and because it was imputed (R² = 0.85), while the second variant was an HD SNP (R² = 0.997).

In the region of the *ANKH* gene, we found 81, 65, 69, and 21 variants in the confidence intervals of QTL for citrate, Mg, K, and Ca, respectively. These represented 82 unique variants, all located either in intronic regions of the *ANKH* gene (59 variants between 58,344,839 and 58,388,849 bp) or in the lncRNA *ENSBTAG00000048498* (23 variants between 58,185,895 and 58,212,187 bp). The same variant was ranked first for both Mg and Ca. Two top-1 variants were located in the lncRNA *ENSBTAG00000048498*—rs109956167 (at 58,204,929 bp) for citrate and rs110048176 (at 58,189,663 bp) for Mg and Ca—while the top-ranked variant for K was located in *ANKH* (rs134021638 at 58,386,888 bp). We propose rs109956167 as the best candidate causal variant because i) it was ranked 1st for the most-affected trait (citrate), ii) it was imputed, in contrast to rs110048176 which was from the HD chip, and iii) it was the only variant present in the top 10 for all traits affected by this region (4th for Mg, 6th for Ca, and 9th for K). We hypothesize that the lncRNA *ENSBTAG00000048498*, which is 83 kbp downstream of the closest gene, *ANKH*, may affect milk mineral content by regulating the expression of *ANKH* and thus the amount of protein available for ion transport.

When we examined the best candidate variants in the regions of *SLC37A1* (rs109459130) and *ANKH* (rs109956167), we found that alleles of these genes did not have significant interactions with respect to the six milk composition traits analyzed in this study and that, combined together, these two regions explained a large part of the genetic variance in milk mineral content. The allele that increased mineral content was the most frequent allele of *SLC37A1* (allele A, frequency = 0.58), while it was the least frequent allele in the *ANKH* region (allele C, frequency = 0.17).

Compared to previous studies, our GWAS-based investigation of milk mineral composition was conducted at the whole-genome level, after imputation with the most recent reference genome and using data from a large population of animals. This approach led to the identification of a large number of genomic regions that explain a great deal of the genetic variance of the traits studied here. Moreover, post-GWAS investigations conducted using different annotation datasets enabled the identification of candidate causal genes and the prioritization of candidate variants in these genomic regions. In most situations, exonic variants were not proposed as the best candidates, although they cannot be excluded in four QTL regions, including *GHR*. In contrast, the best candidate variants often were putative regulatory variants. Although the functional causative effect of the candidate variants remains to be demonstrated, this study highlights a 'short list' of candidate variants, in particular rs109459130 (*SLC37A1*) and rs109956167 (*ANKH*), whose favorable alleles can be feasibly selected to increase the mineral content (Ca, P, Mg, and K) of milk and therefore improve both its nutritional and cheese-making qualities in Montbéliarde cows.

## Methods

**Ethics statement.** Milk samples were analyzed with MIR spectrometry during routine milk recording in commercial herds of Montbéliarde cows in the Franche-Comté region (France). We did not perform any experiments on animals and no ethical approval was required.
On BTA6 at ~45.3 Mbp, we identified a QTL with very strong significant effects on K and to a lesser extent on P, Na, and Ca. Within the confidence interval of this QTL, we found 80 different potential variants, of which the majority (69) were intergenic, 9 were located in an intronic region of *SEL1L3* (*SEL1L family member 3*), and 2 were in the upstream region of *SMIM20* (*small integral membrane protein 20*). However, when we examined only the top 10 variants for the different traits, the number of distinct variants shrank to 18 (9 intergenic, 7 in *SEL1L3*, and 2 in *SMIM20*). Of these, two (rs108972810 at 45,401,485 bp and rs136498639 at 45,401,570 bp) were located in the upstream region of *SMIM20* and overlapped a TFBS. *SEL1L3* was previously identified as a candidate gene in a QTL region with effects on bovine milk protein composition¹⁹. Instead, to the best of our knowledge, this study is the first to propose *SMIM20* as a candidate gene for milk composition. However, the functional link between these genes and milk composition has yet to be established, as we found both genes to be endometrium-specific and underexpressed in mammary glands compared to the 90 other tissues investigated.

For milk mineral composition, the two best functional candidate genes highlighted by our analysis were *SLC37A1* (*solute carrier family 37, member A1*) and *ANKH* (*inorganic pyrophosphate transport regulator*). These two genes were previously found to be overexpressed in mammary glands relative to 17 other types of tissue³¹, and here we found both among the top 10% of mammary-gland-specific genes (Pt = 3.4.10⁻²⁶ and 6.6.10⁻¹⁶, respectively). This suggests that the main function of these genes occurs in epithelial cells of this tissue. *SLC37A1* and *ANKH* both encode transmembrane proteins involved in ion transport and have been found to play a role in inorganic anion transport²⁰. Variants located in both genes explained a large degree of the genetic variation in milk mineral content. Earlier studies also linked *SLC37A1* with mineral content, in particular that of P⁹,¹², while both genes have been implicated in determining milk protein composition in Holstein, Montbéliarde, and Normande cows¹⁹ and cheese-making traits in the same Montbéliarde cows we analyzed here²⁰. These results are consistent with what is known about the association of minerals with casein molecules in micelles; milk mineral composition is strongly related to milk protein composition² and therefore to cheese-making traits³,⁴.

In *SLC37A1*, 7, 31, 9, and 7 variants were located in the confidence intervals of QTL for P, K, Mg, and Na, respectively. These variants (32 of which were unique) were overwhelmingly located in introns of the gene, in a 12.3-kbp-region from 142,826,156 to 142,838,477 bp. Among these 32 variants, we were not able to distinguish the best candidate based on annotation, but two variants, 814 bp apart, were particularly highly ranked for all traits. Specifically, rs109459130, at 142,834,737 bp, was ranked 1st, 2nd, 1st, and 3rd for P, K, Mg, and Na, respectively, while rs109717634, at 142,835,551 bp, was ranked 2nd, 1st, 3rd, and 1st for P, K, Mg, and Na, respectively. Of these, the first (rs109459130) appeared to be the better candidate because it was top-ranked for the most affected trait, P, and because it was imputed (R² = 0.85), while the second variant was an HD SNP (R² = 0.997).

In the region of the *ANKH* gene, we found 81, 65, 69, and 21 variants in the confidence intervals of QTL for citrate, Mg, K, and Ca, respectively. These represented 82 unique variants, all located either in intronic regions of the *ANKH* gene (59 variants between 58,344,839 and 58,388,849 bp) or in the lncRNA *ENSBTAG00000048498* (23 variants between 58,185,895 and 58,212,187 bp). The same variant was ranked first for both Mg and Ca. Two top-1 variants were located in the lncRNA *ENSBTAG00000048498*—rs109956167 (at 58,204,929 bp) for citrate and rs110048176 (at 58,189,663 bp) for Mg and Ca—while the top-ranked variant for K was located in *ANKH* (rs134021638 at 58,386,888 bp). We propose rs109956167 as the best candidate causal variant because i) it was ranked 1st for the most-affected trait (citrate), ii) it was imputed, in contrast to rs110048176 which was from the HD chip, and iii) it was the only variant present in the top 10 for all traits affected by this region (4th for Mg, 6th for Ca, and 9th for K). We hypothesize that the lncRNA *ENSBTAG00000048498*, which is 83 kbp downstream of the closest gene, *ANKH*, may affect milk mineral content by regulating the expression of *ANKH* and thus the amount of protein available for ion transport.

When we examined the best candidate variants in the regions of *SLC37A1* (rs109459130) and *ANKH* (rs109956167), we found that alleles of these genes did not have significant interactions with respect to the six milk composition traits analyzed in this study and that, combined together, these two regions explained a large part of the genetic variance in milk mineral content. The allele that increased mineral content was the most frequent allele of *SLC37A1* (allele A, frequency = 0.58), while it was the least frequent allele in the *ANKH* region (allele C, frequency = 0.17).

Compared to previous studies, our GWAS-based investigation of milk mineral composition was conducted at the whole-genome level, after imputation with the most recent reference genome and using data from a large population of animals. This approach led to the identification of a large number of genomic regions that explain a great deal of the genetic variance of the traits studied here. Moreover, post-GWAS investigations conducted using different annotation datasets enabled the identification of candidate causal genes and the prioritization of candidate variants in these genomic regions. In most situations, exonic variants were not proposed as the best candidates, although they cannot be excluded in four QTL regions, including *GHR*. In contrast, the best candidate variants often were putative regulatory variants. Although the functional causative effect of the candidate variants remains to be demonstrated, this study highlights a 'short list' of candidate variants, in particular rs109459130 (*SLC37A1*) and rs109956167 (*ANKH*), whose favorable alleles can be feasibly selected to increase the mineral content (Ca, P, Mg, and K) of milk and therefore improve both its nutritional and cheese-making qualities in Montbéliarde cows.

## Methods

**Ethics statement.** Milk samples were analyzed with MIR spectrometry during routine milk recording in commercial herds of Montbéliarde cows in the Franche-Comté region (France). We did not perform any experiments on animals and no ethical approval was required.
**Animals and phenotypes.** The original dataset was generated by the From’MIR project and is described in detail by Sanchez et al.¹¹,²⁰. It comprised 6,670,769 mid-infrared (MIR) spectra of milk samples from 410,622 Montbéliarde cows. The concentrations of five milk minerals (Ca, P, Mg, K, and Na) and citrate were predicted from MIR spectra using equations developed by the Optimir project¹⁴,¹⁵ (Table 1). As previously described¹¹,²⁰, to ensure that the dataset was homogeneous (and less computationally demanding), we retained only the first-lactation records with at least seven test-day records per cow for estimating variance components (1,100,238 test-day records from 126,873 cows¹¹). To estimate environmental effects used to derive the phenotypes included in GWAS, at least three test-day records per cow were required, corresponding to a dataset of 1,442,371 test-day records from 189,817 cows²⁰.

**Estimation of genetic parameters.** First lactation data were analyzed using bivariate repeatability animal models applied to all pairs of traits. (Co)variance components were estimated using the AI-REML algorithm as implemented in Wombat software³², with the following linear animal model:

$$y = X\beta + Za + Wp + e_1, \qquad (1)$$

where y is the vector of test-day observations; a ~ N(0, A ⊗ G₀) is the vector of random additive genetic effects, p ~ N(0, I ⊗ P₀) is the vector of random permanent environmental effects, and e₁ ~ N(0, I ⊗ R₀) is the vector of random residual effects. X, Z, and W are incidence matrices, A is the relationship matrix among individuals calculated from the pedigree (traced back over four generations and containing 315,661 animals), and I the identity matrix. G₀, P₀ and R₀ are 2 × 2 matrices of additive genetic, permanent environmental, and residual variances-covariances, respectively. The β vector included the fixed effects of the herd x test-day x spectrometer combination, stage of lactation, and season of calving.

**Genotypes and imputation to whole-genome sequences.** A subset of 19,586 cows for which MIR spectra were available had been genotyped for the purpose of genomic selection using the BovineSNP50 (50 K, 6505 cows) or the EuroG10K BeadChip (Illumina Inc., San Diego, 13,081 cows); the latter contains generic SNPs and a research add-on for causal or predictive SNPs for traits of interest in cattle. Missing genotypes were imputed using FImpute software³³ for the 53,469 autosomal SNPs (50 K and “research” SNPs, 50 K+) that passed all quality control filters (individual call rate > 95%; SNP call rate > 90%; minor allele frequency (MAF) > 1% in at least one major French dairy cattle breed; genotype frequencies in HW equilibrium with P > 10⁻⁴). Whole-genome sequences (WGSs) were then imputed in two steps. In the first step, 777 K high-density and “research” SNPs (HD+) were imputed with FImpute software using a within-breed reference set of 522 Montbéliarde bulls genotyped with the Illumina BovineHD BeadChip (Illumina Inc., San Diego, CA)³⁴. Finally, allele dosages of the WGSs were imputed using Minimac software³⁵ and WGS variants of 1479 Bos taurus animals from the 7th run of the 1000 Bull Genomes Project, representing 17 cattle breeds and including 63 Montbéliarde bulls. This 2-step strategy was found to be the most accurate³⁶ (eg, Boowman & Veerkamp, 2014) because (1) the first step takes advantage of the high number of HD genotypes of influential bulls and is performed with very limited loss in accuracy³⁴ (Hozé et al., 2013); (2) the number of sequenced Montbéliarde bulls was too limited for an accurate imputation to sequence when used alone; and (3) linkage disequilibrium is partially conserved across breeds at the HD level (~ several kb), making use of close breeds sequence data very beneficial. WGS variants were selected following the protocol defined by the 1000 Bull Genomes consortium¹⁶,²³, as described in Bous-saha et al.³⁷. Short reads were filtered for quality and aligned to the ARS-UCD1.2 reference sequence¹⁷, and small genomic variations (SNPs and Indels) were detected using SAMtools 0.0.18³⁸. Raw variants were then filtered to produce a dataset of 25,050,323 variants. The precision of imputation from HD+ to WGS was assessed using R² values calculated with Minimac software³⁵. Only variants with R² ≥ 0.20 and MAF ≥ 0.005 were retained for association analyses, i.e. 12,907,802 variants, with a mean R² of 0.67 and MAF of 0.19.

**GWAS.** We performed single-trait association analyses between all 12,907,802 polymorphic variants and each of six milk mineral and citrate traits, described in Table 1. All association analyses were performed using the mlma option of GCTA software (version 1.24), which applies a mixed linear model that includes the variant to be tested³⁹:

$$y\mathbf{d} = \mathbf{1}m + \mathbf{x}_v\mathbf{b}_v + \mathbf{u} + \varepsilon, \qquad (2)$$

where yd is the vector of so-called yield deviations, i.e. test-day records adjusted for non-genetic effects with the mixed linear model (1) using Genekit software⁴⁰ and averaged per cow; m is the overall mean; bv is the additive fixed effect of the variant to be tested for association; xv is the vector of imputed allele dosages, ranging from 0 to 2; u ~ N(0, Gσ²) is the vector of random polygenic effects, with G the genomic relationship matrix (GRM) calculated using the HD SNP genotypes (which offer both high density and accuracy of imputation), and σ²u the polygenic variance, estimated based on the null model (yd = 1m + u + e) and then assumed as known while testing for the association between each variant and the trait of interest; and ε ~ N(0, Iσ²) is the vector of random residual effects, with I the identity matrix and σ² the residual variance. Association was tested using a t-statistic calculated by dividing the variant effect estimate by its standard error.

In order to correct for multiple testing, the Bonferroni correction was applied to take into account all 12.9 million independent tests. The 5% genome-wide threshold of significance therefore corresponded to a nominal P-value of 4 × 10⁻⁹ (-log₁₀(P) = 8.4) per test. When a given trait was significantly affected by multiple variants, variants that were located less than 2 million base-pairs (Mbp) apart were grouped together to define QTL,
**Animals and phenotypes.** The original dataset was generated by the From’MIR project and is described in detail by Sanchez et al.¹¹,²⁰. It comprised 6,670,769 mid-infrared (MIR) spectra of milk samples from 410,622 Montbéliarde cows. The concentrations of five milk minerals (Ca, P, Mg, K, and Na) and citrate were predicted from MIR spectra using equations developed by the Optimir project¹⁴,¹⁵ (Table 1). As previously described¹¹,²⁰, to ensure that the dataset was homogeneous (and less computationally demanding), we retained only the first-lactation records with at least seven test-day records per cow for estimating variance components (1,100,238 test-day records from 126,873 cows¹¹). To estimate environmental effects used to derive the phenotypes included in GWAS, at least three test-day records per cow were required, corresponding to a dataset of 1,442,371 test-day records from 189,817 cows²⁰.

**Estimation of genetic parameters.** First lactation data were analyzed using bivariate repeatability animal models applied to all pairs of traits. (Co)variance components were estimated using the AI-REML algorithm as implemented in Wombat software³², with the following linear animal model:

$$y = X\beta + Za + Wp + e_1, \qquad (1)$$

where y is the vector of test-day observations; a ~ N(0, A ⊗ G₀) is the vector of random additive genetic effects, p ~ N(0, I ⊗ P₀) is the vector of random permanent environmental effects, and e₁ ~ N(0, I ⊗ R₀) is the vector of random residual effects. X, Z, and W are incidence matrices, A is the relationship matrix among individuals calculated from the pedigree (traced back over four generations and containing 315,661 animals), and I the identity matrix. G₀, P₀ and R₀ are 2 × 2 matrices of additive genetic, permanent environmental, and residual variances-covariances, respectively. The β vector included the fixed effects of the herd x test-day x spectrometer combination, stage of lactation, and season of calving.

**Genotypes and imputation to whole-genome sequences.** A subset of 19,586 cows for which MIR spectra were available had been genotyped for the purpose of genomic selection using the BovineSNP50 (50 K, 6505 cows) or the EuroG10K BeadChip (Illumina Inc., San Diego, 13,081 cows); the latter contains generic SNPs and a research add-on for causal or predictive SNPs for traits of interest in cattle. Missing genotypes were imputed using FImpute software³³ for the 53,469 autosomal SNPs (50 K and “research” SNPs, 50 K+) that passed all quality control filters (individual call rate > 95%; SNP call rate > 90%; minor allele frequency (MAF) > 1% in at least one major French dairy cattle breed; genotype frequencies in HW equilibrium with P > 10⁻⁴). Whole-genome sequences (WGSs) were then imputed in two steps. In the first step, 777 K high-density and “research” SNPs (HD+) were imputed with FImpute software using a within-breed reference set of 522 Montbéliarde bulls genotyped with the Illumina BovineHD BeadChip (Illumina Inc., San Diego, CA)³⁴. Finally, allele dosages of the WGSs were imputed using Minimac software³⁵ and WGS variants of 1479 Bos taurus animals from the 7th run of the 1000 Bull Genomes Project, representing 17 cattle breeds and including 63 Montbéliarde bulls. This 2-step strategy was found to be the most accurate³⁶ (eg, Boowman & Veerkamp, 2014) because (1) the first step takes advantage of the high number of HD genotypes of influential bulls and is performed with very limited loss in accuracy³⁴ (Hozé et al., 2013); (2) the number of sequenced Montbéliarde bulls was too limited for an accurate imputation to sequence when used alone; and (3) linkage disequilibrium is partially conserved across breeds at the HD level (~ several kb), making use of close breeds sequence data very beneficial. WGS variants were selected following the protocol defined by the 1000 Bull Genomes consortium¹⁶,²³, as described in Bous-saha et al.³⁷. Short reads were filtered for quality and aligned to the ARS-UCD1.2 reference sequence¹⁷, and small genomic variations (SNPs and Indels) were detected using SAMtools 0.0.18³⁸. Raw variants were then filtered to produce a dataset of 25,050,323 variants. The precision of imputation from HD+ to WGS was assessed using R² values calculated with Minimac software³⁵. Only variants with R² ≥ 0.20 and MAF ≥ 0.005 were retained for association analyses, i.e. 12,907,802 variants, with a mean R² of 0.67 and MAF of 0.19.

**GWAS.** We performed single-trait association analyses between all 12,907,802 polymorphic variants and each of six milk mineral and citrate traits, described in Table 1. All association analyses were performed using the mlma option of GCTA software (version 1.24), which applies a mixed linear model that includes the variant to be tested³⁹:

$$y\mathbf{d} = \mathbf{1}m + \mathbf{x}_v\mathbf{b}_v + \mathbf{u} + \varepsilon, \qquad (2)$$

where yd is the vector of so-called yield deviations, i.e. test-day records adjusted for non-genetic effects with the mixed linear model (1) using Genekit software⁴⁰ and averaged per cow; m is the overall mean; bv is the additive fixed effect of the variant to be tested for association; xv is the vector of imputed allele dosages, ranging from 0 to 2; u ~ N(0, Gσ²) is the vector of random polygenic effects, with G the genomic relationship matrix (GRM) calculated using the HD SNP genotypes (which offer both high density and accuracy of imputation), and σ² the polygenic variance, estimated based on the null model (yd = 1m + u + e) and then assumed as known while testing for the association between each variant and the trait of interest; and ε ~ N(0, Iσ²) is the vector of random residual effects, with I the identity matrix and σ² the residual variance. Association was tested using a t-statistic calculated by dividing the variant effect estimate by its standard error.

In order to correct for multiple testing, the Bonferroni correction was applied to take into account all 12.9 million independent tests. The 5% genome-wide threshold of significance therefore corresponded to a nominal P-value of 4 × 10⁻⁹ (-log₁₀(P) = 8.4) per test. When a given trait was significantly affected by multiple variants, variants that were located less than 2 million base-pairs (Mbp) apart were grouped together to define QTL,
considering the variants belonged to the same QTL region. The bounds of the confidence intervals (CIs) of each region were then determined based on the positions of variants that were included in the upper third of the QTL peak.

In regions where multiple neighboring QTL were identified (BTA1 and BTA20), conditional analyses were carried out using the cojo option of GCTA³⁹ in order to conclude if multiple significant variants in a genomic region were due to LD with the same causal mutation or to the presence of multiple causal mutations. Association analyses were performed by including in the model the most significant variant as a fixed effect and by testing all variants in these neighboring QTL that were not in strong LD with the conditional variant (r² < 0.9).

**Functional annotations.** Genomic regions and variants were annotated with FAANGMine v1.1 (https://faangmine.elsiklab.missouri.edu/), developed by the Functional Annotation of ANimal Genomes initiative¹⁸ and which integrates the ARS-UCD1.2 bovine reference genome with a variety of external data sources, including RefSeq from NCBI (https://www.ncbi.nlm.nih.gov) and Ensembl (https://www.ensembl.org) gene sets.

The ability of genetic variants to alter transcription factor binding sites (TFBSs) was predicted with a custom script that used TFBS models from the JASPAR (JASPAR CORE 2018 collection⁴¹), HOCOMOCO (version v10⁴²), and TRANSFAC (version v3.2 public⁴³) databases. These databases contain curated sets of transcription factor binding models represented as Position Weight Matrices (PWM), which are derived from published collections of experimentally defined eukaryote TFBSs. Only vertebrate PWMs were downloaded for use in our study.

Gene overexpression or specificity in different tissues was determined using gene expression patterns of 24,616 genes (Ensembl release 94) available in the Cattle Gene Atlas, which contains 723 RNA-seq datasets representing 91 tissues and cell types, classified into 17 biological categories (http://cattlegeneatlas.roslin.ed.ac.uk/). To assess the expression specificity of each gene in a given type of tissue (by excluding tissues in the same biological category), we applied the following linear model as described in Fang et al.⁴⁴:

$$ye = 1me + x_t b_t + zc + ee \quad (3)$$

where ye is the vector of expression level in the tissues, assessed by the scaled log₂FPKM (Fragments Per Kilobase per Million mapped reads); me is the overall mean; x_t is the vector of the variable with value 1 for samples of the tested tissue and -1 for samples outside the same category; b_t is the corresponding tissue effect; z is the incidence matrix related to the corresponding covariables effects c, including age, sex and study effects; ee is the residual effect. Model (3) was implemented adapting R scripts available on the Cattle Gene Atlas website. For each gene, a t-statistic was computed by dividing the tissue effect by its standard error. A gene was considered to be overexpressed in a tissue if the probability associated with the t-statistic (Pt) was lower than 10⁻⁴. To determine the tissue specificity of genes, we ranked the genes in each type of tissue by their t-statistics; the top 10% were considered to be tissue-specific.

**SLC37A1-ANKH genotype interactions.** We tested putative interaction effects between the genes SLC37A1 and ANKH on milk mineral and citrate content with the following mixed linear model:

$$yd = 1mi + Mg_{SLC37A1} + Ng_{ANKH} + Og_{SLC37A1} \times g_{ANKH} + Ps + ei, \quad (4)$$

where yd as defined in (2); mi is the overall mean; gSLC37A1 and gANKH are the fixed effects of the genotypes of the best candidate variant in the SLC37A1 and ANKH regions, respectively; gSLC37A1 x gANKH represents the interaction between the two genotypes; M, N, and O are incidence matrices related to the individual effects of the SLC37A1 and ANKH genotypes and their interaction, respectively; s is the vector of random sire effects and P the corresponding incidence matrix; and ei is the vector of random residual effects. All effects were tested using t-statistics computed in the MIXED procedure of SAS software.

## Data availability

The data (genotypes and phenotypes) that enabled the findings of this study were made available by U MOTEST, CEL25-90, and GENIATEST. However, restrictions apply to the availability of these commercial data: they were used under license for the current study, and are not publicly available.

Received: 3 November 2020; Accepted: 23 March 2021

Published online: 06 April 2021

## References

1. Zamberlin, S., Antunac, N., Havranek, J. & Samarzija, D. Mineral elements in milk and dairy products. *Mljekarstvo* 62, 111–125 (2012).

2. Gaucheron, F. The minerals of milk. *Reprod. Nutr. Dev.* 45, 473–483 (2005).

3. Toffanin, V., De Marchi, M., Lopez-Villalobos, N. & Cassandro, M. Effectiveness of mid-infrared spectroscopy for prediction of the contents of calcium and phosphorus, and titratable acidity of milk and their relationship with milk quality and coagulation properties. *Int. Dairy J.* 41, 68–73 (2015).

4. Malacarne, M. et al. Influence of micellar calcium and phosphorus on rennet coagulation properties of cows milk. *J. Dairy Res.* 81, 129–136 (2014).

5. Toffanin, V. et al. Genetic parameters for milk mineral content and acidity predicted by mid-infrared spectroscopy in Holstein-Friesian cows. *Animal* 9, 775–780 (2015).

6. van Hulzen, K. J. E., Sprong, R. C., van der Meer, R. & van Arendonk, J. A. M. Genetic and nongenetic variation in concentration of selenium, calcium, potassium, zinc, magnesium, and phosphorus in milk of Dutch Holstein-Friesian cows. *J. Dairy Sci.* 92, 5754–5759 (2009).
considering the variants belonged to the same QTL region. The bounds of the confidence intervals (CIs) of each region were then determined based on the positions of variants that were included in the upper third of the QTL peak.

In regions where multiple neighboring QTL were identified (BTA1 and BTA20), conditional analyses were carried out using the cojo option of GCTA³⁹ in order to conclude if multiple significant variants in a genomic region were due to LD with the same causal mutation or to the presence of multiple causal mutations. Association analyses were performed by including in the model the most significant variant as a fixed effect and by testing all variants in these neighboring QTL that were not in strong LD with the conditional variant (r² < 0.9).

**Functional annotations.** Genomic regions and variants were annotated with FAANGMine v1.1 (https://faangmine.elsiklab.missouri.edu/), developed by the Functional Annotation of ANimal Genomes initiative¹⁸ and which integrates the ARS-UCD1.2 bovine reference genome with a variety of external data sources, including RefSeq from NCBI (https://www.ncbi.nlm.nih.gov) and Ensembl (https://www.ensembl.org) gene sets.

The ability of genetic variants to alter transcription factor binding sites (TFBSs) was predicted with a custom script that used TFBS models from the JASPAR (JASPAR CORE 2018 collection⁴¹), HOCOMOCO (version v10⁴²), and TRANSFAC (version v3.2 public⁴³) databases. These databases contain curated sets of transcription factor binding models represented as Position Weight Matrices (PWM), which are derived from published collections of experimentally defined eukaryote TFBSs. Only vertebrate PWMs were downloaded for use in our study.

Gene overexpression or specificity in different tissues was determined using gene expression patterns of 24,616 genes (Ensembl release 94) available in the Cattle Gene Atlas, which contains 723 RNA-seq datasets representing 91 tissues and cell types, classified into 17 biological categories (http://cattlegeneatlas.roslin.ed.ac.uk/). To assess the expression specificity of each gene in a given type of tissue (by excluding tissues in the same biological category), we applied the following linear model as described in Fang et al.⁴⁴:

$$ye = 1me + x_t b_t + zc + ee \quad (3)$$

where ye is the vector of expression level in the tissues, assessed by the scaled log₂FPKM (Fragments Per Kilobase per Million mapped reads); me is the overall mean; x_t is the vector of the variable with value 1 for samples of the tested tissue and -1 for samples outside the same category; b_t is the corresponding tissue effect; z is the incidence matrix related to the corresponding covariables effects c, including age, sex and study effects; ee is the residual effect. Model (3) was implemented adapting R scripts available on the Cattle Gene Atlas website. For each gene, a t-statistic was computed by dividing the tissue effect by its standard error. A gene was considered to be overexpressed in a tissue if the probability associated with the t-statistic (Pt) was lower than 10⁻⁴. To determine the tissue specificity of genes, we ranked the genes in each type of tissue by their t-statistics; the top 10% were considered to be tissue-specific.

**SLC37A1-ANKH genotype interactions.** We tested putative interaction effects between the genes SLC37A1 and ANKH on milk mineral and citrate content with the following mixed linear model:

$$yd = 1mi + Mg_{SLC37A1} + Ng_{ANKH} + Og_{SLC37A1} \times g_{ANKH} + Ps + ei, \quad (4)$$

where yd as defined in (2); mi is the overall mean; gSLC37A1 and gANKH are the fixed effects of the genotypes of the best candidate variant in the SLC37A1 and ANKH regions, respectively; gSLC37A1 x gANKH represents the interaction between the two genotypes; M, N, and O are incidence matrices related to the individual effects of the SLC37A1 and ANKH genotypes and their interaction, respectively; s is the vector of random sire effects and P the corresponding incidence matrix; and ei is the vector of random residual effects. All effects were tested using t-statistics computed in the MIXED procedure of SAS software.

## Data availability

The data (genotypes and phenotypes) that enabled the findings of this study were made available by U MOTEST, CEL25-90, and GENIATEST. However, restrictions apply to the availability of these commercial data: they were used under license for the current study, and are not publicly available.

Received: 3 November 2020; Accepted: 23 March 2021

Published online: 06 April 2021

## References

1. Zamberlin, S., Antunac, N., Havranek, J. & Samarzija, D. Mineral elements in milk and dairy products. *Mljekarstvo* 62, 111–125 (2012).

2. Gaucheron, F. The minerals of milk. *Reprod. Nutr. Dev.* 45, 473–483 (2005).

3. Toffanin, V., De Marchi, M., Lopez-Villalobos, N. & Cassandro, M. Effectiveness of mid-infrared spectroscopy for prediction of the contents of calcium and phosphorus, and titratable acidity of milk and their relationship with milk quality and coagulation properties. *Int. Dairy J.* 41, 68–73 (2015).

4. Malacarne, M. et al. Influence of micellar calcium and phosphorus on rennet coagulation properties of cows milk. *J. Dairy Res.* 81, 129–136 (2014).

5. Toffanin, V. et al. Genetic parameters for milk mineral content and acidity predicted by mid-infrared spectroscopy in Holstein-Friesian cows. *Animal* 9, 775–780 (2015).

6. van Hulzen, K. J. E., Sprong, R. C., van der Meer, R. & van Arendonk, J. A. M. Genetic and nongenetic variation in concentration of selenium, calcium, potassium, zinc, magnesium, and phosphorus in milk of Dutch Holstein-Friesian cows. *J. Dairy Sci.* 92, 5754–5759 (2009).
7. Visentin, G. et al. Genetic (co)variances between milk mineral concentration and chemical composition in lactating Holstein-Friesian dairy cows. *Animal* **13**, 477–486 (2019).

8. Denholm, S. et al. Phenotypic and genetic analysis of milk and serum element concentrations in dairy cows. *J. Dairy Sci.* **102**, 11180–11192 (2019).

9. Buitenhuis, B., Poulsen, N., Larsen, L. & Sehested, J. Estimation of genetic parameters and detection of quantitative trait loci for minerals in Danish Holstein and Danish Jersey milk. *BMC Genet.* **16**, 52 (2015).

10. Bonfatti, V., Vicario, D., Lugo, A. & Carnier, R. Genetic parameters of measures and population-wide infrared predictions of 92 traits describing the fine composition and technological properties of milk in Italian Simmental cattle. *J. Dairy Sci.* **100**, 5526–5540 (2017).

11. Sanchez, M. P. et al. 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.* **101**, 10048–10061 (2018).

12. Kemper, K. E. et al. Leveraging genetically simple traits to identify small-effect variants for complex phenotypes. *BMC Genom.* **17**, 1–9 (2016).

13. De Marchi, M., Toffanin, V., Cassandro, M. & Penasa, M. Invited review: Mid-infrared spectroscopy as phenotyping tool for milk traits. *J. Dairy Sci.* **97**, 1171–1186 (2014).

14. Gengler, N. et al. Capitalizing on fine milk composition for breeding and management of dairy cows. *J. Dairy Sci.* **99**, 4071–4079 (2016).

15. Soyeurt, H. et al. Potential estimation of major mineral contents in cow milk using mid-infrared spectrometry. *J. Dairy Sci.* **92**, 2444–2454 (2009).

16. Bouwman, A. C. et al. Meta-analysis of genome-wide association studies for cattle stature identifies common genes that regulate body size in mammals. *Nat. Genet.* **50**, 362–367 (2018).

17. Rosen, B. D. et al. De novo assembly of the cattle reference genome with single-molecule sequencing. *GigaScience* **9**, 2giaa021 (2020).

18. Giuffra, E., Tuggle, C.K., & Consortium, FAANG. Functional annotation of animal genomes (FAANG): Current achievements and roadmap. *Annu. Rev. Anim. Biosci.* **7**, 65–88 (2019).

19. Sanchez, M. P. et al. 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.* **49**, 68 (2017).

20. Sanchez, M. P. et al. 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.* **51**, 34 (2019).

21. Albert, F. & Kruglyak, L. The role of regulatory variation in complex traits and disease. *Nat. Rev. Genet.* **16**, 197–212 (2015).

22. Littlejohn, M. et al. Expression variants of the lipogenic AGPAT6 gene affect diverse milk composition phenotypes in Bos taurus. *PLoS ONE* **9**, e85757 (2014).

23. Daetwyler, H. D. et al. Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle. *Nat. Genet.* **46**, 858–867 (2014).

24. Pausch, H. et al. Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution. *BMC Genom.* **18**, 1–11 (2017).

25. 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.* **40**, 127–133 (2009).

26. Boichard, D. et al. Genetic determinism of milk composition in fatty acids and proteins in ruminants, and selection potential. *INRAE Prod. Anim.* **27**, 283–298 (2014).

27. Knutsen, T. et al. Unravelling genetic variation underlying de novo-synthesis of bovine milk fatty acids. *Sci. Rep.* **8**, 2179 (2018).

28. Schoonjans, K., Staels, B. & Auwerx, J. The peroxisome proliferator activated receptors (PPARs) and their effects on lipid metabolism and adipocyte differentiation. *Biochim. Biophys. Acta* **1302**, 93–109 (1996).

29. Fougère, H. & Bernard, L. Effect of diets supplemented with starch and corn oil, marine algae, or hydrogenated palm oil on mammary lipogenic gene expression in cows and goats: A comparative study. *J. Dairy Sci.* **102**, 768–779 (2019).

30. Bernard, L., Toral, P. G. & Chilliard, Y. Comparison of mammary lipid metabolism in dairy cows and goats fed diets supplemented with starch, plant oil, or fish oil. *J. Dairy Sci.* **100**, 9338–9351 (2017).

31. Chamberlain, A. et al. Extensive variation between tissues in allele specific expression in an outbred mammal. *BMC Genom.* **16**, 993 (2015).

32. Meyer, K. WOMBAT—A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML). *J. Zhejiang Univ-Sc.* **B** **8**, 815–821 (2007).

33. Sargolzaei, M., Chesnais, J. & Schenkel, F. A new approach for efficient genotype imputation using information from relatives. *BMC Genom.* **15**, 478 (2014).

34. Hoze, C. et al. High-density marker imputation accuracy in sixteen French cattle breeds. *Genet. Sel. Evol.* **45**, 33 (2013).

35. Howie, B., Fuchsberger, C., Stephens, M., Marchini, J. & Abecasis, G. R. Fast and accurate genotype imputation in genome-wide association studies through pre-phasing. *Nat. Genet.* **44**, 955–959 (2012).

36. Bouwman, A. C. & Veerkamp, R. F. Consequences of splitting whole-genome sequencing effort over multiple breeds on imputation accuracy. *BMC Genet.* **15**, 1–9 (2014).

37. Boussaha, M. et al. Construction of a large collection of small genome variations in French dairy and beef breeds using whole-genome sequences. *Genet. Sel. Evol.* **48**, 87 (2016).

38. Li, H. et al. The sequence alignment/map format and SAMtools. *Bioinformatics* **25**, 2078–2079 (2009).

39. Yang, J., Lee, S., Goddard, M. & Visscher, P. GCTA: A Tool for Genome-wide Complex Trait Analysis. *Am. J. Hum. Genet.* **88**, 76–82 (2011).

40. Ducrocq, V. *Genekit, BLUP software* (INRA, 2011).

41. Sandelin, A., Alkema, W., Engstrom, P., Wasserman, W. & Lenhard, B. JASPAR: An open-access database for eukaryotic transcription factor binding profiles. *Nucleic Acids Res.* **32**, D91–D94 (2004).

42. Khamis, A. et al. A novel method for improved accuracy of transcription factor binding site prediction. *Nucleic Acids Res.* **46**, e72 (2018).

43. Knüppel, A., Dietze, P., Lehnberg, W., Frech, K. & Wingender, E. TRANSFAC retrieval program: A network model database of eukaryotic transcription regulating sequences and proteins. *J. Comput. Biol.* **1**, 191–198 (1994).

44. Fang, L. et al. Comprehensive analyses of 723 transcriptomes enhance genetic and biological interpretations for complex traits in cattle. *Genome Res.* **30**, 790–801 (2020).

## Acknowledgements

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é, through the project *From’MIR*. Chip genotypes originated from the genomic selection program of the breeding company UMOTEST. The authors gratefully
7. Visentin, G. et al. Genetic (co)variances between milk mineral concentration and chemical composition in lactating Holstein-Friesian dairy cows. *Animal* **13**, 477–486 (2019).

8. Denholm, S. et al. Phenotypic and genetic analysis of milk and serum element concentrations in dairy cows. *J. Dairy Sci.* **102**, 11180–11192 (2019).

9. Buitenhuis, B., Poulsen, N., Larsen, L. & Sehested, J. Estimation of genetic parameters and detection of quantitative trait loci for minerals in Danish Holstein and Danish Jersey milk. *BMC Genet.* **16**, 52 (2015).

10. Bonfatti, V., Vicario, D., Lugo, A. & Carnier, R. Genetic parameters of measures and population-wide infrared predictions of 92 traits describing the fine composition and technological properties of milk in Italian Simmental cattle. *J. Dairy Sci.* **100**, 5526–5540 (2017).

11. Sanchez, M. P. et al. 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.* **101**, 10048–10061 (2018).

12. Kemper, K. E. et al. Leveraging genetically simple traits to identify small-effect variants for complex phenotypes. *BMC Genom.* **17**, 1–9 (2016).

13. De Marchi, M., Toffanin, V., Cassandro, M. & Penasa, M. Invited review: Mid-infrared spectroscopy as phenotyping tool for milk traits. *J. Dairy Sci.* **97**, 1171–1186 (2014).

14. Gengler, N. et al. Capitalizing on fine milk composition for breeding and management of dairy cows. *J. Dairy Sci.* **99**, 4071–4079 (2016).

15. Soyeurt, H. et al. Potential estimation of major mineral contents in cow milk using mid-infrared spectrometry. *J. Dairy Sci.* **92**, 2444–2454 (2009).

16. Bouwman, A. C. et al. Meta-analysis of genome-wide association studies for cattle stature identifies common genes that regulate body size in mammals. *Nat. Genet.* **50**, 362–367 (2018).

17. Rosen, B. D. et al. De novo assembly of the cattle reference genome with single-molecule sequencing. *GigaScience* **9**, 2giaa021 (2020).

18. Giuffra, E., Tuggle, C.K., & Consortium, FAANG. Functional annotation of animal genomes (FAANG): Current achievements and roadmap. *Annu. Rev. Anim. Biosci.* **7**, 65–88 (2019).

19. Sanchez, M. P. et al. 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.* **49**, 68 (2017).

20. Sanchez, M. P. et al. 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.* **51**, 34 (2019).

21. Albert, F. & Kruglyak, L. The role of regulatory variation in complex traits and disease. *Nat. Rev. Genet.* **16**, 197–212 (2015).

22. Littlejohn, M. et al. Expression variants of the lipogenic AGPAT6 gene affect diverse milk composition phenotypes in Bos taurus. *PLoS ONE* **9**, e85757 (2014).

23. Daetwyler, H. D. et al. Whole-genome sequencing of 234 bulls facilitates mapping of monogenic and complex traits in cattle. *Nat. Genet.* **46**, 858–867 (2014).

24. Pausch, H. et al. Meta-analysis of sequence-based association studies across three cattle breeds reveals 25 QTL for fat and protein percentages in milk at nucleotide resolution. *BMC Genom.* **18**, 1–11 (2017).

25. 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.* **40**, 127–133 (2009).

26. Boichard, D. et al. Genetic determinism of milk composition in fatty acids and proteins in ruminants, and selection potential. *INRAE Prod. Anim.* **27**, 283–298 (2014).

27. Knutsen, T. et al. Unravelling genetic variation underlying de novo-synthesis of bovine milk fatty acids. *Sci. Rep.* **8**, 2179 (2018).

28. Schoonjans, K., Staels, B. & Auwerx, J. The peroxisome proliferator activated receptors (PPARs) and their effects on lipid metabolism and adipocyte differentiation. *Biochim. Biophys. Acta* **1302**, 93–109 (1996).

29. Fougère, H. & Bernard, L. Effect of diets supplemented with starch and corn oil, marine algae, or hydrogenated palm oil on mammary lipogenic gene expression in cows and goats: A comparative study. *J. Dairy Sci.* **102**, 768–779 (2019).

30. Bernard, L., Toral, P. G. & Chilliard, Y. Comparison of mammary lipid metabolism in dairy cows and goats fed diets supplemented with starch, plant oil, or fish oil. *J. Dairy Sci.* **100**, 9338–9351 (2017).

31. Chamberlain, A. et al. Extensive variation between tissues in allele specific expression in an outbred mammal. *BMC Genom.* **16**, 993 (2015).

32. Meyer, K. WOMBAT—A tool for mixed model analyses in quantitative genetics by restricted maximum likelihood (REML). *J. Zhejiang Univ-Sc.* **B** **8**, 815–821 (2007).

33. Sargolzaei, M., Chesnais, J. & Schenkel, F. A new approach for efficient genotype imputation using information from relatives. *BMC Genom.* **15**, 478 (2014).

34. Hoze, C. et al. High-density marker imputation accuracy in sixteen French cattle breeds. *Genet. Sel. Evol.* **45**, 33 (2013).

35. Howie, B., Fuchsberger, C., Stephens, M., Marchini, J. & Abecasis, G. R. Fast and accurate genotype imputation in genome-wide association studies through pre-phasing. *Nat. Genet.* **44**, 955–959 (2012).

36. Bouwman, A. C. & Veerkamp, R. F. Consequences of splitting whole-genome sequencing effort over multiple breeds on imputation accuracy. *BMC Genet.* **15**, 1–9 (2014).

37. Boussaha, M. et al. Construction of a large collection of small genome variations in French dairy and beef breeds using whole-genome sequences. *Genet. Sel. Evol.* **48**, 87 (2016).

38. Li, H. et al. The sequence alignment/map format and SAMtools. *Bioinformatics* **25**, 2078–2079 (2009).

39. Yang, J., Lee, S., Goddard, M. & Visscher, P. GCTA: A Tool for Genome-wide Complex Trait Analysis. *Am. J. Hum. Genet.* **88**, 76–82 (2011).

40. Ducrocq, V. *Genekit, BLUP software* (INRA, 2011).

41. Sandelin, A., Alkema, W., Engstrom, P., Wasserman, W. & Lenhard, B. JASPAR: An open-access database for eukaryotic transcription factor binding profiles. *Nucleic Acids Res.* **32**, D91–D94 (2004).

42. Khamis, A. et al. A novel method for improved accuracy of transcription factor binding site prediction. *Nucleic Acids Res.* **46**, e72 (2018).

43. Knüppel, A., Dietze, P., Lehnberg, W., Frech, K. & Wingender, E. TRANSFAC retrieval program: A network model database of eukaryotic transcription regulating sequences and proteins. *J. Comput. Biol.* **1**, 191–198 (1994).

44. Fang, L. et al. Comprehensive analyses of 723 transcriptomes enhance genetic and biological interpretations for complex traits in cattle. *Genome Res.* **30**, 790–801 (2020).

## Acknowledgements

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é, through the project *From’MIR*. Chip genotypes originated from the genomic selection program of the breeding company U MOTEST. The authors gratefully
acknowledge the breeders who participated in the *From'MIR* project; colleagues from the Conseil-Elevage 25-90 who coordinated data collection; and the contribution of the 1000 Bull Genomes consortium.

**Author contributions**

M.P.S. estimated genetic parameters; performed imputation, GWAS, and post-GWAS analyses; and drafted the first version of the manuscript. D.R. and M.C. managed rSNP analyses. M.Bo. managed sequence analyses of the 1000 Bull Genomes Project. C.H. contributed to imputation analyses. M.Br., A.D.B., P.G., and D.B. designed and managed the *From'MIR* project. D.B. supervised all analyses. All authors have read and approved the manuscript.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-021-87078-1.

**Correspondence** and requests for materials should be addressed to M.-P.S.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021
acknowledge the breeders who participated in the *From'MIR* project; colleagues from the Conseil-Elevage 25-90 who coordinated data collection; and the contribution of the 1000 Bull Genomes consortium.

**Author contributions**

M.P.S. estimated genetic parameters; performed imputation, GWAS, and post-GWAS analyses; and drafted the first version of the manuscript. D.R. and M.C. managed rSNP analyses. M.Bo. managed sequence analyses of the 1000 Bull Genomes Project. C.H. contributed to imputation analyses. M.Br., A.D.B., P.G., and D.B. designed and managed the *From'MIR* project. D.B. supervised all analyses. All authors have read and approved the manuscript.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1038/s41598-021-87078-1.

**Correspondence** and requests for materials should be addressed to M.-P.S.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_8.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2021