The DGAT1 gene encodes an enzyme responsible for catalysing the terminal reaction in mammary triglyceride synthesis, and underpins a well-known pleiotropic quantitative trait locus (QTL) with a large influence on milk composition phenotypes. Since first described over 15 years ago, a protein-coding variant K232A has been assumed as the causative variant underlying these effects, following in-vitro studies that demonstrated differing levels of triglyceride synthesis between the two protein isoforms. In the current study, we used a large RNAseq dataset to re-examine the underlying mechanisms of this large milk production QTL, and hereby report novel expression-based functions of the chr14 g.1802265AA>GC variant that encodes the DGAT1 K232A substitution. Using expression QTL (eQTL) mapping, we demonstrate a highly-significant mammary eQTL for DGAT1, where the K232A mutation appears as one of the top associated variants for this effect. By conducting in vitro expression and splicing experiments in bovine mammary cell culture, we further show modulation of splicing efficiency by this mutation, likely through disruption of an exon splice enhancer as a consequence of the allele encoding the 232A variant. Although the relative contributions of the enzymatic and transcription-based mechanisms now attributed to K232A remain unclear, these results suggest that transcriptional impacts contribute to the diversity of lactation effects observed at this locus.

• Genotypes
  • Animals were genotyped with the Illumina BovineHD BeadChip by GeneSeek. A subset of 115 SNPs in the 1 Mbp interval centred on DGAT1 K232A (chr14:1302265–2302265) were assessed for this study.
  • RNAseq-derived genotypes for K232A (chr14:1802265G>A SNP) were also included in the analyses, as this variant is the first base of the DGAT1 K232A MNP.
  • All genotypes were recoded using PLINK to 0, 1 or 2 to represent the number of alternative alleles for each marker (i.e. 0, 1, and 2 to represent the homozygous reference, heterozygous, and homozygous alternative genotypes, respectively).
  • In addition to these genotypes, 3,128 imputed whole-genome sequence (WGS) derived variants in the 1 Mbp interval of interest were used. These were imputed with Beagle v4, using a reference population of 556 animals (not included in this dataset).
• RNA Seq splicing phenotypes
  • RNA sequence data were mapped to the UMD3.1 genome using Tophat2, locating an average of 88.9 million read-pairs per sample.
  • The splicing efficiency phenotypes for each individual RefSeq DGAT1 junction were calculated as the ratio of exonic reads to intronic reads corresponding to the junction (of spliced and unspliced reads, respectively). Reads were considered exonic if they bridged the splicing junction i.e. mapped to the 3′ end of the preceding exon and the 5′ end of the following exon. Reads were considered intronic or unspliced if they mapped to the 3′ end of the preceding exon and through the intron-exon boundary into the intron.
• Association analyses
  • Association analyses for splicing efficiency were also conducted using ASReml-R to fit pedigree-based mixed models. This analysis fitted 115 BovineHD SNPs plus K232A genotypes derived from the RNAseq data, with phenotypes defined as described above.