Identification of single nucleotide polymorphisms of PIK3R1 and DUSP1 genes and their genetic associations with milk production traits in dairy cows

Abstract

Background

Previously, phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1) and dual specificity phosphatase 1 (DUSP1) were identified as promising candidate genes for milk production traits due to their being differentially expressed between the dry period and the peak of lactation in livers of dairy cows. Hence, in this study, the single nucleotide polymorphisms (SNPs) of PIK3R1 and DUSP1 genes were identified and their genetic associations with milk yield, fat yield, fat percentage, protein yield, and protein percentage, were investigated using 1067 Chinese Holstein cows from 40 sire families.

Results

By re-sequencing the entire coding region and 2000 bp of the 5′ and 3′ flanking regions of the two genes, one SNP in the 5′ untranslated region (UTR), three in the 3′ UTR, and two in the 3′ flanking region of PIK3R1 were identified, and one in the 5′ flanking region, one in the 3′ UTR, and two in the 3′ flanking region of DUSP1 were found. Subsequent single-locus association analyses showed that five SNPs in PIK3R1, rs42590258, rs210389799, rs208819656, rs41255622, rs133655926, and rs211408208, and four SNPs in DUSP1, rs207593520, rs208460068, rs209154772, and rs210000760, were significantly associated with milk, fat and protein yields in the first or second lactation (P values ≤ 0.0001 and 0.0461). In addition, by the Haploview 4.2 software, the six and four SNPs in PIK3R1 and DUSP1 respectively formed one haplotype block, and the haplotype-based association analyses showed significant associations between their haplotype combinations and the milk traits in both two lactations (P values ≤ 0.0001 and 0.0364). One SNP, rs207593520(T/G), was predicted to alter the transcription factor binding sites (TFBSs) in the 5′ flanking region of DUSP1. Further, the dual-luciferase assay showed that the transcription activity of allele T in rs207593520 was significantly higher than that of allele G, suggesting the activation of transcriptional activity of DUSP1 gene by allele T of rs207593520. Thus, the rs207593520 SNP was highlighted as a potential causal mutation that should be further verified.

Conclusions

We demonstrated novel and significant genetic effects of the PIK3R1 and DUSP1 genes on milk production traits in dairy cows, and our findings provide information for use in dairy cattle breeding.

Background

Genomic selection has been widely applied in dairy cattle breeding. The evaluation system with DNA marker technology and genomics has increased the rate of genetic progress for economic traits [1]. Zhang et al. illustrated that the available quantitative trait locus (QTL) lists detected by hundreds of genome-wide association studies (GWASs) and QTL mapping studies improved the performance of genomic prediction in dairy cattle [2]. To date, a large number of QTLs and genetic associations have been reported for milk traits (http://www.animalgenome.org/cgi-bin/QTLdb/index). Nowadays, RNA sequencing (RNA-Seq) has been proved to be an effective tool to identify crucial functional genes for complex traits in humans, domestic animals, and plants [3,4,5]. In a previous study, the liver transcriptomes of Chinese Holstein cows in the dry period, early lactation, and peak of lactation, were analyzed and the expression of phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1; P value = 0.0009) and dual specificity phosphatase 1 (DUSP1; P value = 0.00005) were found significantly decreased and increased in the peak of lactation compared to that in the dry period, respectively. Moreover, the PIK3R1 gene was linked with metabolic gene ontology (GO) terms and pathways, including protein phosphatase binding, protein transport, positive regulation of glucose import, and AMPK (adenosine 5′-monophosphate (AMP)-activated protein kinase), insulin, PI3K-Akt [phosphatidylinositol 3′-kinase (PI3K)-Akt], mTOR [mammalian (mechanistic) target of rapamycin], and Jak-STAT (janus kinase/signal transducers and activators of transcription) signaling pathways, and the DUSP1 gene was linked with inactivation of mitogen-activated protein kinase (MAPK) activity, protein binding, protein dephosphorylation, and MAPK signaling pathway [6], implying that the two genes were involved in milk metabolisms.

PIK3R1 is the regulatory subunit 1 of PI3K that plays an important role in the metabolic actions of insulin in the PI3K signaling pathway. DUSP1, which dephosphorylates c-Jun N-terminal kinase and p38 MAPK, is a negative regulator of MAPK involved with lipid, glucose and energy metabolisms, mitochondrial biogenesis, immune, and various diseases [7,8,9,10,11,12]. In addition, the PIK3R1 gene is located on chr.20:14.0466 cM within a distance of 1.95 ~ 3.97 cM to the reported QTLs that have large effects on protein yield [13, 14], and it is 0.38 ~ 2.80 Mb away from five SNPs significantly associated with milk production traits, ARS-BFGL-NGS-16696, BTA-51542-no-rs, ARS-BFGL-NGS-45127, ARS-BFGL-NGS-29910, and Hapmap50995-BTA-51556, identified by the previous GWAS [15]. The DUSP1 gene (chr.20: 5.47596 cM) was found to be within the reported QTL regions that were confirmed to have large genetic effects on fat yield [14] and protein yield [13], and within 0.044 ~ 3.50 Mb from the four SNPs, ARS-BFGL-NGS-48030, Hapmap54098-rs29010434, Hapmap49207-BTA-51446, and Hapmap36217-SCAFFOLD290026_21689, significantly associated with milk traits [15]. These data suggest that PIK3R1 and DUSP1 might be potential candidate genes for milk production traits in dairy cows. Hence, in this study, the single nucleotide polymorphisms (SNPs) of the two genes were identified by re-sequencing and then their genetic effects on milk yield, fat yield, fat percentage, protein yield, and protein percentage were investigated using single-locus and haplotype-based association analyses.

Materials and methods

Animal, sample and phenotypic data collection

Chinese Holstein cows were maintained with the same feeding conditions in 22 dairy farms belonging to the Sanyuan Lvhe Dairy Farming Centre (Beijing, China), and 1067 cows from 40 sire families were selected for the study. The semen of 40 sires and blood samples of 1067 cows respectively were collected for SNP discovery and association analysis. DNAs were extracted from the semen and blood samples using a salt-out procedure and TIANamp Blood DNA Kits (Tiangen, Beijing, China), respectively. The quantity and quality of the extracted DNA samples were measured by using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Hudson, NH, USA) and gel electrophoresis respectively. The phenotypic data for 305-day milk yield, fat yield, fat percentage, protein yield, and protein percentage were provided by the Beijing Dairy Cattle Centre (http://www.bdcc.com.cn/; Beijing, China). The descriptive statistics of the phenotypic values for milk production traits in the first and second lactations are presented in Table 1.

Table 1 Descriptive statistics of the phenotypic values for milk production traits

SNP identification and genotyping

The Primer3 (http://bioinfo.ut.ee/primer3-0.4.0/) was used with the bovine reference genome sequence (NC_037347.1) template to design primers that would amplify the entire coding region and 2000 bp of the 5′ and 3′ flanking regions of the PIK3R1 and DUSP1 genes (Additional file 1). Primers were synthesized by the Beijing Genomics Institute (BGI, Beijing, China). Forty semen DNAs with equal concentrations (2.5 ng) for each DNA, were randomly mixed into two pools of 20 sires each and used for all the polymerase chain reactions (PCR). PCR was performed with the conditions presented in Additional file 1. After the amplification, the purified PCR products were sequenced using an ABI3730XL DNA analyser (Applied Biosystems, Foster City, CA, USA), and the sequencing data were analyzed by using CHROMAS (version 2.23) to identify the potential SNPs. The identified SNPs were then individually genotyped using Sequenom MassArray for all the 1067 cows by matrix-assisted laser desorption/ionization time of flight mass spectrometry (MALDI-TOF MS, Sequenom MassARRAY, Bioyong Technologies Inc., HK).

Linkage disequilibrium (LD) estimation and association analyses

The extent of LD between the identified SNPs were estimated using Haploview 4.2 (Broad Institute of MIT and Harvard, Cambridge, MA, USA). The single-locus and haplotype-based association analyses in the first and second lactations for milk production traits were performed using the mixed procedure of SAS 9.13 software with the following animal model: y = μ + HYS + b × M + G + a + e, where, y is the phenotypic value of each trait for each cow; μ is the overall mean; HYS is the fixed effect of farm, year and season of calving; b is the regression coefficient of covariant M; M is the fixed effect of calving month; G is the genotype or haplotype combination effect; a is the individual random additive genetic effect, distributed as $\mathrm{N}\ \left(0,\mathbf{A}{\updelta}_a^2\right)$, with the additive genetic variance ${\updelta}_a^2$; and e is the random residual, distributed as $\mathrm{N}\ \left(0,\mathbf{I}{\updelta}_{\mathrm{e}}^2\right)$, with identity matrix I and residual error variance ${\updelta}_{\mathrm{e}}^2$. Each trait was analyzed separately and each SNP/haplotype block was also fitted separately. Bonferroni correction was applied, and the significant level was equal to the raw P value divided by number of genotypes or haplotype combinations. Furthermore, the additive (a), dominant (d), and substitution (α) effects were calculated using the following formulas: $\mathrm{a}=\frac{\mathrm{AA}-\mathrm{BB}}{2};\mathrm{d}=\mathrm{AB}-\frac{\mathrm{AA}+\mathrm{BB}}{2};\upalpha =\mathrm{a}+\mathrm{d}\ \left(\mathrm{q}-\mathrm{p}\right)$, where, AA, BB, and AB are the least square means of the milk production traits in the corresponding genotypes, p is the frequency of allele A, and q is the frequency of allele B [16].

Transcription factor binding site (TFBS) prediction and dual-luciferase assay

The MatInspector (https://david.ncifcrf.gov/home.jsp) was used to predict the changes of the TFBSs (Matrix similarity threshold, MST > 0.90) caused by the SNPs in the 5′ flanking or UTR region of the two genes. The fragments containing the binding sites of transcription factor in the 5′ flanking region of the DUSP1 gene were synthesized by Genewiz company (Suzhou, China) and cloned into the pGL4.14 Luciferase Assay Vector (Promega, Madison, WI). The plasmid constructs were sequenced to confirm the integrity of each insertion, and purified with an Endo-free Plasmid Maxi Kit (ComWin Biotech, Beijing, China). The human embryonic kidney 293 T (HEK 293 T) cells were grown in Dulbecco’s modified Eagle’s medium (DMEM; Gibco, Life Technologies, Carlsbad, CA) supplemented with 10% fetal bovine serum (FBS; Gibco), and maintained at 37 °C in a humidified incubator with 5% CO₂. Cells were seeded in the 24-well plates at approximately 2 × 10⁵ cells per well before transfection. For each well, 500 ng of the constructed plasmid was co-transfected along with 10 ng of pRL-TK Renilla luciferase reporter vector (Promega) using Lipofectamine 3000 (Invitrogen, CA, USA) according to the manufacturer’s protocol. All the experiments were performed in triplicate. The cells were harvested at 48 h after transfection and the activities of firefly and renilla luciferases were measured using a Dual-Luciferase Reporter Assay System (Promega) on a Modulus microplate multimode reader (Turner Biosystems, CA, USA). The normalized luciferase data (firefly/renilla) were used to calculate the average statistics of the replicates.

Results

Identification of polymorphisms

The entire coding region and 2000 bp of 5′ and 3′ flanking regions of the PIK3R1 and DUSP1 genes were re-sequenced. Six SNPs were identified in PIK3R1, including rs42590258 in the 5′ untranslated region (UTR), rs210389799, rs208819656, and rs41255622 in the 3′ UTR, and rs133655926 and rs211408208 in the 3′ flanking region. Four SNPs were found in DUSP1, including rs207593520 in the 5′ flanking region, rs208460068 in the 3′ UTR, and rs209154772 and rs210000760 in the 3′ flanking region (Table 2). The number of cows with different genotypes in the first and second lactations are presented in Table 2, and the allelic and genotypic frequencies are also shown.

Table 2 Detailed information about the identified SNPs

Single-locus association analyses with five milk production traits

For PIK3R1 gene, the single-locus association results (Table 3) showed that the rs42590258 was significantly associated with milk yield (P values = 0.0373 and 0.0233) and fat yield (P values = 0.0461 and 0.0065) in the first and second lactations, and protein percentage in the first lactation (P value = 0.0131), respectively. The SNP, rs208819656, was strongly associated with protein yield in the both lactations (P values = 0.0258 and 0.0117), fat yield in the first lactation (P value = 0.0272), and milk yield in the second lactation (P value = 0.0207), respectively. The rs41255622 and rs133655926 were significantly associated with milk, fat and protein yields in the first lactation (P values: 0.0032 ~ 0.0437). The rs211408208 had a strong association with fat yield in the first lactation (P value = 0.0111). While, the rs210389799 had no significant association with the five milk traits (P values > 0.05). Interestingly, there were no association between the six SNPs of PIK3R1 and the fat and protein percentage traits, except rs42590258 was associated with the protein percentage in the first lactation (P value = 0.0131). The allele additive, dominant, and substitution effects of the six SNPs of PIK3R1 gene were also calculated, and their significant associations with milk, fat and protein yields were found (P values < 0.05; Additional file 2).

Table 3 Associations of the SNPs in PIK3R1 and DUSP1 genes with milk production traits in two lactations in Chinese Holstein (LSM ± SE)

As for the DUSP1 gene, the genetic associations between the four identified SNPs and five milk production traits were analyzed. The results showed that the four SNPs were mainly significantly associated with milk, fat and protein yields (P values ≤ 0.0001 and 0.0418), and only one SNP, rs207593520, had strongly association with protein percentage in the first lactation (P value = 0.026). In the first lactation, rs208460068, rs209154772, and rs210000760 were significantly associated with milk and protein yields (P values: 0.0052 ~ 0.0338), while, rs207593520 had associations with fat yield (P value = 0.0201) and protein percentage (P value = 0.026). In the second lactation, the four SNPs, rs207593520, rs208460068, rs209154772, and rs210000760, were strongly associated with milk, fat and protein yields (P values ≤ 0.0001 and 0.0418; Table 3). Further, the additive, dominant and substitution effects of the four SNPs were analyzed, and their significant associations with milk, fat and protein yields were found (P values < 0.05; Additional file 2).

Haplotype-based association analyses with five milk traits

The extent of LD between the six identified SNPs in PIK3R1 were estimated using Haploview 4.2, and six SNPs were found highly linked (D′ > 0.89; Fig. 1) in one block. The haplotype block was formed by five haplotypes, H1 (GCTTTA), H2 (GTATTG), H3 (ATACCG), H4 (GCTTTG), and H5 (GTATTA), with the frequency of 34.7%, 18.7%, 26.8%, 17.6%, and 1.3%, respectively. The haplotype-based association analysis showed that the haplotype block was significantly associated with milk yield (P value < 0.0001), fat yield (P value = 0.0002), protein yield (P value < 0.0001), and protein percentage (P value = 0.0002) in the first lactation, and milk yield (P value = 0.0006), fat percentage (P value = 0.0305), protein yield (P value = 0.003), and protein percentage (P value = 0.0142) in the second lactation, respectively (Table 4).

Fig. 1

Linkage disequilibrium estimated among the SNPs in DUSP1 (a; D′ =0.99 ~ 1.00) and PIK3R1 (b; D′ =0.89 ~ 1.00). The blocks indicate haplotype blocks and the text above the horizontal numbers is the SNP names. The values in boxes are pairwise SNP correlations (D′), while bright red boxes without numbers indicate complete LD (D′ = 1)

Table 4 Associations of haplotype blocks with milk production traits in two lactations in Chinese Holstein (LSM ± SE)

As for DUSP1, the four SNPs were found highly linked in one haplotype block with D´ > 0.99 (Fig. 1). The haplotype block was consist of three haplotypes, H1 (CCGG), H2 (CCGT), and H3 (GTAT), with the frequency of 35.5%, 32.9%, and 30.6%, respectively. The haplotype block was significantly associated with milk yield (P value = 0.0364), protein yield (P value = 0.0054) and percentage (P value = 0.0131) in the first lactation, and milk, fat and protein yields (P values < 0.0001) in the second lactation, respectively. While, no genetic association was observed between the haplotype block and fat percentage in the both two lactations (P values > 0.1; Table 4).

Transcriptional activity of DUSP1 increased by allele T of rs207593520

The MatInspector was used to predict the changes of TFBSs due to the SNPs in the 5′ flanking or UTR region of the two genes, and the allele T of rs207593520 in DUSP1 was found to create the binding sites for the transcription factor MYB proto-oncogene like 1 (MYBL1; MST = 0.91), and the allele G invented the binding sites for kruppel like factor 12 (KLF12; MST = 0.96; Fig. 2a). Further, two plasmids containing allele T or G in the rs207593520 were synthesized for dual-luciferase assay to observe the changes of the transcriptional activity of DUSP1. As the results shown in Fig. 2b, the luciferase activities of two constructs were significantly higher than that of the empty vector (PGL4.14) and blank control (P value < 0.0001), confirming the regulatory role of rs207593520. Moreover, the relative luciferase activity of allele T was significantly higher than that of allele G (P value < 0.0001), implying that the allele T of rs207593520 in DUSP1 might have higher transcriptional activity than the allele G (Fig. 2b).

Fig. 2

Dual-luciferase activity assay. a Sketches of recombinant plasmids with rs207593520 (T/G) in the 5′ flanking region of DUSP1 gene. Underlined nucleotides represent the transcription factor binding site sequences of MYBL1 or KLF12 (in blue), and the nucleotide in red was the SNP. b Luciferase activity analysis of the recombinant plasmids in HEK 293 T cells. PGL4.14 was the empty vector, Blank was the blank cell, and *** P value < 0.0001

Discussion

Our previous RNA-Seq work considered PIK3R1 and DUSP1 as the candidate genes for milk production traits, and this follow-up investigation first demonstrated that the polymorphisms of the two genes were both mainly significantly associated with milk, fat and protein yields. To our knowledge, the PIK3R1 occupies a center role in the insulin signaling pathway, and it was reported to be associated with alterations in glucose and insulin homeostasis [17]. It can recruit protein kinase A to the lipid droplet in conveying endogenous glucocorticoid-induced lipolysis [18]. Studies also revealed that PIK3R1 influences the serum leptin and body fat [19], apolipoprotein B, and low density lipoprotein cholesterol [20] in female. DUSP1 is belongs to the MKP phosphatase family, and it regulates the MAPK signaling pathway by the dephosphorylation. Studies have uncovered an important regulatory role for DUSP1 in hepatic lipid metabolism [7, 8, 10]. Mice lacking DUSP1 were resistant to the acquisition of a fatty liver suggesting that DUSP1 negatively regulates hepatic fatty acid oxidation [10]. Lawan et al. [9, 21] have reported the contribution of DUSP1 for glucose and energy metabolisms. These data collectively illustrate that PIK3R1 and DUSP1 participate in substance metabolisms, especially lipid metabolisms, and their polymorphisms were found to be significantly associated with milk, fat and protein yields by the association analyses in this study. Additionally, the haplotype analyses were generally used to the genetic variation studies [22, 23]. Our haplotype-based association analyses showed that the six and four SNPs of PIK3R1 and DUSP1 were respectively highly linked, and the haplotype block of the two genes were all significantly associated with milk traits in Holstein cows, which were consistent with the genetic associations of SNPs with the milk traits.

In the present study, the allele T of rs207593520 in the 5′ flanking region of DUSP1 was predicted to invent the TFBSs for MYBL1 and the allele C for KLF12, and further, the transcriptional activity of DUSP1 was found significantly increased by the allele T of rs207593520 in the dual-luciferase assay. It is generally known that the SNPs in TFBSs could lead to allele-specific binding of transcription factors thereby activating or suppressing the gene expression [24,25,26]. MYB proteins are nuclear DNA-binding proteins that act as transcriptional transactivators of many genes [27], and MYBL1 has been reported as a master regulator of meiotic genes that are involved in multiple meiotic processes [28]. The transcription factor MYBL1 activates the murine tissue-specific lactate dehydrogenase expression by binding the cAMP-responsive element site [29]. In addition, MYBL1 can act as a transcriptional repressor to suppress the expression of multiple anthocyanin pigment pathway genes [30]. Transcription factor KLF12 is a member of KLFs family, which regulates gene transcription through binding to the CACCC sequence of target genes [31, 32]. Studies also showed that it can bind to the promoter regions of target genes and represses their expression [32,33,34,35]. KLF12 acts to negatively regulate the expression of the decidual marker genes decidual prolactin and insulin like growth factor binding protein 1 in human endometrial stromal cells [32]. KLF12 binds to the promoter region of leukemia inhibitory factor and directly represses its transcription [35]. While, KLF12 can directly activate the expression of early growth response protein 1 to promote the colorectal cancer growth [36]. These data suggest that the transcription factors MYBL1 and KLF12 can activate or repress the expression of their target genes, however, MYBL1 mainly acts as an activator, and KLF12 prefers to be a repressor. Based on our results, we speculated that MYBL1 might activate the expression of DUSP1 by binding the TFBSs caused by allele T of rs207593520, thereby regulating the milk production, and KLF12 might inhibit the DUSP1 expression through binding the TFBSs caused by the allele G to affect the milk production. Hence, the SNP, rs207593520, may be a potential causal mutation for milk yield traits in dairy cows because of its ability to change the transcriptional activity of DUSP1, and the further functional studies are required to validate its role.

Genomic selection is widely used in dairy cattle breeding, and the development of efficient SNP markers can improve the accuracy of the selection. Studies have shown that the SNPs in the functional genes significantly influenced the milk production traits in dairy cattle [37,38,39,40]. Thus, the SNPs with large genetic effects on milk traits could be used as markers to increase the selection efficiency in specific dairy cattle populations. Certainly, the significant SNPs of PIK3R1 and DUSP1 genes identified in this study could also be used for genomic selection in dairy cattle.

Conclusion

In conclusion, this is the first study to reveal the significant genetic effects of PIK3R1 and DUSP1 genes on milk production traits in dairy cows, and the valuable SNPs could be used for the genomic selection in dairy cattle. In addition, the rs207593520 in DUSP1 was highlighted as a functional mutation for milk traits that could change the transcriptional activity of DUSP1 gene. Further, the functional validation experiments should be performed to reveal the molecular regulatory mechanisms of PIK3R1 and DUSP1 genes on milk synthesis metabolism.