A post-GWAS confirming effects of PRKG1 gene on milk fatty acids in a Chinese Holstein dairy population

Background We previously conducted a genome-wide association study (GWAS) strategy for milk fatty acids in Chinese Holstein, and identified 83 genome-wide significant single nucleotide polymorphisms (SNPs) and 314 suggestive significant SNPs. Among them, two SNPs, BTB-01077939 and BTA-11275-no-rs associated with C10:0, C12:0, and C14 index (P = 0.000014 ~ 0.000024), were within and close to (0.85 Mb) protein kinase, cGMP-dependent, type І (PRKG1) gene on BTA26, respectively. PRKG1 gene plays a key role in lipolysis to release fatty acids and glycerol through the hydrolysis of triacyglycerol in adipocytes. We herein considered it as a promising candidate for milk fatty acids. The purpose of this study was to investigate whether PRKG1 had effects on milk fatty acids. Results By direct sequencing the PCR products of pooled DNA, we identified a total of six SNPs, including one in 5′ flanking region, four in 3′ untranslated region (UTR), and one in 3′ flanking region. The single-locus association analysis was carried out, and showed that the six SNPs mainly had significant associations with C6:0, C8:0 and C17:1 (P < 0.0001 ~ 0.0035). In addition, we observed a haplotype block formed by g.6903810G > A and g.6904047G > T with Haploview 4.1, and it was strongly associated with C8:0, C10:0, C16:1, C17:1, C20:0 and C16 index (P = < 0.0001 ~ 0.0123). The SNP, g.8344262A > T, was predicted to alter the binding site (BS) of transcription factor (TF) GAGA box with Genomatix software, and the subsequent luciferase assay verified that it really changed the transcriptional activity of PRKG1 gene (P = 0.0009). Conclusion In conclusion, to our best of knowledge, we are the first who identified the significant effects of PRKG1 on milk fatty acids in dairy cattle. Electronic supplementary material The online version of this article (10.1186/s12863-019-0755-7) contains supplementary material, which is available to authorized users.


Background
Dairy products are well known for vital nutrients providing high quality protein and energy in human diet [1][2][3]. The most important economic traits of milk production in dairy cattle include milk yield, fat and protein yield, and fat and protein percentage [4]. Milk fat contains a lot of fatty acids composed of saturated fatty acid (SFA) and unsaturated fatty acid (UFA), and it determines the physiological and sensory properties of the milk [5]. SFA increases the risk of cardiovascular diseases, while UFA decreases the risk [6][7][8][9][10][11]. For example, C12:0, C14:0 and C16:0 have adverse effects on lower-density lipoprotein cholesterol, that is a risk factor for cardiovascular diseases [10], and substituting n-6 and n-3 polyunsaturated fatty acids for SFAs decreases cardiovascular diseases morbidity and mortality [11]. Fatty acids are regulated by a huge network of genes encoding transcription and translational regulators in living organisms [12], and the heritability of SFA (0.14~0.33) and UFA (0.08~0. 29) have been reported [13][14][15][16][17].
In dairy cattle, some promising candidate genes and QTL regions for milk fatty acids have been identified in previous Genome-wide association studies (GWASs), such as fatty acid synthase (FASN; on BTA19), diacylglycerol O-acyltransferase 1 (DGAT1; BTA14), stearoyl-CoA desaturase (SCD; on BTA26), 22,833,168 bp to 26,284,743 bp on BTA9 and 1,588,879 bp to 2,764, 862 bp on BTA14 [18][19][20][21][22][23]. Our previous GWAS [23] discovered 83 genome-wide significant single nucleotide polymorphisms (SNPs) and 314 suggestive significant SNPs associated with milk fatty acids in Chinese Holstein cows, in which, one SNP, BTB-01077939 for C10:0 (P = 0.000014) and C14 index (P = 0.000014), was located within the protein kinase, cGMP-dependent, type І (PRKG1) gene, and the other SNP, BTA-11275-no-rs for C12:0 (P = 0.000024), was close to the PRKG1 gene with a distance of 0.85 Mb. In addition, Li et al. reported that a QTL region (1.00 Mbp~27.94 Mbp) on BTA26 is significantly associated with milk fatty acids by performing a joint GWAS based on Chinese and Danish Holstein populations [24]. PRKG1 gene is located on BTA26 (6, 901,760~8,343,635 bp) spanning about 1442 kb and includes 20 exons. It regulates the lipolysis in adipocytes to release fatty acids and glycerol by the hydrolysis of triacyglycerol. It was reported that PRKG1 gene was involved in cGMP-PKG signaling pathway to inhibit rat brown adipocyte proliferation [25]. To date, no study has reported the associations of PRKG1 gene with milk fatty acids in dairy cattle. Hence, the objective of this study was to detect whether the PRKG1 gene had effects on milk fatty acids. We herein searched potential SNPs in PRKG1 gene, and examined associations of the identified SNPs with 24 traits in Chinese Holstein cows.
Further, we verified the impact of one regulatory SNP on transcriptional activity of PRKG1 with dual-luciferase assay.

SNPs identification
By screening the entire coding region and 5′ and 3′ flanking regions, we identified six SNPs (Table 1) in PRKG1 gene, including g.8344262A > T in 5′ flanking region, g.6904047G > T, g.6903810G > A, g.6903365C > A and g.6902878 T > G in 3′ untranslated region (UTR), and g.6901713 T > G in 3′ flanking region. The genotypic and allele frequencies of the six SNPs in PRKG1 gene were shown in Table 1.
Further, by estimating the LD among the SNPs of PRKG1 and SCD genes with a distance of 12.79 Mbp, we did not found haplotype block between the two genes ( Fig. 2), indicating that the significant effects of PRKG1 gene on milk fatty acids were not induced by LD with SCD gene.
Change of transcriptional activity caused by g.8344262A > T We predicted the change of TFBS caused by the SNP in the 5′ flanking region of PRKG1 gene using Genomatix software, and found that the allele T of g.8344262A > T created a TFBS for GAGA-Box (GAGA).
To detect whether g.8344262A > T changed the transcriptional activity of PRKG1, we synthesized two constructs with A and T of g.8344262A > T, respectively (Fig. 3a). We measured the luciferases activities of firefly and renilla, and showed the results in Fig. 3b. Luciferase activities of the two constructs were significantly higher than that of the blank control (P ≤ 0.0005) and empty vector (PGL4.14; P ≤ 0.0006), suggesting that g.8344262A > T might have the transcriptional activity. The allele T of g.8344262A > T had significantly lower luciferase activity than the allele A (P = 0.0009), implying that g.8344262A > T could alter the transcriptional activity of PRKG1 gene.

Discussion
In our initial GWAS [23], PRKG1 was considered as one of the promising candidate gene on milk fatty acids in Chinese Holstein. It was also reported that PRKG1 gene was associated with fatty acid composition in swine [26], and the PRKG1 knockout mice had lower triglyceride stores in the brown adipose tissue [27]. Here, we first verified that PRKG1 gene had significant effects on medium-chain saturated fatty acids in dairy cattle, especially C8:0.
SCD gene was also a promising candidate gene on BTA26 for milk fatty acids [23], and its effects had been confirmed [28]. In a joint GWAS based on Chinese and Danish Holstein populations, Li et al. also identified a significant QTL region for milk fatty acids (10.002 7.94 Mbp on BTA26), which included SCD and PRKG1 [24], and the SCD was in downstream of PRKG1 with a distance of 12.79 Mbp. In this study, it was shown that no LD among the SNPs of PRKG1 and SCD was observed, indicating that the effects of PRKG1 on milk fatty acid traits were independent from SCD.
From the KEGG database, we found that PRKG1 was involved in the cGMP-PKG signaling pathway (ko04022) and regulation of lipolysis in adipocytes (ko04923). In rat, cGMP signaling inhibited brown adipocyte proliferation and thereby promoted brown adipocyte differentiation [25]. The brown adipose tissue from PRKG1 knockout mice decreased triglyceride stores, suggesting an increase in the ratio of pre-adipocytes to adipocytes and fewer fully differentiated brown adipocytes [27]. In swine, the RNA-Seq analysis identified that the PRKG1 gene was the differentially-expressed in muscle between high and low groups for fatty acid composition traits [29]. Considering the significant effect of PRKG1 on milk fatty acid in the present study, we deduced that the gene might have important regulatory function for milk fatty acid metabolism in dairy cattle.
Gene expression is commonly controlled by TFs that are bound to specific sequence elements and the regular regions of the genome [30,31]. TFBSs are the biding sites (BS) targeted by a DNA binding protein [32], and its disruption is associated with phenotypic diversity [33,34]. In this study, the allele T of g.8344262A > T was predicted to create a BS for TF GAGA, and it significantly lowered the transcriptional activity of PRKG1. GAGA is a drosophila transcription factor involved in many nuclear activities, and can enhance transcription by stabilizing pre-initiation complex and promoting reinitiation [35]. Interactions of GAGA-binding proteins with the GAGA of the V1bR promoter activate V1bR gene expression, and provide a potential mechanism for physiological regulation of V1bR transcription [36]. The significant associations of g.8344262A > T with milk fatty acids, and its impact on transcriptional activity of PRKG1 gene, suggested that g.8344262A > T might be a potentially causal mutation regulating the PRKG1 expression by changing the BS for the TF GAGA to affect the formation of milk fatty acids in dairy cattle.

Conclusions
According to our previous GWAS, we considered PRKG1 gene as a promising candidate for milk fatty acids in dairy cattle. In the present study, we further confirmed the effects of PRKG1 on milk fatty acids, and showed that the gene mainly impact on medium-chain saturated fatty acid traits. In addition, we revealed that Table 3 Additive g.8344262A > T might be a potentially causal mutation altering the transcriptional activity due to the change of a BS for TF GAGA. Our findings might be helpful for the marker-assisted selection in dairy cattle.

SNP identification and genotyping
We extracted semen DNAs of 44 Holstein bulls who were sires of the above-mentioned cows using the saltout procedure, and extracted blood DNAs of the 1065 Chinese Holstein cows with the TIANamp Blood DNA Kit (Tiangen, Beijing, China) according to the manufacturer's instructions. We measured the quantity and quality of extracted DNA using a NanoDrop™ ND-2000 Spectrophotometer (Thermo Scientific, Hudson, DE, USA) and 2% agarose gel electrophoresis, respectively.
Based on the genomic sequence of the bovine PRKG1 (Gene ID: 282004), we designed 33 pairs of primers (Additional file 1: Table S1) corresponding the entire exons and 3000 bp of 5′ and 3′ flanking regions using the Primer 3.0 (http://primer3.wi.mit.edu/). The primers were synthesized in Beijing Genomics Institute (Beijing Genomic Institute, Beijing, China). We diluted the genomic DNA of each bull into the concentration of 50 ng/ μL, and randomly pooled them into two equal pools  For the SNPs identified in PRKG1 gene, we used the matrix-assisted laser desorption/ionization time of flight mass spectrometry (MALDI-TOF MS, Sequenom Mas-sARRAY, Bioyong Technologies Inc. HK) to perform the genotyping for the 1065 Chinese Holstein cows.

Prediction of the transcription factor binding site (TFBS)
We used Genomatix software suite v3.9 (http://www.genomatix.de/cgi-bin/welcome/welcome.pl?s= d1b5c9a9015b02bb3b1a806f9c03293f) to predict the change of TFBS caused by the SNP g.8344262A > T in 5′ flanking region of PRKG1 gene. In the prediction, we input 30 bp of up/down-stream sequences of g.8344262A > T, namely AGTTTAATATTTATGAAATGTCTCTCTCT-C[A]CACACACACACACACACACTCACACGCACA and AGTTTAATATTTATGAAATGTCTCTCTCTC[T]-CACACACACACACACACACTCACACGCACA, to research the different transcription factor (TF) bound for allele A and T.

Recombinant plasmid construction and luciferase assay
In this study, we used the luciferase assay to verify whether the SNP g.8344262A > T changed the transcriptional activity of PRKG1 gene. We synthesized (Genewiz, Suzhou, China) two fragments (Fig. 3a), A and T of g.8344262A > T, including NheI and HindIII restriction sites at the 5′ to 3′ termini, respectively, and cloned them into the pGL4.14 Luciferase Assay Vector (Promega, Madison, USA). We sequenced these two plasmid constructs to confirm the integrity of the insertions. Then, we purified all the plasmids using the Endo-free Plasmid DNA Mini Kit II (OMEGA, omega bio-tek, Norcross, Georgia, USA). We cultured Human Embryonic Kidney (HEK)-293 T cells in Dulbecco's modified Eagle's medium (DMEM; Gibco, Life Technologies, USA) containing 10% heatinactivated fetal bovine serum (FBS; Gibco) at 5% CO 2 and 37°C. We seeded approximately 2 × 10 5 cells per well in the 24-well plates, and then used Lipofectamine 2000 (Invitrogen, CA, USA) to transfect the cells according to the manufacturer's protocol. We transfected 500 ng constructed plasmid DNA along with 10 ng of pRL-TK renilla luciferase reporter vector (Promega) into each well. The experiments were conducted in three replicates for each construct.
We harvested the cells about 48 h after transfection, and measured the activities of firefly and renilla luciferases using a Dual-Luciferase Reporter Assay System (Promega, Madison, USA) on a Modulus microplate multimode reader (Turner Biosystems, CA, USA). We used the average statistics of three replicates as the normalized luciferase data (firefly/renilla).

Estimation of the linkage disequilibrium (LD)
We estimated the LD among the SNPs identified in PRKG1 gene in this study using Haploview 4.1 (Broad Institute of MIT and Harvard, Cambridge, MA, USA), and identified the haplotype block. In the process, 95% confidence bounds on D' were generated and each comparison was called "strong LD", "inconclusive" or "strong recombination". If 95% of informative (i.e. noninconclusive) comparisons were "strong LD", a block would be created [38]. In addition, we used the r 2 to represent the correlation coefficient between two loci.