A post-GWAS confirming the genetic effects and functional polymorphisms of AGPAT3 gene on milk fatty acids in dairy cattle

Abstract

Background

People are paying more attention to the healthy and balanced diet with the improvement of their living standards. Milk fatty acids (FAs) have been reported that they were related to some atherosclerosis and coronary heart diseases in human. In our previous genome-wide association study (GWAS) on milk FAs in dairy cattle, 83 genome-wide significant single nucleotide polymorphisms (SNPs) were detected. Among them, two SNPs, ARS-BFGL-NGS-109493 and BTA-56389-no-rs associated with C18index (P = 0.0459), were located in the upstream of 1-acylglycerol-3-phosphate O-acyltransferase 3 (AGPAT3) gene. AGPAT3 is involved in glycerol-lipid, glycerol-phospholipid metabolism and phospholipase D signaling pathways. Hence, it was inferred as a candidate gene for milk FAs. The aim of this study was to further confirm the genetic effects of the AGPAT3 gene on milk FA traits in dairy cattle.

Results

Through re-sequencing the complete coding region, and 3000 bp of 5′ and 3′ regulatory regions of the AGPAT3 gene, a total of 17 SNPs were identified, including four in 5′ regulatory region, one in 5′ untranslated region (UTR), three in introns, one in 3′ UTR, and eight in 3′ regulatory region. By the linkage disequilibrium (LD) analysis with Haploview4.1 software, two haplotype blocks were observed that were formed by four and 12 identified SNPs, respectively. Using SAS9.2, we performed single locus-based and haplotype-based association analysis on 24 milk FAs in 1065 Chinese Holstein cows, and discovered that all the SNPs and the haplotype blocks were significantly associated with C6:0, C8:0 and C10:0 (P < 0.0001–0.0384). Further, with Genomatix, we predicted that four SNPs in 5′ regulatory region (g.146702957G > A, g.146704373A > G, g.146704618A > G and g.146704699G > A) changed the transcription factor binding sites (TFBSs) for transcription factors SMARCA3, REX1, VMYB, BRACH, NKX26, ZBED4, SP1, USF1, ARNT and FOXA1. Out of them, two SNPs were validated to impact transcriptional activity by performing luciferase assay that the alleles A of both SNPs, g.146704373A > G and g.146704618A > G, increased the transcriptional activities of AGPAT3 promoter compared with alleles G (P = 0.0004).

Conclusions

In conclusion, our findings first demonstrated the significant genetic associations of the AGPAT3 gene with milk FAs in dairy cattle, and two potential causal mutations were detected.

Introduction

Milk fat is one of critical breeding objectives in dairy cattle. It is comprised of triglyceride (> 95%), diglyceride (2%), phospholipids (1%), cholesterol (0.05%) and small amount of free fatty acids (FAs) (~ 0.1%) [1]. The main components of triglyceride are glycerin and FAs, in which, the FAs act as precursors for the formation of other aroma components, such as esters and alcohols [2]. For the various milk fatty acid traits in Holstein cows, the estimated heritability values have been reported to be 0.14–0.33 for saturated fatty acids (SFAs) and 0.08–0.29 for unsaturated fatty acids (UFAs) [3,4,5,6,7].

Genome-wide association study (GWAS) is a commonly used strategy to identify potential genetic variants underlying important complex traits in human and domestic animals. So far, some candidate genes and QTL regions for milk production traits have been detected with GWA studies in dairy cattle, such as DLGAP1, AP2B1, SCD, BTA11 (1.59–3.37 Mb), and BTA3 (70.34–73.69 Mb) [8,9,10,11,12,13]. In our previous GWAS for milk FAs in Chinese Holstein cows, 83 genome-wide significant single nucleotide polymorphisms (SNPs) were detected in total [12], in which, two SNPs (ARS-BFGL-NGS-109493 and BTA-56389-no-rs) associated with C18index (P = 0.0459), were located in the upstream of 1-acylglycerol-3-phosphate O-acyltransferase 3 (AGPAT3) gene. In addition, we performed a joint GWAS for milk FAs in combined Chinese and Danish Holstein populations and found that a chromosome-wide significant QTL region of 146.29–146.31 Mb on BTA1 was associated with C18:0 [13]. The AGPAT3 gene was nearby this region with approximately 400 kb. 1-acylglycerol-sn-glycero 3-phosphate acyltransferase (AGPAT), encoded by the AGPAT3 gene, is one of the isoforms of AGPATs [14] and is involved in the glycerolipid (ko00561) and glycerophospholipid metabolisms (ko00564), and phospholipase D signaling pathway (ko04072). Mammalian AGPAT catalyzed the acylation of lysophosphatidic acid to form the phosphatidic acid that was the precursor of all glycerplipids [14]. Therefore, it was implied that the AGPAT3 gene was a promising candidate gene for milk FA traits in dairy cattle. The purpose of the present study was to further detect whether the AGPAT3 gene had significant genetic effects on milk FAs in a Chinese Holstein cow population.

Materials and methods

Animals and phenotypic data

In this study, a total of 1065 Chinese Holstein cows were used as descripted in a previous research [15], which milk samples were collected in Beijing Dairy Cattle Center (www.bdcc.com.cn) to measure milk FA contents. With the gas chromatography method, a total of 16 milk FAs (C6:0, C8:0, C10:0, C11:0, C12:0, C13:0, C14:0, C14:1, C15:0, C16:0, C16:1, C17:0, C17:1, C18:0, C18:1cis-9 and C20:0) were measured as the weight proportion of total fat weight [12]. With the phenotypes, we calculated five indices based on the formula$\frac{cis\hbox{-} 9\mathrm{unsaturated}}{cis\hbox{-} 9\mathrm{unsaturated}+\mathrm{saturated}}\times 100$ [16]: C14index = $\frac{\mathrm{C}14:1}{\mathrm{C}14:1+\mathrm{C}14:0}\times 100$, C16index = $\frac{\mathrm{C}16:1}{\mathrm{C}16:1+\mathrm{C}16:0}\times 100$, C17index = $\frac{\mathrm{C}17:1}{\mathrm{C}17:1+\mathrm{C}17:0}\times 100$, C18index = $\frac{cis\hbox{-} 9}{\mathrm{C}18:1 cis\hbox{-} 9+\mathrm{C}18:0}\times 100$, and total index = $\frac{\mathrm{C}14:1+\mathrm{C}16:1+\mathrm{C}17:1+\mathrm{C}18:1 cis\hbox{-} 9}{\mathrm{C}14:1+\mathrm{C}14:0+\mathrm{C}16:1+\mathrm{C}16:0+\mathrm{C}17:1+\mathrm{C}17:0+\mathrm{C}18:1 cis\hbox{-} 9+\mathrm{C}18:0}\times 100$. In addition, the summarized SFA and UFA, and SFA/UFA were obtained as well.

SNP identification and genotyping

Based on the genomic sequence of bovine AGPAT3 gene (Gene ID: 506607), 14 pairs of primers (Table S1) were designed by the Primer 3 version 4.0 (http://bioinfo.ut.ee/primer3-0.4.0/) and were synthesized in the Beijing Genomics Institute (Beijing, China) to amplify all the exons with partial adjacent intron region, and 3000 bp of 5′ and 3′ regulatory regions. As previously descripted [15], two DNA pools were constructed and the polymerase chain reaction (PCR) amplifications were performed with each DNA pool as template. To identify potential polymorphisms, the PCR amplification products were bi-directionally sequenced with an ABI3730XL DNA analyzer (Applied Biosystems, Foster, CA, USA). Then, the identified SNPs were genotyped for the 1065 cows by the matrix-assisted laser-desorption/ionization time of flight mass spectrometry (MALDI-TOF MS, Sequenom MassARRAY, Bioyong Technologies Inc., HK).

Linkage disequilibrium (LD) and association analyses

We estimated the LD among the identified SNPs of AGPAT3 gene with Haploview 4.1 (Broad Institute, Cambridge, MA, USA).

For association analysis, the 1065 cows were traced back to three-generation pedigrees to construct the kinship matrix (A-matrix) by SAS 9.2 (SAS institute, Cary, NC, USA), so that 3335 individuals were totally included. Single-locus and haplotype-based associations with 24 kind of milk FAs were performed by the following mixed animal model with SAS 9.2:

$$ {Y}_{ijklm}=\upmu +{G}_i+{h}_j+{l}_k+{a}_l+\mathrm{b}\times {M}_m+{e}_{ijklm} $$

Here, Y_ijklm is the phenotypic value of each milk fatty acid trait; μ is the overall mean; G_i is the fixed effect corresponding to the genotype or haplotype combination of individual i; h_j (j = 1–23) and l_k (k = 1–4) were the fixed effect of farm j and stage of lactation ll, respectively; a_l is the random polygenic effect; M_m (m = 1–293) is the fixed effect of age at calving m; b is the regression coefficient of covariate M; and e_ijklm is the random residual. Further, we calculated the additive effect (a), dominant effect (d), and allele substitution effect (α) according to $a=\frac{AA- BB}{2},d= AB-\frac{AA+ BB}{2}, and\alpha =a+d\left(q-p\right)$ [17]. Here, AA, AB and BB were the least square means of milk FAs corresponding to the genotypes, and p and q were the frequencies of allele A and B, respectively.

Prediction of changes of transcription factor binding sites (TFBSs) caused by the SNPs in 5′ regulatory region

We used the Genomatix software suite v3.9 (http://www.genomatix.de/cgi-bin/welcome/welcome.pl?s=d1b5c9a9015b02bb3b1a806f9c03293f) [18] to predict whether the four SNPs in 5′ regulatory region of AGPAT3 (g. 146702957G > A, g.146704373A > G, g.146704618A > G, and g.146704699G > A) changed the TFBSs.

Recombinant plasmid construction and luciferase assay

To detect the allele-specific effects of the SNPs g.146702957G > A, g.146704373A > G, g.146704618A > G, and g.146704699G > A on the transcriptional activity of AGPAT3 gene, eight luciferase reporter gene fragments (G and A of g. 146702957G > A; A and G of g.146704373A > G; A and G of g.146704618A > G; and G and A of g.146704699G > A) corresponding to the eight alleles of the four SNPs (Fig. 1a) were designed and synthesized (Genewiz, Suzhou, China). The eight fragments with the KpnI and Nhel restriction sites at the 5′ and 3′ termini, respectively, were cloned into the pGL4.14 luciferase assay vector (Promega, Madison, USA). In addition, all the plasmids were purified by the Endo-free Plasmid DNA Mini Kit II (OMEGA, omega bio-tek, Norcross, Georgia, USA), and were sequenced to confirm the integrity of each construct’s insertion.

Fig. 1

Luciferase assay. a Sketches of recombinant plasmids with g.146702957G > A, g.146704373A > G, g.146704618A > G, and g.146704699G > A in the 5′ flanking region of AGPAT3gene. The nucleotides in red highlight refer to the SNP. b Luciferase assay analysis of the recombinant plasmids in HEK293 cells. Blank: Blank cells. PGL4.14: Empty vector. ^**: P < 0.01

The human embryonic kidney (HEK) 293 T cells were cultured with Dulbecco′s modified Eagle′s medium (Gibco, Life Technologies) and 10% fetal bovine serum (Gibco) at 37 °C in a humidified incubator containing 5% CO₂. Before transfection, about 2 × 10⁵ cells were seeded in each 24-well plate. For eight luciferase reporter gene fragments of g. 146702957G > A, g.146704373A > G, g.146704618A > G and g.146704699G > A, 500 ng constructed plasmid was co-transfected along with 10 ng pRL-TK Renilla luciferase reporter vector (Promega) into each well. All the experiments were performed in three replicates for each construct. Approximate 48 h after transfection, the cells were harvested and the activity of both firefly and Renilla luciferases were measured with a Dual-Luciferase Reporter Assay System (Promega) on a Modulus microplate multimode reader (Turner Biosystems, CA, USA). The average statistic of three replicates were calculated as the normalized luciferase data (Firefly/Renilla).

Results

Identification of SNPs

A total of 17 SNPs of the AGPAT3 gene was detected in this study (Table 1), which consisted of four (g.146702957G > A, g.146704373A > G, g.146704618A > G and g.146704699G > A) in 5′ flanking region, one (g.146705692G > A) in 5′ untranslated region (UTR), three (g.146725085 T > C, g.146726096A > G and g.146729107A > C) in introns, one (g.146735090G > T) in 3′ UTR, and eight (g.146737188C > T, g.146737545G > A, g.146737748 T > C, g.146737849C > T, g.146737879 T > G, g.146737916 T > C, g.146737946C > T and g.146738055G > A) in 3′ flanking region. The genotype and allele frequencies of the identified SNPs were presented in Table 1.

Table 1 Information of 17 SNPs of AGPAT3 gene with genotypic and allelic frequencies

Estimation of LD among the identified SNPs of AGPAT3

We used the haploview 4.1 to estimate the LD among these 17 SNPs, and observed two haplotype blocks (Fig. 2) that was formed by four and 12 SNPs, respectively. The haplotype block 1 included four haplotype combinations, namely, H1: GAAG (38%), H2: GAAA (32.2%), H3: AGGG (26.6%), and H4: GAGG (3%), and the haplotype block 2 had six haplotype combinations: H1 = GTAAGCGTCTTC, H2 = GCACGTACTGCT, H3 = GCAATCGTCTTC, H4 = ACACGCGTCTTC, H5 = GTGATCGTCTTC, and H6 = GCAAGCGTCTTC with their frequencies of 20%, 39.8%, 13.4%, 13.4%, 7.9% and 4.1%.

Fig. 2

Linkage disequilibrium (LD) among the 17 SNPs of AGPAT3 gene. The blocks indicate haplotype blocks, and the text above the horizontal number is the SNP names. The values in boxes are pairwise SNP correlations (D´), while the bright red boxes without numbers inffer complete LD (D´ = 1). The boxes have the greater LD with the brigther red

Associations between AGPAT3 and milk FAs

The associations of the 17 SNPs with 24 milk FAs were summarized in Table 2. Among these SNPs, 17 were strongly associated with C6:0 (P < 0.0001–0.0004) and C8:0 (P < 0.0001–0.0384); 14 were significantly associated with total index (P < 0.0001–0.0318); ten were significantly associated with C10:0 (P = 0.0016–0.0151); nine were strongly associated with C17:1 (P < 0.0001–0.0149); seven were significantly associated with C20:0 (P < 0.0001–0.0072); five had significant associations with C14:0 (P < 0.0001–0.0477); five were strongly associated with C17index (P = 0.0006–0.0389); five had strong associations with C18:1cis-9 (P < 0.0001–0.0258); three had significant associations with C18:0 (P = 0.0020–0.0246); three had strong associations with SFA (P < 0.0001–0.0434); two were significantly associated with C17:0 (P = 0.0212–0.0413); two were significantly associated with UFA (P < 0.0001 and P = 0.0386); one was strongly associated with C18index (P = 0.0249); and one had significant association with SFA/UFA (P = 0.0005). However, no significant association was found with C11:0, C12:0, C13:0, C14:1, C15:0, C16:0, C16:1, C14index and C16index (P > 0.05).

Table 2 Association between 17 SNPs and milk fatty acid traits in Chinese Holstein cows (LSM ± SE)

Further, the additive (a), dominant (d), and allele substitution effects (α) of the 17 SNPs on each kind of fatty acid were calculated. Results showed that the 17 SNPs exhibited significant additive, dominant, and allele substitution effects on C6:0, C8:0, C10:0, C14:0, C16:0, C16:1, C17:0, C18:0, C18:1cis-9, C18index, C20:0, C14index, C16index, C17index, SFA, UFA, and total index (Table S2; P < 0.05). For C11:0, C12:0, C13:0, C14:1 and C15:0, none of significant additive, dominant, and allele substitution effects was found (P > 0.05).

Also, association analysis on two haplotype blocks with 24 milk FAs was performed (Table 3). The haplotype blcok1 was significantly associated with C6:0, C8:0, C10:0, C14:0, C18:0, C20:0, C17index and total index (P < 0.0001–0.0245), and the block 2 was strongly associated with C6:0, C8:0, C10:0, C14:0, C18:0, C17:1, C18:1cis-9, C18index, C20:0, C16index, C17index, SFA, UFA and total index (P < 0.0001–0.0498; Table 3). While, no significant association was found for C11:0, C12:0, C13:0, C14:1, C15:0, C16:0, C16:1, C14index and SFA/UFA (P > 0.05).

Table 3 Association between haplotype blocks in AGPAT3 and milk fatty acid traits in Chinese Holstein cows (LSM ± SE)

Prediction of TFBSs changing caused by the SNPs in 5′ regulatory region

By performing Genomatix software suite v3.9, it was predicted that that four SNPs in the 5′ regulatory region of AGPAT3 gene, g.146702957G > A, g.146704373A > G, g.146704618A > G and g.146704699G > A altered the binding sites of some transcription factors (Table 4). The allele A of g.146702957G > A created a TFBS for SMARCA3 (SWI/SNF related, matrix associated, actin dependent regulator of chromatin, subfamily a, member 3) and REX1 (REX1 transcription factor; zinc finger protein 42), respectively, and the allele G created a TFBS for VMYB (v-Myb, variant of AMV v-myb). The alleles A and G of g.146704373A > G created a TFBS for BRACH (Brachyury) and NKX26 (NK2 homeobox 6, Csx2), respectively. The allele G of g.146704618A > G created two TFBSs for ZBED4 (Zinc finger, BED-type containing 4; GC-box binding sites) and SP1 (Stimulating protein 1, ubiquitous zinc finger transcription factor). The allele G of g.146704699G > A created two TFBSs for USF1 (Upstream stimulating factor 1) and ARNT (AhR nuclear translocator homodimers), and the allele A created a TFBS for FOXA1 (Forkhead box protein A1, hepatocyte nuclear factor 3-alpha (HNF-3-alpha)).

Table 4 Changes of transcription factor binding site (TFBS) caused by the SNP in the 5′untranslated (UTR) and flanking regions of AGPAT3

Exploring for luciferase activity altered by the SNPs in 5′ regulatory region

To validate the TFBS prediction results, the luciferase assay was further performed for the four SNPs (g.146702957G > A, g.146704373A > G, g.146704618A > G and g.146704699G > A) (Fig. 1b). We observed that the luciferase activities of six constructs containing g.146704373A > G, g.146704618A > G, and g.146704699G > A, were significantly higher than that of the pGL4.14 empty vector (P < 0.0007) and blank control (P < 0.0008), while g.146702957G > A did not (P > 0.05). Further, the luciferase activity of alleles A of g.146704373A > G and g.146704618A > G were significantly higher than that of their alleles G (P = 0.0004; Fig. 1b). The luciferase activity of allele G of g.146704699G > A was higher than that of allele A, while not significant (P > 0.05). These results indicated that the transcriptional activity of the AGPAT3 gene significantly altered by g.146704373A > G and g.146704618A > G might be the reasons strongly impacted on FAs.

Discussion

This study was a follow-up investigation for our previous GWAS on milk FAs in Chinese Holstein [12]. AGPAT3 is involved in pathways related to lipid metabolism (ko00561, ko00564 and ko04072). In human, docosapentaenoic acid as the substrate of AGPAT3 protein transfers a fatty acid in sn-2 position of lysophosphatic acid, a step in the phospholipid biosynthesis pathway [19]. Here, we detected that the AGPAT3 gene mainly impacted the medium-chain milk FAs in dairy cattle.

Mammalian AGPAT catalyzed the acylation of lysophosphatidic acid to form the phosphatidic acid, which was the precursor of all glycerplipids [14]. For the AGPAT families, AGPAT1 had significant association with milk FA CLA [20], and AGPAT6 was strongly associated with C14:0, C16:0, C10:1, C12:1, C14:1 and C16:1 [21]. In our previous GWA studies [12, 13], AGPAT3 gene was identified as a candidate for two milk FAs, C18index and C18:0. In this study, using an independent Chinese Holstein population that was different from the previous GWA studies, we also observed that AGPAT3 showed a significant genetic effect on C18index and C18:0. In addition, our results revealed that the AGPAT3 had strong associations with C6:0, C8:0 and C10:0. Overall, the previous GWASs and this study suggested that AGPAT3 gene had significantly genetic effects on milk FAs.

Sequences-specific binding of transcription factors to the regulatory regions on the DNA is a key regulatory mechanism that affects gene expression and hence heritable phenotypic variation [22, 23]. Eukaryotic regulatory sequences, including enhancers and promoters, are typically between a hundred and several thousand base pairs in length, and can harbor many TFBSs [24]. It is essential to understand the evolution dynamics of transcription factor binding for understanding the evolution of gene regulation [25]. In this study, by prediction, g.146704373A > G changed the bindings of transcription factors (TFs) BRACH and NKX26, and g.146704618A > G altered the bindings of TFs ZBED4 and SP1. Further, we used the luciferase assay to verify that the alleles A of g.146704373A > G and g.146704618A > G strongly increased the transcription activity of the AGPAT3 gene than the alleles G. Previous studies showed that BRACH as a regulatory factor directly activated downstream mesoderm-specific genes to exert its mesoderm-inducing effects [26], and NKX26 restrained the transcription activity of Cx40 through the F151L missense mutation to impact the heart development [27]. ZBED4 could act as a co-repressor of nuclear hormone receptors (NHRs) by its LXXLL motifys in cones [28]. Through interfering with the recruitment of SP1 to ZNF132 promoter region, methylation of SP1-binding site can inhibit ZNF132 transcriptional expression to impact the tumor in the development of esophageal squamous cell carcinoma [29]. These reports have indicated that the TFs BRACH, NKX26, ZBED4 and SP1 could activate or repress the expression of their target genes. Based on our association analysis, the cows with the AA genotypes of g.146704373A > G and g.146704618A > G of AGPAT3, yielded significantly lower contents of C6:0 and C8:0 than those with GG genotypes. According to above, we deduced that the BRACH as a TF might activate AGPAT3 gene transcription activity by binding to the allele A of g.146704373A > G thereby reducing the contents of C6:0 and C8:0, while, the transcription factors NKX26, ZBED4 and SP1 might have the contrary effects.

Nowadays, genomic selection is the main implication for dairy cattle breeding, where the genomic chips are used. Among the SNP markers in these chips, most of them were collected from the current SNP database and almost evenly distributed across the whole genome. Hence, g.146704373A > G and g.146704618A > G of AGPAT3 as the potentially causal mutations could be put into the SNP chip instead of used in marker selection to increase selection efficiency in some specific dairy cattle populations to improve the contents of milk FAs.

Conclusion

In conclusion, through a post-GWAS approach, our study firstly indicated there were significant genetic associations between the AGPAT3 gene and milk FAs in dairy cattle. Further, we found that two SNPs in 5′ regulatory region (g.146704373A > G and g.146704618A > G) changed the transcriptional activity of AGPAT3 implying their potential causal function. These findings provided important molecular information for dairy cattle breeding.