Using QTL to Identify Genes and Pathways Underlying the Regulation and Production of Milk Components in Cattle

1. Introduction

The composition of milk is complex, featuring an emulsion of fat globules and a colloidal dispersion of casein micelles in an aqueous solution of lactose (and other carbohydrates), whey proteins, and minerals. Although milk from different species contains the same basic constituents, their proportions can vary greatly. In cattle, the typically average percentages (g/ 100 $\mathrm{g}$) of fat, caseins, whey proteins, and lactose are 3.9%, 2.6%, 0.6% and 4.6% respectively; in humans, the corresponding percentages are 4.5%, 0.4%, 0.5% and 7.1% [1]. Even more extreme differences can be seen in other species, with some seal species, for example, expressing little to no lactose, resulting in highly concentrated milk showing fat percentages of 50% [2].

Less extreme differences in milk composition are also visible within species. In many cases, the differences in composition among individuals are under partial genetic control. Regions of the genome where the genotypes of genetic variants are associated with phenotypes such as milk composition are known as quantitative trait loci (QTL). When they can be identified, the causative genes underlying these QTL can help elucidate the pathways involved in milk production. The aim of this review is to describe some of the major pathways required for milk production in terms of the QTL and genes that have helped to identify them, as well as to note some of the methods that can be used to identify causative genes underlying QTL. This will focus primarily on cattle (Bos taurus) as a model species; however, some comparisons with human milk composition are also presented.

2. QTL for Major Pathways Involved in Milk Production

Many QTL have been observed for both milk yield and milk composition traits. The genes attributed to these QTL have a variety of functions. Some, such as the hormones prolactin and growth hormone, and their associated signalling pathways, are involved in mammogenesis: the development of the mammary gland during puberty and pregnancy. Other pathways, such as fat and protein synthesis, and ion channels, are involved directly in milk production. The following sections list some of the pathways involved in mammogenesis and lactation, as identified by QTL for milk production where the candidate causal gene sits within those pathways.

2.1. Milk Proteins

QTL have been mapped to many milk proteins, i.e., those expressed directly into milk. In cattle, the largest proportion of milk protein content (80% [3]) consists of the four casein proteins, encoded by a cluster of genes mapping to BTA6: CSN1S1, encoding ${\alpha }_{S1}$-casein, CSN1S2 for ${\alpha }_{S2}$-casein, CSN2 for $\beta$-casein, and CSN3 for $\kappa$-casein. The casein proteins aggregate into miscelles in the milk, sequestering calcium phosphate as a nutrient source for the neonate. $\kappa$-casein seems to be particularly important for successful lactation, as knockout mice deficient in $\kappa$-casein fail to lactate, due to destabilisation of the casein miscelles [4]. Like $\beta$-lactoglobulin, the various casein proteins all show a range of coding variants, although, with the exception of CSN1S1*G[5], these have not been intrinsically linked to lower rates of gene expression like the $\beta$-lactoglobulin B variants discussed above. However, there is evidence that the SNP responsible for the A2 variant of $\beta$-casein is associated with both milk yield and protein yield [6], and QTL for milk protein phenotypes have also been detected at this locus in other studies [7,8]. These QTL can overlap with, but are not in linkage disequilibrium with, other nearby QTL assigned to the GC gene [7,9,10]: this gene encodes the protein Group-Specific Component (Vitamin D Binding), and maps around one megabase from the casein gene cluster. Outside the mammary gland, casein expression in human CD14⁺ monocytes has been observed to upregulate expression of of the cytokines granulocyte-macrophage colony-stimulating factor (GM-CSF), interleukin-1$\beta$ (IL-1$\beta$), and IL-6 via the p38-MAPK pathway [11,12]. This suggests a possible role for caseins in regulating the innate immune response in the neonate intestine following milk consumption.

The major whey protein in cattle, accounting for around 50% of whey protein [3] and 10% of total protein, is $\beta$-lactoglobulin, encoded by the gene PAEP on BTA11. A QTL has been detected mapping to this gene, with pleiotropic effects on milk fat yield, protein yield, and volume [7,13]. Other studies looking at individual milk proteins have also shown that genetic variants associated with low $\beta$-lactoglobulin concentrations yield higher concentrations of $\alpha$-, $\beta$-, and $\kappa$-caseins [14]. Variation in milk $\beta$-lactoglobulin concentrations are primarily driven by differences between the A and B protein variants [15,16], with the B variant showing lower expression than the A variant. Other variants mapping within the B variant background, such as B^* [17] and B′ [18], cause further reductions in $\beta$-lactoglobulin expression. Milk with low $\beta$-lactoglobulin concentrations has potential uses in infant formula (as human milk lacks this protein), and the associated higher concentrations of caseins would also be expected to provide better properties for cheese making.

The dominant whey protein in humans, and second in cattle (18.5% of total whey protein [3]), is $\alpha$-lactalbumin, encoded by the gene LALBA on BTA5. QTL associated with milk protein concentration have been identified at this locus in several dairy populations and breeds [19,20,21,22]. The $\alpha$-lactalbumin protein also comprises a regulatory subunit of lactose synthase complex, so it is no surprise that a QTL for milk lactose concentration has also been identified at this locus [23]. The importance of this is shown in $\alpha$-lactalbumin knockout mice, which produce highly viscous milk, with otherwise normal fat and protein composition, that cannot be extracted by the pups from the mammary gland [24]. It has been shown in vitro that multimeric $\alpha$-lactalbumin is cytotoxic to stem cells and transformed cell lines [25], suggesting that $\alpha$-lactalbumin may also have a protective function in either the mammary gland or neotate digestive system.

Another whey protein is the iron-binding antimicrobial protein Lactoferrin, encoded by the LTF gene on BTA22. In cattle, this protein exists in bovine milk at low but variable concentrations (average $115.4$ $\mu$$\mathrm{g}$/$\mathrm{m}$$\mathrm{l}$ in [26] but ranging from 31.8 to 485.6), but reaches much higher concentrations in human milk at around 2 $\mathrm{g}$/$\mathrm{l}$ in mature milk [27]. Interestingly, another antimicrobial protein, lysozyme, which is present at lower levels in human milk [3], is also relatively depleted in cattle [28], with concentrations of lysozyme in Bovidae reported to be as low as $1/{1000}^{\mathrm{th}}$ of other mammalian species [29]. Although lactoferrin expression can be induced by mastitic infection [30] it is also under partial genetic control [31], and genetic variants affecting expression have been identified at the LTF locus [32]. QTL at this locus have also been associated with casein number and lactose concentration [33].

A third antibacterial protein present in milk is lactoperoxidase, encoded by the LPO gene on BTA19. In contrast to lysozyme and lactoferrin, lactoperoxidase activity has been observed at around $20\times$ higher in bovine milk compared to that of humans [34]. A trans-eQTL for LPO, overlapping a QTL for milk protein concentration, has been identified on BTA20 at the locus of the C6 and C7 genes [35]: these two genes encode proteins that comprise part of the complement pathway in the innate immune system, suggesting that lactoperoxidase may be co-regulated with with this system.

2.2. Fat Synthesis Pathways

Fat is one of the major components of milk, forming membrane-bound droplets known as milk fat globules (MFG), mostly in the form of triglycerides. MFG membranes also contain a range of proteins, such as butyrophilin, adipophilin, mucin, lactadherin, lactoferrin, and xanthine oxidase [36]. Two of these, butyrophilin (BTN1A1) and xanthine oxidase (XOR), are required for enveloping the MFGs with the apical cell membrane and therefore for secretion of the MFG into milk [37]. Heterozygous knockout mice for XOR are unable to maintain lactation [38]. Many of the QTL for fat yield or milk concentration in cattle map to genes encoding fat synthesis or metabolic enzymes. One of the most prominent QTL detected in cattle [39] for milk volume, fat, and protein phenotypes in cattle maps to the DGAT1 locus on BTA14 [40]. This gene encodes the enzyme diacylglycerol O-acyltransferase 1, responsible for the final stage triglyceride fat synthesis [41]. The causative variant for this QTL is a non-conservative lysine to alanine substitution at position 232 [40]; more recently, this same variant has been shown to cause an expression QTL for the DGAT1 gene by disrupting an exon splice enhancer, which in turn alters the splicing efficiency of several introns in the transcript [42]. Another enzyme sitting earlier in the same trigyceride synthesis pathway is glycerol-3-phosphate acyltransferase 4, encoded by the GPAT4 gene (formerly known as AGPAT6) on BTA27. Like DGAT1, a highly pleiotropic QTL has been detected at this locus for many milk phenotypes, including volume, fat, protein, and lactose traits [43]. The enzyme lipin 1 sits in a related pathway, catalysing the conversion of phosphatidate into diacylglycerol, as well as functioning as a tanscriptional coactivator for genes involved in fatty acid oxidation [44]. QTL for milk protein and casein concentrations [33], as well as milk yield [45], have been detected at the LPIN1 locus, which encodes this enzyme, on BTA11 in cattle.

Before they can be assembled into triglycerides, fatty acids first need to be obtained either from the diet or from de novo synthesis. Several genes involved in this latter pathway have also shown QTL for milk-related phenotypes. The rate-limiting step in fatty acid synthesis is the carboxylation of acetyl-CoA to malonyl-CoA, catalysed by the enzyme Acetyl-CoA carboxylase (ACC). The alpha form of this enzyme is encoded by the gene ACACA on BTA19, and QTL for fatty acid composition [46] and somatic cell score (a proxy phenotype for mastitis) [33] have been mapped to this locus. Another important enzyme in the fatty acid synthesis pathway is fatty acid synthase (FAS), encoded by the gene FASN on BTA19, dimers of which are responsible for synthesising the saturated C16 fatty acid palmitic acid from acetyl-CoA and malonyl-CoA [47]. QTL attributed to this gene have been identified for milk fat yield [7], fat concentration [8,48], and fatty acid composition [49] in cattle.

The fatty acid synthase gene will produce only saturated fatty acids. To generate monounsaturated or polyunsaturated fatty acids, desaturase enzymes are required. One cluster of fatty acid desaturase genes on BTA29 includes the two gene FADS1 and FADS2; these enzymes are responsible for the final, rate-limiting step in omega-3 and -6 fatty acid synthesis [50,51]. Variants mapping in cattle to these two genes have been associated with concentrations of a range of polyunsaturated fatty acids in milk [52], and similar associations are also observed in human milk [53,54]. A third desaturase enzyme involved in fatty acid synthesis is stearoyl-CoA desaturase, encoded by the gene SCD on BTA26 and responsible for oxidising the C16 and C18 saturated fatty acid compounds palmitoyl- and stearoyl-CoA into the monounsaturated compounds palmitoleoyl- and oleoyl-CoA respectively [55]. This enzyme is also important in regulating metabolism: knockout mice, as well as exhibiting lower levels of tissue triglycerides and low-density lipoproteins [56], also show an increase in insulin signalling and glucose uptake in muscle tissue [57]. In cattle, a QTL assigned to SCD has been identified for milk fat yield [7,48].

2.3. Hormones and Signalling

Hormones, and the receptors and signalling pathways they activate, are important in most if not all biological functions, and milk production is no exception. For example, knocking out the SCD gene discussed in the previous section leads to an increase in tyrosine phosphorylation of the insulin receptor, with downstream effects on the PI3K-Akt signalling pathway, leading to increased levels of the glucose transporter GLUT4 in the plasma membrane and increased glucose uptake in muscle [57]. These knockout mice also showed lower levels of the hormone leptin, which is involved in regulating energy intake and partitioning. In cattle, leptin (encoded by the gene LEP on BTA4) has been associated with both milk yield and feed intake [58]. No milk phenotype QTL have been reported for the leptin receptor (LEPR on BTA3), though an association has been reported with body size [59]. In both humans [60] and mice [61], leptin resistance is associated with obesity. Another gene associated with obesity [62] is FTO, encoding the enzyme FTO alpha-ketoglutarate dependent dioxygenase. This enzyme is involved in DNA repair; specifically, it demethylates 3-methylthymidine [63]. It can also demethylate bases in RNA, including 3-methyluracil and 6-methyladenosine [64]. Via this latter mechanism, FTO can inhibit adipogenesis by demethylating cyclin A2 and cyclin-dependent kinase 2 mRNA, reducing expression of these genes and thereby prolonging the cell cycle [65]. In cattle, variants at the FTO locus on BTA18 have been associated with milk fat yield [66].

One hormone of particular importance for lactation is prolactin, a peptide hormone secreted by the anterior pituitary gland. In cattle, this peptide is encoded by the PRL gene on BTA23. The importance of this hormone in cattle was underlined by the discovery of a dominant missense mutation that caused, among other phenotypes, a failure to lactate [67]. Prolactin promotes development of the mammary gland during pregnancy in conjunction with progesterone to generate ductal branching and alveolar buds [68]. Prolactin is detected by cells using prolactin receptor (PRLR; BTA20) and acts via the PRLR/JAK2/STAT5 signalling pathway to promote mammary gland development and milk protein expression [69], as well as inducing expression of the enzyme UDP-glucose pyrophosphorylase 2 (UGP2) and the transporter UDP-galactose transporter 2 (SLC35A2), thereby promoting the synthesis of lactose [70]. In parallel, prolactin receptor also acts via the PI3K/Akt pathway to downregulate repressors of PRLR/JAK2/STAT5 signalling [71], with Akt1 also upregulating fat synthesis and glucose uptake into the cell for lactose production [72,73]. More recently, it has been shown that Akt signalling induces developing mammary epithelial cells to express prolactin, which in turn acts in an autocrine manner via STAT5 to cause terminal differentiation of the mammary epithelium [69,74]. Given the importance of these pathways, it is not surprising that QTL for milk yield have been widely identified at the PRLR locus [33,75,76], as well as for somatic cell score [76]. Likewise, many studies have reported milk related QTL at the STAT5 locus on BTA19 (genes STAT5A and STAT5B) [7,22,23,49,77].

Beyond the leptin and prolactin pathways, other hormones have also been linked to lactogenesis. One example is growth hormone, also known as somatotropin, which comprises part of the same family of protein hormones as prolactin [78]. Growth hormone acts in the lactating animal to partition nutrients and energy towards the mammary gland [79], though to be mediated by increased serum insulin-like growth factor (IGF-1) levels [80]. In the cow, growth hormone is encoded by the gene GH1 on BTA19, and its receptor, encoded by GHR, maps to BTA20. Both the GH1 locus [33], and especially the GHR locus [8,58,81,82,83], have been associated with milk, fat, and protein yield phenotypes in a range of cattle populations.

Another important family of signalling molecules is the interleukin family of cytokines, a group of primarily immunomodulatory proteins that operate in a paracrine or autocrine manner to affect cell growth and differentiation during immune responses [84]. The gene CSF2RB on BTA5 encodes the $\beta$-chain of GM-CSF (alos known as CSF2) receptor, and also forms a common subunit with the receptors for interleukin-3 (IL-3) and IL-5. QTL for milk, fat, and protein volume, as well as fat and protein concentrations, have been identified at this locus [7,20,85,86], with CSF2RB suggested as the most likely candidate causative gene, although the neighbouring genes NCF4 [85,86] and TST [8] have also been put forward as candidates. The GM-CSF receptor operates via the JAK/STAT signalling pathway, specifically, JAK2 and STAT5 [84], the same pathway as used by prolactin signalling.

In contrast to STAT5, which is linked to increased milk protein and lactose expression, STAT3 has been identified as a mediator of involution and apoptosis in the mammary gland [87], primarily activated by the cytokine leukaemia inhibitory factor (LIF) [88,89]. STAT3 is also involved in leptin signalling, interacting with the long form of the leptin receptor [90]. Mice where either STAT3 or LIF are knocked out show delayed involution and reduced levels of apoptosis [88], raising the possibility of improving efficiency in dairy herds, especially those milking once-a-day, by breeding for animals lacking one or more of these genes. It is difficult to assign QTL unambiguously to the STAT3 locus, as it maps on BTA directly between the two STAT5 genes STAT5A and STAT5B. Yet another STAT protein, STAT1, is believed to be involved in the development of the mammary gland [91]. A QTL mapping to the STAT1 locus on BTA2 has been associated with milk, fat, and protein yield [91].

One group of proteins that is upregulated by the JAK/STAT pathway is the suppressors of cytokine signalling (SOCS) gene family [92]. Proteins in this family down-regulate JAK/STAT-mediated signalling. For example, SOCS1 expression is induced by prolactin signalling, and in turn binds to JAK2, inhibiting its association with STAT5 and dampening signal transmission [68,93]. SOCS2 negatively regulates GH signalling, with knockout mice showing increased weight and skeletal dimensions [94]. Because SOCS proteins ultimately downregulate milk protein expression via prolactin and growth hormone signalling pathways, there is the potential that knocking the genes coding them out could improve lactation phenotypes in dairy animals. For example, mice with homozygous Socs1 knockout genotypes show enhanced alveolar development [93], and a point mutation in the ovine Socs2 gene has been associated with higher milk production in dairy sheep, albeit with increased susceptibility to mastitis [95]. In cattle, several different SOCS genes have been associated with milk volume, fat, and protein phenotypes [92].

2.4. Transporters and Ion Channels

Another important mechanism affecting milk production involves trans-membrane transport and ion channels. All major milk components need to either be produced within the mammary epithelial cells or transported across them from the blood, and in both cases need to cross the apical membrane into the lumen. While some small molecules like urea can simply diffuse across the membrane, in most cases either passive or active transport channels are required.

The volume of water excreted into the milk is driven by osmotic pressure, in turn created by exporting lactose and ions across the cell membrane against their concentration gradients. This means that QTL for milk phenotypes frequently map to genes encoding transporters. One important group is the sugar transporters. The gene SLC37A1 (on BTA1) encodes a glucose-6-phosphate : inorganic phosphate antiporter [96], responsible for importing glucose into the cell, and QTL for milk volume have been detected at this locus in several populations [20,85,97]. Another glucose transporter, SWEET1, is encoded by the gene SLC50A1 on BTA3. SWEET1 has been observed in the Golgi in mammary cells, and is possibly responsible for importing glucose into the Golgi for lactose synthesis [98]. QTL for lactose concentration [23] and protein concentration [85] have been identified near this gene, though the latter QTL has been assigned to GBA1, which encodes the lysosomal protein glucosylceramidase beta 1. Other glucose transporters have also been implicated in milk production, such as GLUT1 (SLC2A1 on BTA3) [99], the ${\mathrm{Na}}^{+}$/glucose co-transporter SGLT1 (SLC5A1 on BTA17) [100], and GLUT12 (SLC2A12 on BTA9) [100].

While the mammary epithelial cells are using osmotic pressure to secrete milk, it is important that they maintain their own cell volumes correctly. One mechanism by which they can do this is using voltage-regulated anion channels (VRACs), which help regulate cell volume by exporting Cl⁻ ions, as well as small organic anions such as taurine [101,102]. VRACs are comprised of heteromers of leucine-rich repeat containing 8 (LRRC8) proteins A to E, encoded by the genes LRRC8A–LRRC8E. In cattle, the locus on BTA3 containing LRRC8B–LRRC8D has been associated with milk lactose concentration [23,77], with LRRC8C suggested as the likely causative gene on the basis of gene expression data. Interestingly, LRRC8C has been associated with adipocyte differentiation (under the name fad158) [103], and is also present in the membranes of MFGs [104], suggesting it may also have a role in the storage or export of fat.

As stated above, it is important that water be able to move across the cell membrane to osmotic balance as milk is expressed. However, the lipid bilayer of the cell membrane is impermeable to water, requiring channels to facilitate the crossing. These channels are provided by aquaporin proteins (AQPs). At least seven aquaporins are expressed in the mammary gland [105,106], where they are believed to play a role in gland development, and in transporting water to the lactating gland for milk synthesis and secretion. For example, AQP1 is expressed in the capillary endothelial cells, and may be involved in oestrogen-mediated angiogenesis in the developing mammary gland [107], while AQP5 is expressed in mammary epithelial cells, and the protein pores are moved from the cytoplasm to the apical cell membrane under the regulation of prolactin [108], suggesting a role in milk production. Genetic effects mapping to aquaporin genes have been reported, such as a study in sheep [109] that mapped QTL for milk fat and protein concentrations to a window on OAR3 that contains the genes AQP2, AQP5, and AQP6, albeit alongside LALBA, which also presents a strong candidate causal gene for milk-related traits.

Another important class of transporter is the potassium channel, a type of widely expressed ion channel found in the majority of cell types. A number of different families of potassium channel have been identified. The largest family, the voltage-gated potassium channels (${\mathrm{K}}_{\mathrm{v}}$, reviewed in [110]), responds to voltage changes in the cell’s membrane potential. One channel, ${\mathrm{K}}_{\mathrm{v}}3.3$ is encoded by the gene KCNC3 on BTA18, has been associated with milk yield in cattle [22]. Another QTL for both milk volume and fat yield maps to the gene KCNS2, encoding ${\mathrm{K}}_{\mathrm{v}}9.2$, on BTA14 [83]. A third gene in the same family, KCNH4, encodes the channel ${\mathrm{K}}_{\mathrm{v}}12.3$, and is a potential candidate for a lactose concentration QTL on BTA19 [23], although the QTL also encompasses the genes STAT3, STAT5A, and STAT5B, which are also candidates. A second large family is the potassium inwardly rectifying channel family (${\mathrm{K}}_{\mathrm{ir}}$, reviewed in [111]), comprising lipid-gated channels that are activated by phosphatidylinositol 4,5-bisphosphate (${\mathrm{PIP}}_2$). In cattle, a locus on BTA19 has been associated with QTL for fat and protein concentrations [22], lactose concentration [23,77], and milk yield [7]. This locus contains two genes encoding ${\mathrm{K}}_{\mathrm{ir}}$ channels: KCNJ2 (${\mathrm{K}}_{\mathrm{ir}}2.1$) and KCNJ16 (${\mathrm{K}}_{\mathrm{ir}}5.1$), and both genes have been proposed as candidate causatives underlying the QTL. A third, smaller, family of channels is the two-pore-domain potassium channels (${\mathrm{K}}_{2\mathrm{P}}$), known as leak channels [112]. A QTL for milk fat, protein, and volume has been identified on BTA26 at the locus of the KCNK18 gene, encoding the potassium channel ${\mathrm{K}}_{2\mathrm{P}}18.1$[83]. A fourth family is the calcium- and sodium-activated potassium channels [113], of which the most well-known member is ${\mathrm{K}}_{\mathrm{Ca}}1.1$, also known as BK. The activity of ${\mathrm{K}}_{\mathrm{Ca}}1.1$ is modulated by auxiliary $\beta$ and $\gamma$ subunits [114,115], including $\gamma$1, encoded by the gene LRRC26 on BTA11, where a QTL for fat concentration has been detected [22].

The channels discussed so far are limited to moving solutes along an existing concentration gradient. Another class of transporter requires energy in the form of ATP to concentrate solutes to establish a gradient or membrane potential. Many of these belong to the ATP binding cassette family (ABC). One important example from this family is ATP binding cassette subfamily G member 2 (ABCG2). Initially identified as a xenobiotic drug transporter [116], ABCG2 also functions as a riboflavin transporter into the milk [117] and a urate transporter in the kidneys [118]. QTL for several milk phenotypes have been mapped to the ABCG2 locus on BTA6, including milk yield [7,20,83,119], fat and protein concentration [8,77,119], lactose concentration [23,77] and ${\alpha }_{\mathrm{S}1}$-CN concentration [97]. The causative variant at this QTL is generally believed to be a tyrosine to serine substitution (Y581S), identified by Cohen-Zinder et al. [119]; however the adjacent gene SPP1, encoding the protein osteopontin (involved in bone remodelling), has also been proposed as a candidate [120,121]. A second member of the ABC family is SUR2 (encoded by ABCC9 on BTA5), which forms a component of the ATP-sensitive potassium channel ${\mathrm{K}}_{\mathrm{ATP}}$, alongside SUR1 (ABCC8) and the inward rectifying channels ${\mathrm{K}}_{\mathrm{ir}}6.1$ (KCNJ8) and ${\mathrm{K}}_{\mathrm{ir}}6.2$ (KCNJ11) [122], with the exact composition varying by tissue [123]. The channel is inhibited by ATP and activated by MgADP. The ${\mathrm{K}}_{\mathrm{ATP}}$ transporter is important for glucose level sensing to control insulin release in pancreatic beta cells [124]: at low glucose levels, ATP levels are low and ADP levels raised, so the channel is open, while at high glucose levels, ATP closes the channel. This polarises the plasma membrane, a change that is detected by voltage-gated calcium channels, which then open, causing calcium to enter the cell and trigger the release of insulin. In cattle, a QTL for milk fat yield has been identified at the ABCC9 locus [22,82]. A third ATP-binding transporter is the calcium transporter SERCA2, encoded by the gene ATP2A2 on BTA17. This transporter pumps ${\mathrm{Ca}}^{2+}$ from the cytosol into the endoplasmic reticulum, whence it can be exported into milk [125]. A QTL at this locus has been detected for milk and protein yield, as well as for milk calcium [35].

Another ion transporter, showing widely reported associations with milk phenotypes, is inorganic pyrophosphate transport regulator, encoded by the gene ANKH on BTA20. The ANKH transporter controls extracellular mineralisation by regulating the levels of inorganic pyrophosphate in the extracellular matrix. In humans and mice, mutations in this transporter have been associated with arthritis and bone growth disorders linked to tissue calcification [126,127]. This transporter has previously been assumed to be a pyrophosphate transporter; however, recent work has shown that the channel in fact transports ATP, and that the production of extracellular pyrophosphate from ATP requires the enzyme ectonucleotide pyrophosphatase/phosphodiesterase 1 (ENPP1) [128]. In cattle, QTL mapping to the ANKH locus have been identified for milk yield [7,20], milk lactose concentration [23,77], and $\alpha$-lactalbumin concentration [97]. Another recent study [129] discovered that the ANK protein cycles between the plasma membrane and trans-Golgi network using clathrin-coated vesicles mediated by clathrin adaptors AP1 and AP2. The phosphatidylinositol binding clathrin assembly protein also interacts with AP2 [130], binding to the signalling molecule phosphatidylinositol and recruiting the AP2 complex to form clathrin-coated pits. This protein is encoded by the gene PICALM on BTA29, and is associated with QTL for ${\alpha }_{\mathrm{S}1}$-CN [97] and lactose concentration [23,77].

3. Identifying Candidate Causative Genes

The previous section may have given the impression that identifying the causative genes or even variants underlying QTL is easy. However, in many cases, there will be either no obvious candidate genes under the QTL, possibly because the causative variant sits in a long-range regulatory element for a distant gene. In other cases, there may be many potential candidate genes, with different selection methods highlighting different candidates. For example, the window from 42.2 to 42.5 Mbp on BTA19 envelopes several candidate genes for milk phenotypes. A 2015 study by Raven et al. [85] highlighted the genes GHDC, STAT5A, and STAT5B on the basis of differential expression and enrichment of significant variants, and proposed STAT5A as causative on the basis of knockout studies in mice. In contrast, later work by our group [23,77] used eQTL data and missense variants in strong LD with the top QTL variant to highlight the genes DHX58, GHDC, KAT2A, KCNH4, STAT5A, and STAT5B, with the two STAT5 genes again proposed as the most likely candidates. As discussed above, potassium channels such as that encoded by KCNH4 have been associated with milk at a number of loci, while the histone deacetylase and transcription activator KAT2A2 is also a possible candidate, based on the high levels of gene expression required in the lactating mammary gland. It is possible that multiple QTL are segregating at this locus, and therefore that more than one gene is causative.

Another region with two strong candidate causal genes maps between 15.3 and 15.6 Mbp on BTA3. At this locus, Raven et al. [85] identified a QTL for milk protein concentration, and proposed the epithelial mucin gene MUC1 as the best candidate. The mucin 1 protein coded by this gene is considered a “metabolic master regulator” [131], regulating tyrosine-kinase signalling and the expression of metabolic genes. Work by our group [23,77] as identified a QTL for milk lactose concentration at the same locus, and proposed the gene SLC50A1, encoding a sugar transporter. Again, it is possible that these are two separate QTL, and that both genes are causative for the corresponding phenotypes; however, another possibility is that the QTL is pleiotropic, and only one QTL is present controlling both phenotypes. Distinguishing between these two possibilities is difficult or impossible using purely statistical or bioinformatic means, and additional experiments are likely to be required.

3.1. Molecular Phenotypes

Widely used phenotypes such as fat or protein are effectively aggregating signals from a large number of milk components: individual fatty acids and other lipids for fat, and a number of casein and whey proteins for protein, for example. It is likely that these different components are not all under the same genetic regulation, and some correlations may even be negative. For example, Stoop et al [132] observed genetic correlations of as low as $-0.84$ between milk fatty acid measures (C14:0 and C16:0), as well as correlations over $0.9$. This suggests that using finer composition measures, such as individual fatty acids or proteins, should give cleaner, and possibly a larger number of, genetic signals compared to complex phenotypes such as fat or protein. This approach has been successfully applied to use milk minerals and individual protein measurements to identify novel QTL and to highlight SLC37A1 as involved in lactation [35]. These phenotypes can be considered “molecular”, as measurements for a single molecule, for example using gas or liquid chromatography with mass spectroscopy, are providing the phenotype.

Molecular phenotypes in milk are often measured using Fourier-transform mid-infrared spectroscopy, as this method is commonly used for commercial herd testing, and is cheaper than to perform on a large scale than spectroscopy. This method uses the absorbance of around 900 different frequencies of infrared light (known as wavenumbers), then uses these values in a model to predict the phenotypes of interest. In commercial herd testing, the predicted phenotypes are typically fat, protein, and lactose concentrations; however, models have been developed for a range of other phenotypes, such as individual fatty acids [133] and proteins [134], as well as phenotypes more remote to milk such as methane production [135] and fertility [136]. However, because FT-MIR phenotypes are predicted rather than measured, their utility for QTL discovery can be variable. For example, our recent work [137] showed that both HPLC measured and FT-MIR predicted phenotypes could detect significant cis-QTL for $\alpha$-casein, $\kappa$-casein, and $\beta$-lactoglobulin. However, FT-MIR was unable to identify the highly significant cis-QTL for lactoferrin that was identified using HPLC, and identified a QTL at the locus of the fat synthesis gene DGAT1 for $\alpha$-casein; this latter may be due to signal crossover from fat, which is highly correlated with protein in milk. This form of crossover, if present, would complicate the task of identifying pleiotropy.

In addition to predicting concentrations of milk components, it is also possible to use the individual wavenumber data itself as phenotypes [77]: different wavenumbers represent different chemical bonds that absorb MIR light at the corresponding frequencies, so wavenumber phenotypes are effectively measuring the concentrations in milk of compounds that contain specific chemical bonds. Wavenumber phenotypes have been shown to detect stronger and more numerous genetic signals compared to the predicted phenotypes such as fat and protein [77].

Examining the expression level of each gene in a relevant tissue sample provides an alternative to measuring or estimating protein concentrations. In milk, for example, expression levels of the genes encoding milk proteins in the lactating mammary epithelial tissue can give a proxy measurement for milk proteins. Traditionally, this could be done using techniques such as qPCR or expression microarrays; more recently, using RNA sequencing (RNA-seq) to sequence and count (relatively) mRNA molecules has become common. Because RNA-seq captures data for every expressed gene, it is possible to search for expression QTL (eQTL), as well as allele-specific expression (ASE), for every gene, including those not expressed in the milk, such as fat synthesis enzymes and hormone receptors. The presence of a QTL for a given gene can provide evidence to the causality (or otherwise) of the gene at an overlapping QTL: when the QTL is caused by an underlying eQTL, we expect that the causal variant or variants will be strongly associated with both the QTL and eQTL, with associations for other variants decaying away in proportion to linkage disequilibrium with the causal variants. This pattern can be identified by examining the correlations across variants between the QTL and eQTL variant effects [138] or p-values (on a logarithmic scale) [23,43], or correlations between local genomic estimated breeding values (GEBVs) and gene expression [139]. Other methods used to associate eQTL with GWAS results include Transcriptome-Wide Association Scan (TWAS) [140], Mendelian Randomisation [141,142] and Bayesian colocalisation methods [139,143].

Molecular phenotypes and their QTL have assisted in highlighting candidate causative genes in cases where no obvious candidates existed. For example, the gene MGST1 on BTA5 encodes microsomal glutathione S-transferase 1, belonging to a family of detoxification enzymes [144], with no obvious role in milk production. Nevertheless, QTL for milk traits mapping to this locus have been reported in many different studies [20,85,97,145,146], with the neighbouring gene EPS8 (Epidermal growth factor receptor kinase substrate 8) sometimes proposed as a candidate. However, gene expression data from lactating bovine mammary tissue have shown that the milk QTL at this locus co-segregate with an eQTL for MGST1, whereas EPS8 is barely expressed in this tissue. Similarly, a QTL for milk fat percentage on the distal end of BTA11 has been linked to the gene encoding the ABO blood group, also named ABO (alpha 1-3-N-acetylgalactosaminyltransferase and alpha 1-3-galactosyltransferase), using RNA-seq data to show that a splice donor site mutation causes aberrant splicing of the ABO transcript. in turn causing an eQTL that co-segregates with the fat percentage QTL [77].

3.2. Chromatin Structure Phenotypes

Chromatin is the name given to the compound comprising DNA wound around nucleosomes, themselves composed of histone proteins. The three-dimensional folded structure of chromatin can have an important impact on gene expression. On a proximal scale, the structure of the chromatin region surrounding a gene can be in open or closed configurations, with the former providing access for transcription factors and RNA polymerase components to reach the promoter region and trigger gene expression, while the latter, closed, configuration blocks expression. Variants sitting within an open chromatin region in the appropriate cell type are more likely to have a regulatory impact on gene expression [138], and therefore are more likely to present good candidate causal variants for a given trait, compared to other non-coding variants.

One method commonly used to study chromatin state is chromatin immunoprecipitation (ChIP) [147], which can be used to identify transcription factor binding sites or histone modifications. Histone proteins feature long tails, and a wide variety of post-translational modifications can be made to these tails to alter the chromatin state. For example, the modifications H3K4me2 (histone 3, lysine 4 dimethylation), H3K4Me3, and K3K27ac (lysine 27 acetylation) are associated with open chromatin in active promoters and enhancers, while closed, repressed regions can be marked by H3K9me2 or H3K27me3 [148]. More recent versions of ChIP, using DNA sequencing techniques to target the whole genome in a single experiment, include ChIP-seq [149,150] and CUT&RUN-seq [151]. Recent work [138] has identified cis-QTL for ChIP-seq phenotypes, and demonstrated that these QTL frequently exhibit strong correlation with nearby eQTL, and that eQTL tag variants are significantly enriched within ChIP-seq identified open chromatin windows; this illustrates the utility of of this type of data for understanding regulatory variation in the genome.

Open chromatin is accessible to transcription factors and other proteins required for gene expression. However, it is also more accessible than other regions to nuclease enzymes, and this fact is used in another approach to identifying open chromatin regions. One commonly-used nuclease is DNase 1, which can be used in conjunction with high-throughput sequencing in the technique DNase-seq [152]. The hypersensitivity sites identified using this technique highlight the positions of regulatory elements such as promoters, silencers, and enhancers. A similar technique, using the enzyme micrococcal nuclease, is called MNase-seq [153,154], and identifies the positions of nucleosomes. As nucleosomes are scarcer in open chromatin, these data give an inverse signal to hypersensitivity sites identified by DNase-seq. A more recent technique called Assay for Transposase Accessible Chromatin (ATAC-seq) [155] uses a modified hyperactive transposase enzyme Tn5 to fragment and load sequencing adaptors to open chromatin regions. ATAC-seq typically requires fewer cells and less time to perform compared to DNase-seq, although single-cell protocols have now been developed for both methods [156,157].

On a larger scale, the 3D folded structure allows distal enhancer and silencer elements to enter close contact with the gene to perform their respective functions on gene expression [158]. Identifying this structure allows linkages between genes and distal regulatory elements to be identified. These sorts of long-range interactions can be studied using chromasome conformation capture (3C) [159], and related methods such as 3C-on-chip (4C) [160] and 3C carbon copy (5C) [161]. These older techniques can study only small parts of the genome. The more recent technique Hi-C [162] could originally study the structure of the whole genome only at low resolution ($\approx 1\mathrm{Mbp}$), but later refinements of the method [163,164] have allowed for improved resolution and signal-to-noise ratios.

3.3. From Candidate Genes to Causative Variants

The results of GWAS and similar experiments are typically genomic intervals containing several associated variants that are in strong LD with one another. Although candidate causal genes may be selected using known pathways, gene expression data, or chromatin structure as discussed above, for some applications it is useful or necessary to known the causal variant underlying the QTL. These applications include improving the accuracy of genomic selection, performing genetic testing, or creating gene-edited animals [165]. However, identifying the causal variant using purely statistical or bioinformatics approaches is typically not possible, as variants in strong LD can not be distinguished from one other statistically.

In some cases, variant annotation (using tools such as SnpEff [166] or Ensembl’s Variant Effect Predictor [167]) can highlight missense or nonsense mutations that will make stronger candidates [39], but many QTL are driven by regulatory effects rather than coding ones. This means that other tools will be needed to resolve QTL with underlying regulatory effects. One commonly used technique for discovering cis-regulatory elements (CREs) is transcription factor binding site (TFBS) prediction. Several tools have been developed for this, generally using information from databases containing binding site motifs for large numbers of transcription factors, such as TRANSFAC [168] and JASPAR [169]. A second approach to investigating regulatory effects is to use a reporter assay, where the effects of putative promoter variants can be tested against the expression of a reporter gene such as GFP. This method has recently been scaled up to test thousands of variants simultaneously using a massively parallel reporter assay (MPRA) [170]. A third method is to use CRISPR-Cas9 or other gene editing technologies to test putative regulatory sequences, either by deleting them using NHEJ, or by editing in alleles of interest using HDR, then observing the resulting effect on gene expression, often using single-cell RNA-seq. Methods incorporating CRISPR include Perturb-seq [171], CROP-seq [172], and HCR-Flowfish [173].

4. Conclusions

Over the last twenty years since the seminal work of Grisart et al [40] showed that the DGAT1 gene underlies the large milk QTL on BTA14, a large number of milk related QTL have been identified in cattle, with candidate causative genes proposed for many of them. In coming years, we can expect that sequence-resolution data sets will continue both to grow, and to diversify to include additional breeds from around the world. These larger, more diverse, data sets will likely empower the discovery of many novel QTL. Although proving the causality of genes and variants underlying these QTL will likely remain difficult, the use of molecular phenotypes, massively parallel reporter assays, and CRISPR will extend the list of proven causatives. This new information should provide greater insight into the genes and pathways underlying the initiation and maintenance of lactation in mammals.