Skip to main content

Identification of genetic associations and functional SNPs of bovine KLF6 gene on milk production traits in Chinese holstein



Our previous research identified the Kruppel like factor 6 (KLF6) gene as a prospective candidate for milk production traits in dairy cattle. The expression of KLF6 in the livers of Holstein cows during the peak of lactation was significantly higher than that during the dry and early lactation periods. Notably, it plays an essential role in activating peroxisome proliferator-activated receptor α (PPARα) signaling pathways. The primary aim of this study was to further substantiate whether the KLF6 gene has significant genetic effects on milk traits in dairy cattle.


Through direct sequencing of PCR products with pooled DNA, we totally identified 12 single nucleotide polymorphisms (SNPs) within the KLF6 gene. The set of SNPs encompasses 7 located in 5′ flanking region, 2 located in exon 2 and 3 located in 3′ untranslated region (UTR). Of these, the g.44601035G > A is a missense mutation that resulting in the replacement of arginine (CGG) with glutamine (CAG), consequently leading to alterations in the secondary structure of the KLF6 protein, as predicted by SOPMA. The remaining 7 regulatory SNPs significantly impacted the transcriptional activity of KLF6 following mutation (P < 0.005), manifesting as changes in transcription factor binding sites. Additionally, 4 SNPs located in both the UTR and exons were predicted to influence the secondary structure of KLF6 mRNA using the RNAfold web server. Furthermore, we performed the genotype-phenotype association analysis using SAS 9.2 which found all the 12 SNPs were significantly correlated to milk yield, fat yield, fat percentage, protein yield and protein percentage within both the first and second lactations (P < 0.0001 ~ 0.0441). Also, with Haploview 4.2 software, we found the 12 SNPs linked closely and formed a haplotype block, which was strongly associated with five milk traits (P < 0.0001 ~ 0.0203).


In summary, our study represented the KLF6 gene has significant impacts on milk yield and composition traits in dairy cattle. Among the identified SNPs, 7 were implicated in modulating milk traits by impacting transcriptional activity, 4 by altering mRNA secondary structure, and 1 by affecting the protein secondary structure of KLF6. These findings provided valuable molecular insights for genomic selection program of dairy cattle.

Peer Review reports


In recent years, researchers have implemented multiple strategies, such as genome-wide association study (GWAS) [1,2,3,4,5], RNA sequencing (RNA-seq) [6,7,8,9], whole genome bisulfite sequencing (WGBS) [10,11,12,13,14], signature selection [15,16,17,18], to explore the genetic mechanisms underlying complex traits in human, mouse and domestic animals, and numerous functional genes and genetic associations have been identified. In our previous RNA-seq study, we found a significant upregulation of the Kruppel-like factor 6 (KLF6) gene in the livers of Holstein cows during the peak lactation period in comparison to both the dry and early lactation periods [19]. The KLF6 protein participates in lipid metabolic processes by activating peroxisome proliferator-activated receptor α (PPARα) signaling pathway. Notably, KLF6 was located near the peak regions of the known quantitative trait loci (QTL) for milk fat yield (0.09 Mb), milk fat percentage (0.28 Mb), milk protein yield (0.45 Mb) and protein percentage (1.38 Mb) ( Therefore, we considered the KLF6 gene as a promising candidate for milk production traits.

The KLF6 gene is located on bovine chromosome 13 and consists of 2 exons, with a total length of 8460 base pairs (bp). The KLF6 protein belongs to the specificity proteins/Krüppel-like factor (SP/KLF) transcription factors family, known for its regulation of various genes involved in fundamental cellular and biological processes. These processes encompass vascular remodeling [20], cell differentiation, proliferation, cycle and apoptosis [21,22,23], tissue development and growth [24,25,26,27] in humans.

Overexpression of KLF6 has been demonstrated to increase the number of lipid droplets and elevate the expression of peroxisome proliferator-activated receptor-γ (PPARγ) and CCAAT/enhancer binding protein α (C/EBPα) in pig adipose-derived stem cells. Conversely, treatment with small interfering RNAs (siRNAs) targeting KLF6 resulted in the opposite consequence [28]. The knocked-down of KLF6 negatively influences both adipogenesis and differentiation in primary adipocyte cells [29]. Moreover, polymorphisms associated with the KLF6 gene have shown significant correlations with intramuscular fat traits in Qinchuan cattle [30]. Nevertheless, no prior investigations have explored the influence of the KLF6 gene on milk traits in dairy cattle. Thus, the primary objective of the present study was to ascertain the genetic effects of the KLF6 gene on milk yield and composition traits and to identify potential functional genetic variations in Chinese Holstein cattle population.


SNPs identification

We extracted genomic DNA from 1123 blood samples of Chinese Holstein cows (Fig. 1). By sequencing the PCR products of the KLF6 gene for the entire coding sequence and 2000 bp of the 5’ and 3’ flanking regions, we identified a total of 12 SNPs, including seven in 5′ flanking region, two in exon 2 and three in 3′ UTR. Among these, the g.44601035G > A was a missense mutation that changed from arginine to glutamine of the amino acid. The genotype and allele frequencies of the identified SNPs were presented in Table 1.

Fig. 1
figure 1

Agarose gel electrophoresis result of 10 DNA samples extracted from the blood samples of cows. Lane 1 is λDNA /HindIII Marker

Table 1 Detailed information of 12 SNPs identified in the KLF6 gene

Associations between SNPs and the five milk production traits

Using SAS 9.2, we assessed the genetic associations between the 12 identified SNPs and five milk traits in Holstein cows, including milk yield, fat yield, protein yield, fat percentage and protein percentage (Table 2).

Table 2 Associations of the SNPs in KLF6 gene with milk production traits in two lactations in Chinese Holstein (LSM ± SE).

In the first lactation, we observed all the 12 SNPs in KLF6 showed significant associations with milk protein yield (P < 0.0001 ~ 0.0084), ten with milk fat yield (P < 0.0001 ~ 0.0383), nine with milk protein percent (P < 0.0001 ~ 0.0441), seven with milk fat percent (P < 0.0001 ~ 0.0335) and six with milk yield (P < 0.0001 ~ 0.0221). In the second lactation, all the 12 SNPs in KLF6 were significantly associated with milk yield (P < 0.0001 ~ 0.0183) and protein yield (P < 0.0001 ~ 0.0221), ten with fat yield (P < 0.0001 ~ 0.096), nine with fat percent (P < 0.0001 ~ 0.0003) and seven with protein percent (P < 0.0001 ~ 0.0021).

In addition, as the results showed in Table S1, the additive, dominant, and substitution effects of the 12 SNPs on milk traits were significant as well (P < 0.05). The phenotypic variance ratios explained by these SNPs of five milk traits ranged from less than 0.0001–2.1301%, and the principal explanations were milk yield and protein yield (Table S2).

Associations between haplotype and the five milk traits

We estimated the degree of linkage disequilibrium (LD) among the 12 SNPs in KLF6 and inferred one haplotype block using Haploview4.2 (D′ = 0.829-1.0; Fig. 2). Regarding the block, the frequencies of the haplotype H1 (AGGCACACGCGG), H2 (GAACATGTATGA), H3 (GAATCCGCGCGA), H4 (GAATCCGCGCTA), H5 (GAACATGCGTGA) and H6 (GAACACGCGTGA) were 38.5%, 23%, 13.5%, 12.1%, 7.9% and 3.1%, respectively.

Fig. 2
figure 2

Linkage disequilibrium estimated among the 12 SNPs in KLF6 (D′=0.829 ~ 1.0). The text above the horizontal numbers is the SNP names. The values within boxes are pairwise SNP correlations (D′), while bright red boxes without numbers indicate complete LD

Similarly, the haplotype-based association analysis revealed that the haplotypes of KLF6 were significantly associated with milk yield (P = < 0.0001), fat yield (P = 0.0203), protein yield (P = < 0.0001) and protein percent (P = 0.003) in the first lactation, and all the five milk traits in the second lactation (P = < 0.0001 ~ 0.009) (Table 3).

Table 3 Associations of haplotypes with milk production traits in two lactations in Chinese Holstein (LSM ± SE).

Functional variation prediction and verification caused by SNPs

We used Jaspar software to predict the transcription factor binding sites (TFBSs) change of the seven SNPs (g.44,595,327 A > G, g.44595580G > A, g.44595695G > A, g.44,595,808 C > T, g.44,596,221 A > C, g.44,596,271 C > T and g.44,596,874 A > G) in the 5′ flanking region of KLF6, and found that all the SNPs changed TFBSs (Table 4).

Table 4 Changes of TFBSs caused by the SNPs in the 5′flanking region of KLF6.

Further, we performed the luciferase assay to validate the above prediction results (Fig. 3). The luciferase activities of the construct G of g.44,595,327 A > G, A of g.44595580G > A, A of g.44595695G > A, T of g.44,595,808 C > T, C of g.44,596,221 A > C, T of g.44,596,271 C > T, G of g.44,596,874 A > G were observed significantly higher than those of the blank control (P < 0.0001), empty vector pGL4.14 (P < 0.0001), and the construct A of g.44,595,327 A > G, G of g.44595580G > A, G of g.44595695G > A, C of g.44,595,808 C > T, A of g.44,596,221 A > C, C of g.44,596,271 C > T, A of g.44,596,874 A > G, respectively (P < 0.0001 ~ 0.0015). These results indicated the seven SNPs in the 5′ flanking region regulated the transcriptional activity of KLF6 gene.

Fig. 3
figure 3

Luciferase assay result of the recombinant plasmids in HEK293 cells. Blank: blank cells. pGL4.14 + pRL-TK: empty vector. The nucleotides in red highlight referred to the mutation compared to the first plasmids. **: P < 0.01

We utilized the RNAfold web server to predict the secondary structures of mRNA for the 5 SNPs in UTR and exon regions of KLF6 gene (Table 5). When T replaced C in g.44,600,440 and A replaced G in g.44,603,687, the minimum free energy (MFE) of mRNA secondary structures decreased, and KLF6 got more stable. Conversely, when T replaced C in g.44,601,887 and G in g.44,602,809, the MFE of mRNA secondary structures increased, and KLF6 got more unstable. These observations suggested that the 4 SNPs might change the KLF6 mRNA secondary structure to affect the KLF6 expression.

Table 5 The minimum free energy (MFE) values of optimal secondary structure of KLF6 mRNA.

Furthermore, we employed the SOPMA software to predict the protein secondary structure for the missense mutation (g.44601035G > A) in exon 2 of the KLF6 gene. The analysis revealed that the beta turn was changed from 6.47 to 5.83% and the random coil from 60.84 to 61.49%. However, according to the prediction result from PROVEAN, there was no observed change in the protein function (score = -0.038).


Our previous liver transcriptome study of Holstein cows identified KLF6 gene as a promising candidate affecting milk production traits. Not only did it exhibit differential expression across various lactation periods, but it was also found to be intricately linked to the PPARα signaling pathway and located near the known QTLs for milk yield and composition traits. In this study, through single locus and haplotype-based association analysis, we confirmed the significant genetic effects of KLF6 gene on milk yield, milk fat and protein traits in dairy cattle.

The KLF6 gene was highly expressed in various tissues and played a critical role in adipocyte differentiation by activating PPARα that regulated hepatic steatosis, lipoprotein synthesis, hepatic gluconeogenesis and fatty acid transport proteins in mouse [31,32,33]. KLF6 also promoted preadipocyte differentiation by suppressing a factor named delta-like1 (DLK1), which acted by maintaining the preadipocyte state and preventing adipocyte differentiation in mouse preadipocytes [34, 35]. These studies demonstrated the involvement of KLF6 gene in lipid metabolism, and were consistent with the results in this study that KLF6 gene has significant effects on milk fat traits. In addition, we found KLF6 gene was significantly associated with milk protein traits as well. This is probably due to the high genetic correlation between milk fat and protein traits in dairy cattle [36,37,38].

Transcription factors, which belong to the family of sequence-specific DNA-binding proteins, regulate gene expression across diverse organisms [39, 40]. Specifically, SNPs located within transcription factor binding sites (TFBSs) can contribute to the allele-specific binding of transcription factors. In the present study, we identified 7 SNPs within 5’ regulatory region of the KLF6 gene that altered TFBSs, subsequently activating the KLF expression. These TFBS modifications corresponded to transcription factors GATA2, GATA1, GATA3, NR1I3, FOXL1, SOX15, ATF1, SOX10 and KLF12 which have been shown to promote the expression of target genes in human and mouse in several prior studies [41,42,43,44,45,46,47,48]. Conversely, transcription factors SOX18, EVX2, EMX2, EVX1, FOXC1, NPAS4, KLF4, HIC2 and ERG have been reported to inhibit the expression of their target genes [49,50,51,52,53,54,55,56]. In this study, for the g.44,595,808 C > T and g.44,596,221 A > C, the positive genetic effects of alleles T and C on milk protein may be due to the activated KLF6 expression induced by transcription factors NR1I3 and ATF1, respectively. The regulatory roles of these SNPs on KLF6 expression thereby impacting milk traits need more in-depth investigations to be validate.

Furthermore, in our study, the missense mutation in exon 2, g.44601035G > A, was predicted to impact the secondary structure of the KLF6 protein. As the backbone of advanced RNA function, the secondary structure of mRNA plays a pivotal role in various biological processes, including protein folding and transport, initiation and extension of translation process, regulation of translation rate, and directly affects the stability of mRNA itself [57, 58]. In addition, the secondary structures of KLF6 mRNA corresponding to allele T of g.44,600,440 C > T, allele C of g.44,601,887 C > T, allele G of g.44602809G > T and allele A of g.44603687G > A exhibited increased stability in comparison to the alleles C, T, T and G, respectively. The enhanced stability of the mRNA secondary structure associated with allele C of g.44,601,887 C > T may provide insight into its active genetic effects on milk yield, fat yield, and protein traits.

Studies have shown that the incorporation of functional gene information associated with substantial genetic effects on target breeding traits can improve the accuracy of genomic evaluation [59, 60]. These significant SNPs identified in this study could be used as molecular markers for genetic improvement programs of dairy cattle through genomic selection.


Based on our preceding RNA-seq investigation that identified KLF6 genes as a promising candidate for milk traits in dairy cattle, this study first demonstrated the significant genetic effects of KLF6 gene on milk yield and composition traits. A total of 12 SNPs were identified, 7 of them altered the transcriptional activity of KLF6, 4 changed the KLF6 mRNA secondary structures, and 1 missense mutation changed KLF6 protein secondary structure. These SNPs could be used as valuable genetic markers for molecular breeding program by incorporating them into genomic selection chips in dairy cattle.

Materials and methods

Animals and phenotypic data

We gathered the blood samples of 1123 Chinese Holstein cows from 80 sire families from Qingyuanlvze and Shuangfeng dairy farms (Baoding, Hebei, China) in Hebei province. The pedigree and phenotypic data of five milk traits, containing milk yield, milk fat yield, milk fat percentage, milk protein yield, and milk protein percentage were provided by the Hebei Province Animal Husbandry and Fine Breeds Work Station (Shijiazhuang, Hebei, China). Descriptive statistics of the phenotypic values for milk traits in the first and second lactations were presented in Table 6.

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

SNP identification and genotyping

We extracted genomic DNA from whole blood samples of 1123 cows using the TIANamp Blood DNA Kit (Tiangen, Beijing, China) and measured the quantity and quality of DNA samples with NanoDrop 2000 spectrophotometer (Thermo Scientific, Hudson, NH, USA) and agarose gel electrophoresis (1.5%), respectively. Eighteen pairs of primers were designed with Primer3web version 4.1.0 ( and synthesized by the Beijing Genomics Institute (Beijing, China) to amplify the entire coding sequence and 2000 bp of the 5’ and 3’ flanking regions for the bovine KLF6 gene in accordance with its genomic sequence (GenBank IDs: NC 037340.1) (Table S3).

We randomly selected 111 blood DNA samples from the 80 sire families to construct five DNA pools at the same concentration of 50 ng/µL, including one DNA pool with 23 samples and four pools with 22 samples in each. Based on the DNA pools as templates for PCR amplification (Table S4), the PCR products were bi-directionally sequenced on ABI3730XL DNA analyser (Applied Biosystems, Foster City, CA, USA) and aligned to the bovine reference sequences (ARS-UCD1.2) with BLAST ( to identify potential SNPs. The matrix-assisted laser desorption/ionization-time of flight mass spectrometry assay (MALDI-TOF MS, Sequenom MassARRAY, Agena, San Diego, USA) was performed for subsequent individual genotyping of the identified SNPs for the 1123 Holstein cows.

Linkage disequilibrium (LD) estimation

We calculated the extent of LD among the 12 identified SNPs in KLF6 with Haploview 4.2 (Broad Institute, Cambridge, MA, USA). The D’ or r2 value was used to measure the degree of LD. The Haplotypes with frequencies of less than 0.05 were discarded.

Association analyses between SNPs/ haplotypes on milk production traits

The pedigrees of 1123 individuals were traced back three generations for the association studies between the detected SNPs or haplotypes and five milk traits. SAS 9.2 mixed procedure was performed with the following model: Y = µ + hys + b × M + G + a + e, where Y is the phenotypic value of five milk traits of each cow; µ is the overall mean; hys is the fixed effect of farm (1–2 for the two farms, respectively), calving year (1–7 for the years of 2013–2019, respectively), and calving season (1 for April to May; 2 for June to August; 3 for September to November, and 4 for December to March); b is the regression coefficient of covariant M; M is calving age as a covariant; G is the genotype or haplotype combination effect; a is the individual random additive genetic effect, distributed as \({\text{N}}(0,{\text{A}}\delta _{\text{a}}^2)\), with the additive genetic variance \({{\delta }}_{\text{a}}^{2}\); and e is the random residual effect, distributed as \(\text{N} \left(0, \text{I}{{\delta }}_{\text{e}}^{2}\right)\), with identity matrix I and residual error variance \({{\delta }}_{\text{e}}^{2}\). Multiple tests were implemented by Bonferroni correction, with the significance level equal to the original P value divided by the number of genotype or haplotype combinations.

In addition, the additive (a), dominant (d), and substitution (α) effects were also computed as follows: \(\text{a}=\frac{\text{A}\text{A}-\text{B}\text{B}}{\begin{array}{c}2\\ \end{array}}\),\(\text{d}=\text{A}\text{B}-\frac{\text{A}\text{A}+\text{B}\text{B}}{\begin{array}{c}2\\ \end{array}}\),\({\alpha }=\text{a}+\text{d}(\text{q}-\text{p})\), where AA, AB, and BB were represented the least square means of milk traits attributes matching to genotypes, while p was the allele frequency of A and q was the allele frequency of B.

The phenotypic variance ratio of five milk traits explained by SNPs was calculated as follow: phenotypic variance ratio = \(\,{{2pq{\alpha ^2}} \mathord{\left/{\vphantom {{2pq{\alpha ^2}} {\sigma _p^2}}} \right.\kern-\nulldelimiterspace} {\sigma _p^2}}\), where p and q were the allele frequencies of A and B, α was substitution effects and \({\sigma }_{p}^{2}\) was the phenotypic variance of the target trait.

Biological function prediction

Using the Jaspar software ( to determine whether SNPs in 5′ flanking region of the KLF6 gene altered transcription factor binding sites (TFBSs) (relative score ≥ 0.85).

To predict changes in mRNA secondary structures for SNPs in UTR and exon regions, RNAfold Web Server ( was used. The minimum free energy (MFE) of the optimal secondary structure reflects the stability of mRNA structure. A lower MFE value indicates greater stability in the mRNA structure.

To predict the changes in protein secondary structure caused by missense mutation in the coding regions of genes, we applied SOPMA ( sopma.html) with the following parameters: similarity threshold: 8; number of states: 4-Helix, sheet, Turn, Coil; output width: 70.

Subsequently, we implemented PROVEAN software ( to confirm whether the missense mutation changed protein function. The threshold value for PROVEAN is set at -2.5, below which alterations in protein function are considered.

Construction of recombinant plasmid, cell culture and luciferase assay

To examine the effects of the SNPs that were predicted to alter the transcription factor binding sites in the 5′ flanking region of KLF6, a total of 8 luciferase reporter gene fragments (Fig. 4) were designed and synthesized in Hitrobio (Beijing, China). These fragments contained 2000 bp of the 5′ flanking sequences of KLF6, with XhoI and HindIII restriction sites at the 5′ and 3′ termini, and were cloned into the pGL4.14 Luciferase Assay Vector (Promega, Madison, USA). These plasmids were subsequently purified using the EndoFree Mini Plasmid Kit II (Tiangen, Beijing, China), and their integrity of each construct’s insertions were sequenced to verify.

Human Embryonic Kidney (HEK)-293 T cells (Procell, Wuhan, China) were cultured at 37 °C and 5% CO2 with Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 10% heat-inactivated fetal bovine serum (Biosun, Shanghai, China). Approximately 1 × 105 cells were seeded into 48-well plates and co-transfected using Lipofectamine 3000 (Invitrogen, CA, USA) with 200 ng of the constructed plasmid DNA and 20 ng of pRL-TK renilla luciferase reporter vector (Promega, WI, USA) per well. Three replicates were conducted for each construct. After 36 h, the 293T cells were collected and then the activities of firefly and renilla luciferases were measured using the DualLuciferase Reporter Assay System (Promega, WI, USA) with a multifunctional microplate detection system (BioTek, NY, USA). The average of three replicates was then calculated to obtain the normalized luciferase data (firefly/renilla).

Fig. 4
figure 4

Sketches of recombinant plasmids. The nucleotides in red highlight referred to the mutation compared to the first plasmids

Data Availability

The information of genetic polymorphisms loci can be inquired on European Variation Archive (EVA; using the RS IDs (rs41692335, rs110464810, rs110289079, rs208019372, rs211266130, rs29024529, rs41692337, rs209109676, rs211273884, rs41692344, rs208700974 and rs41692345) in Table 2. The datasets generated during and analyzed during the current study are available in the article and its additional files.



Additive effect


Coding sequence


CCAAT enhancer-binding proteins alpha


CCAAT enhancer-binding proteins beta/gamma


Dominant effect

DLK1 :

Delta-like 1

KLF6 :

Kruppel Like Factor 6


Linkage disequilibrium


Minimum free energy


Polymerase chain reaction


Peroxisome proliferator-activated receptor alpha


Peroxisome proliferator-activated receptor-gamma


Peroxisome proliferator-activated receptor-gamma 2


Quantitative trait loci


Single nucleotide polymorphism


Specificity proteins/Krüppel-like factor


Sterol regulatory element-binding protein 1


Transcription factor


Transcription factor binding site


Untranslated region


Substitution effect


  1. Sun BB, Kurki MI, Foley CN, Mechakra A, Chen C, Marshall E, Wilk JB, Chahine M, Chevalier P, Christé G, Palotie A, Daly MJ, Runz H, Biogen BT. FinnGen: genetic associations of protein-coding variants in human Disease. Nature. 2022;603(7899):95–102.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Pan Z, Yao Y, Yin H, Cai Z, Wang Y, Bai L, Kern C, Halstead M, Chanthavixay G, Trakooljul N, Wimmers K, Sahana G, Su G, Lund MS, Fredholm M, Karlskov-Mortensen P, Ernst CW, Ross P, Tuggle CK, Fang L, Zhou H. Pig genome functional annotation enhances the biological interpretation of complex traits and human Disease. Nat Commun. 2021;12(1):5848.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Zhou Z, Li M, Cheng H, Fan W, Yuan Z, Gao Q, Xu Y, Guo Z, Zhang Y, Hu J, Liu H, Liu D, Chen W, Zheng Z, Jiang Y, Wen Z, Liu Y, Chen H, Xie M, Zhang Q, Huang W, Wang W, Hou S, Jiang Y. An intercross population study reveals genes associated with body size and plumage color in ducks. Nat Commun. 2018;9(1):3971–4.

    Google Scholar 

  4. Li X, Yang J, Shen M, Xie X, Liu G, Xu Y, Lv F, Yang H, Yang Y, Liu C, Zhou P, Wan P, Zhang Y, Gao L, Yang J, Pi W, Ren Y, Shen Z, Wang F, Deng J, Xu S, Salehian-Dehkordi H, Hehua E, Esmailizadeh A, Dehghani-Qanatqestani M, Štěpánek O, Weimann C, Erhardt G, Amane A, Mwacharo JM, Han J, Hanotte O, Lenstra JA, Kantanen J, Coltman DW, Kijas JW, Bruford MW, Periasamy K, Wang X, Li M, Sub CPC, Toxicologie OH, Pharmacology A. Whole-genome resequencing of wild and domestic sheep identifies genes associated with morphological and agronomic traits. Nat Commun. 2020;11(1):2815.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Liu S, Gao Y, Canela-Xandri O, Wang S, Yu Y, Cai W, Li B, Xiang R, Chamberlain AJ, Pairo-Castineira E, Mellow D, Rawlik K, Xia K, Yao C, Navarro Y, Rocha P, Li D, Yan X, Li Z, Rosen C, Van Tassell BD, Vanraden CP, Zhang PM, Ma S, Cole L, Liu JB, Tenesa GE, Fang A. A multi-tissue atlas of regulatory variants in cattle. Nat Genet. 2022;54(9):1438–47.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Lachmann A, Torre D, Keenan AB, Jagodnik KM, Lee HJ, Wang L, Silverstein MC, Ma’Ayan A. Massive mining of publicly available RNA-seq data from human and mouse. Nat Commun. 2018;9(1):1366.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Zhao Y, Hou Y, Xu Y, Luan Y, Zhou H, Qi X, Hu M, Wang D, Wang Z, Fu Y, Li J, Zhang S, Chen J, Han J, Li X, Zhao S. A compendium and comparative epigenomics analysis of cis-regulatory elements in the pig genome. Nat Commun. 2021;12(1):2217.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Thiede BR, Mann ZF, Chang W, Ku YC, Son YK, Lovett M, Kelley MW, Corwin JT. Retinoic acid signalling regulates the development of tonotopically patterned hair cells in the chicken cochlea. Nat Commun. 2014;5:3840.

    Article  CAS  PubMed  Google Scholar 

  9. Dong Y, Xie M, Jiang Y, Xiao N, Du X, Zhang W, Tosser-Klopp G, Wang J, Yang S, Liang J, Chen W, Chen J, Zeng P, Hou Y, Bian C, Pan S, Li Y, Liu X, Wang W, Servin B, Sayre B, Zhu B, Sweeney D, Moore R, Nie W, Shen Y, Zhao R, Zhang G, Li J, Faraut T, Womack J, Zhang Y, Kijas J, Cockett N, Xu X, Zhao S, Wang J, Wang W. Sequencing and automated whole-genome optical mapping of the genome of a domestic goat (Capra hircus). Nat Biotechnol. 2013;31(2):135–41.

    Article  CAS  PubMed  Google Scholar 

  10. MacKay H, Scott CA, Duryea JD, Baker MS, Laritsky E, Elson AE, Garland TJ, Fiorotto ML, Chen R, Li Y, Coarfa C, Simerly RB, Waterland RA. DNA methylation in AgRP neurons regulates voluntary exercise behavior in mice. Nat Commun. 2019;10(1):5364.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Zhang D, Wu S, Zhang X, Ren S, Tang Z, Gao F. Coordinated transcriptional and post-transcriptional epigenetic regulation during skeletal muscle development and growth in pigs. J Anim Sci Biotechno. 2022;13(1):146.

    Article  CAS  Google Scholar 

  12. Zhou Y, Liu S, Hu Y, Fang L, Gao Y, Xia H, Schroeder SG, Rosen BD, Connor EE, Li C, Baldwin RL, Cole JB, Van Tassell CP, Yang L, Ma L, Liu GE. Comparative whole genome DNA methylation profiling across cattle tissues reveals global and tissue-specific methylation patterns. Bmc Biol. 2020;18(1):1–85.

    Article  Google Scholar 

  13. Fan Y, Liang Y, Deng K, Zhang Z, Zhang G, Zhang Y, Wang F. Analysis of DNA methylation profiles during sheep skeletal muscle development using whole-genome bisulfite sequencing. BMC Genomics. 2020;21(1):327.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Zhang J, Han B, Zheng W, Lin S, Li H, Gao Y, Sun D. Genome-wide DNA methylation Profile in Jejunum reveals the potential genes Associated with paratuberculosis in dairy cattle. Front Genet. 2021;12:735147.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Mohamadipoor SL, Mohammadabadi M, Amiri GZ, Babenko O, Stavetska R, Kalashnik O, Kucher D, Kochuk-Yashchenko O, Asadollahpour NH. Signature selection analysis reveals candidate genes associated with production traits in Iranian sheep breeds. Bmc Vet Res. 2021;17(1):369.

    Article  Google Scholar 

  16. Witt KE, Huerta-Sánchez E. Convergent evolution in human and domesticate adaptation to high-altitude environments. Philosophical Trans Royal Soc Lond Ser B Biol Sci. 2019;374(1777):1–9.

    Google Scholar 

  17. Yasumizu Y, Sakaue S, Konuma T, Suzuki K, Matsuda K, Murakami Y, Kubo M, Palamara PF, Kamatani Y, Okada Y, Satta Y. Genome-wide natural selection signatures are linked to genetic risk of modern phenotypes in the Japanese Population. Mol Biol Evol. 2020;37(5):1306–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Almeida OAC, Moreira GCM, Rezende FM, Boschiero C, de Oliveira Peixoto J, Ibelli AMG, Ledur MC, de Novais FJ, Coutinho LL. Identification of selection signatures involved in performance traits in a paternal broiler line. BMC Genomics. 2019;20(1):449.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Liang R, Han B, Li Q, Yuan Y, Li J, Sun D. Using RNA sequencing to identify putative competing endogenous RNAs (ceRNAs) potentially regulating fat metabolism in bovine liver. Sci Rep-Uk. 2017;1(7):6396.

    Article  Google Scholar 

  20. Botella LM, Sanchez-Elsner T, Sanz-Rodriguez F, Kojima S, Shimada J, Guerrero-Esteo M, Cooreman MP, Ratziu V, Langa C, Vary CP, Ramirez JR, Friedman S, Bernabeu C. Transcriptional activation of endoglin and transforming growth factor-beta signaling components by cooperative interaction between Sp1 and KLF6: their potential role in the response to vascular injury. Blood. 2002;100(12):4001–10.

    Article  CAS  PubMed  Google Scholar 

  21. Li D, Yea S, Li S, Chen Z, Narla G, Banck M, Laborda J, Tan S, Friedman JM, Friedman SL, Walsh MJ. Kruppel-like factor-6 promotes preadipocyte differentiation through histone deacetylase 3-dependent repression of DLK1. J Biol Chem. 2005;280(29):26941–52.

    Article  CAS  PubMed  Google Scholar 

  22. Inuzuka H, Nanbu-Wakao R, Masuho Y, Muramatsu M, Tojo H, Wakao H. Differential regulation of immediate early gene expression in preadipocyte cells through multiple signaling pathways. Biochem Bioph Res Co. 1999;265(3):664–8.

    Article  CAS  Google Scholar 

  23. Miele L, Beale G, Patman G, Nobili V, Leathart J, Grieco A, Abate M, Friedman SL, Narla G, Bugianesi E, Day CP, Reeves HL. The Kruppel-Like factor 6 genotype is Associated with Fibrosis in nonalcoholic fatty Liver Disease. Gastroenterology. 2008;135(1):282–91.

    Article  CAS  PubMed  Google Scholar 

  24. Park JH, Eliyahu E, Narla G, DiFeo A, Martignetti JA, Schuchman EH. KLF6 is one transcription factor involved in regulating acid ceramidase gene expression. Biochim Biophys Acta. 2005;1732(1–3):82–7.

    Article  CAS  PubMed  Google Scholar 

  25. Ghaleb AM, Katz JP, Kaestner KH, Du JX, Yang VW. Kruppel-like factor 4 exhibits antiapoptotic activity following gamma-radiation-induced DNA damage. Oncogene. 2007;26(16):2365–73.

    Article  CAS  PubMed  Google Scholar 

  26. Suske G, Bruford E, Philipsen S. Mammalian SP/KLF transcription factors: bring in the family. Genomics (San Diego Calif). 2005;85(5):551–6.

    CAS  Google Scholar 

  27. Kaczynski J, Cook T, Urrutia R. Sp1- and kruppel-like transcription factors. Genome Biol. 2003;4(2):206.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Li D, Yea S, Li S, Chen Z, Narla G, Banck M, Laborda J, Tan S, Friedman JM, Friedman SL, Walsh MJ. Krüppel-like Factor-6 promotes Preadipocyte differentiation through histone deacetylase 3-dependent repression of DLK1. J Biol Chem. 2005;280(29):26941–52.

    Article  CAS  PubMed  Google Scholar 

  29. Raza S, Khan R, Cheng G, Long F, Bing S, Easa AA, Schreurs NM, Pant SD, Zhang W, Li A, Zan L. RNA-Seq reveals the potential molecular mechanisms of bovine KLF6 gene in the regulation of adipogenesis. Int J Biol Macromol. 2022;195:198–206.

    Article  CAS  PubMed  Google Scholar 

  30. Raza S, Khan R, Schreurs NM, Guo H, Gui LS, Mei C, Zan L. Expression of the bovine KLF6 gene polymorphisms and their association with carcass and body measures in Qinchuan cattle (Bos Taurus). Genomics. 2020;112(1):423–31.

    Article  CAS  PubMed  Google Scholar 

  31. Qi W, Chen X, Holian J, Tan CYR, Kelly DJ, Pollock CA. Transcription factors Krüppel-Like factor 6 and peroxisome proliferator-activated Receptor-γ Mediate High glucose-Induced Thioredoxin-Interacting protein. Am J Pathol. 2009;175(5):1858–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Escalona-Nandez I, Guerrero-Escalera D, Estanes-Hernandez A, Ortiz-Ortega V, Tovar AR, Perez-Monter C. The activation of peroxisome proliferator-activated receptor gamma is regulated by Kruppel-like transcription factors 6 & 9 under steatotic conditions. Biochem Bioph Res Co. 2015;458(4):751–6.

    Article  CAS  Google Scholar 

  33. Bechmann LP, Vetter D, Ishida J, Hannivoort RA, Lang UE, Kocabayoglu P, Fiel MI, Muñoz U, Patman GL, Ge F, Yakar S, Li X, Agius L, Lee Y, Zhang W, Hui KY, Televantou D, Schwartz GJ, LeRoith D, Berk PD, Nagai R, Suzuki T, Reeves HL, Friedman SL. Post-transcriptional activation of PPAR alpha by KLF6 in hepatic steatosis. J Hepatol. 2013;58(5):1000–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Wang Y, Kim K, Kim J, Sul HS. Pref-1, a preadipocyte secreted factor that inhibits adipogenesis. J Nutr. 2006;136(12):2953–6.

    Article  CAS  PubMed  Google Scholar 

  35. Wu Z, Wang S. Role of kruppel-like transcription factors in adipogenesis. Dev Biol. 2013;373(2):235–43.

    Article  CAS  PubMed  Google Scholar 

  36. Gibson KD, Dechow CD. Genetic parameters for yield, fitness, and type traits in US Brown Swiss dairy cattle. J Dairy Sci. 2018;101(2):1251–7.

    Article  CAS  PubMed  Google Scholar 

  37. Oliveira Junior GA, Schenkel FS, Alcantara L, Houlahan K, Lynch C, Baes CF. Estimated genetic parameters for all genetically evaluated traits in Canadian holsteins. J Dairy Sci. 2021;104(8):9002–15.

    Article  CAS  PubMed  Google Scholar 

  38. Haile-Mariam M, Pryce JE. Variances and correlations of milk production, fertility, longevity, and type traits over time in Australian Holstein cattle. J Dairy Sci. 2015;98(10):7364–79.

    Article  CAS  PubMed  Google Scholar 

  39. Lelli KM, Slattery M, Mann RS. Disentangling the many layers of eukaryotic transcriptional regulation. Annu Rev Genet. 2012;46:43–68.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Lambert SA, Jolma A, Campitelli LF, Das PK, Yin Y, Albu M, Chen X, Taipale J, Hughes TR, Weirauch MT. Hum Transcription Factors Cell. 2018;4(172):650–65.

    Google Scholar 

  41. Chen P, Liu X, Liu Y, Bao X, Wu Q. ARHGAP18 is upregulated by transcription factor GATA1 promotes the Proliferation and Invasion in Hepatocellular Carcinoma. Appl Biochem Biotech; 2023.

  42. Chen C, Zhang L, Ruan Z. GATA3 encapsulated by Tumor-Associated macrophage-derived extracellular vesicles promotes Immune Escape and Chemotherapy Resistance of Ovarian Cancer cells by upregulating the CD24/Siglec-10 Axis. Mol Pharmaceut. 2023;20(2):971–86.

    Article  CAS  Google Scholar 

  43. Barosso IR, Miszczuk GS, Ciriaci N, Andermatten RB, Maidagan PM, Razori V, Taborda DR, Roma MG, Crocenzi FA, Sanchez PE. Activation of insulin-like growth factor 1 receptor participates downstream of GPR30 in estradiol-17beta-D-glucuronide-induced cholestasis in rats. Arch Toxicol. 2018;92(2):729–44.

    Article  CAS  PubMed  Google Scholar 

  44. Ding Y, Feng Y, Huang Z, Zhang Y, Li X, Liu R, Li H, Wang T, Ding Y, Jia Z, Yang J. SOX15 transcriptionally increases the function of AOC1 to modulate ferroptosis and progression in Prostate cancer. Cell Death Dis. 2022;13(8):673.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Madison BB, McKenna LB, Dolson D, Epstein DJ, Kaestner KH. FoxF1 and FoxL1 link hedgehog signaling and the control of epithelial proliferation in the developing stomach and intestine. J Biol Chem. 2009;284(9):5936–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Liu L, Xiao Y, Huang W, Liu S, Huang L, Zhong J, Jia P, Liu W. ATF1/miR-214-5p/ITGA7 axis promotes osteoclastogenesis to alter OVX-induced bone absorption. Mol Med (Cambridge Mass). 2022;28(1):56.

    Article  CAS  Google Scholar 

  47. Seberg HE, Van Otterloo E, Cornell RA. Beyond MITF: multiple transcription factors directly regulate the cellular phenotype in melanocytes and Melanoma. Pigment Cell and Melanoma Research. 2017;30(5):454–66.

    Article  CAS  PubMed  Google Scholar 

  48. Suda S, Rai T, Sohara E, Sasaki S, Uchida S. Postnatal expression of KLF12 in the inner medullary collecting ducts of kidney and its trans-activation of UT-A1 urea transporter promoter. Biochem Bioph Res Co. 2006;344(1):246–52.

    Article  CAS  Google Scholar 

  49. Hoeth M, Niederleithner H, Hofer-Warbinek R, Bilban M, Mayer H, Resch U, Lemberger C, Wagner O, Hofer E, Petzelbauer P, de Martin R, Reitsma PH. The transcription factor SOX18 regulates the expression of matrix metalloproteinase 7 and guidance molecules in human endothelial cells. PLoS ONE. 2012;7(1):e30982.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Yamagishi T, Ozawa M, Ohtsuka C, Ohyama-Goto R, Kondo T. Evx2-Hoxd13 intergenic region restricts enhancer association to Hoxd13 promoter. PLoS ONE. 2007;2(1):e175.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Li J, Mo M, Chen Z, Chen Z, Sheng Q, Mu H, Zhang F, Zhang Y, Zhi X, Li H, He B, Zhou H, Katoh M. Adenoviral delivery of the EMX2 gene suppresses growth in human gastric cancer. PLoS ONE. 2012;7(9):e45970.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Briata P, Ilengo C, Van DeWerken R, Corte G. Mapping of a potent transcriptional repression region of the human homeodomain protein EVX1. Febs Lett. 1997;402(2):131–5.

    Article  CAS  PubMed  Google Scholar 

  53. Speckmann T, Sabatini PV, Nian C, Smith RG, Lynn FC. Npas4 transcription factor expression is regulated by Calcium Signaling pathways and prevents Tacrolimus-induced cytotoxicity in pancreatic Beta Cells*. J Biol Chem. 2016;291(6):2682–95.

    Article  CAS  PubMed  Google Scholar 

  54. Zhu M, Zhang N, He S. Transcription factor KLF4 modulates microRNA-106a that targets Smad7 in gastric cancer. Pathol - Res Pract. 2019;215(8):152467.

    Article  CAS  PubMed  Google Scholar 

  55. Huang P, Peslak SA, Ren R, Khandros E, Qin K, Keller CA, Giardine B, Bell HW, Lan X, Sharma M, Horton JR, Abdulmalik O, Chou ST, Shi J, Crossley M, Hardison RC, Cheng X, Blobel GA. HIC2 controls developmental hemoglobin switching by repressing BCL11A transcription. Nat Genet. 2022;54(9):1417–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Sperone A, Dryden NH, Birdsey GM, Madden L, Johns M, Evans PC, Mason JC, Haskard DO, Boyle JJ, Paleolog EM, Randi AM. The transcription factor erg inhibits vascular inflammation by repressing NF-kappaB activation and proinflammatory gene expression in endothelial cells. Arterioscler Thromb Vasc Biol. 2011;31(1):142–50.

    Article  CAS  PubMed  Google Scholar 

  57. Wan Y, Kertesz M, Spitale RC, Segal E, Chang HY. Understanding the transcriptome through RNA structure. Nat Rev Genet. 2011;12(9):641–55.

    Article  CAS  PubMed  Google Scholar 

  58. Dethoff EA, Chugh J, Mustoe AM, Al-Hashimi HM. Functional complexity and regulation through RNA dynamics. Nature. 2012;482(7385):322–30.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  59. Zhang Z, Ober U, Erbe M, Zhang H, Gao N, He J, Li J, Simianer H. Improving the accuracy of whole genome prediction for complex traits using the results of genome wide association studies. PLoS ONE. 2014;9(3):e93017.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Brondum RF, Su G, Janss L, Sahana G, Guldbrandtsen B, Boichard D, Lund MS. Quantitative trait loci markers derived from whole genome sequence data increases the reliability of genomic prediction. J Dairy Sci. 2015;98(6):4107–16.

    Article  CAS  PubMed  Google Scholar 

  61. Percie Du Sert N, Hurst V, Ahluwalia A, Alam S, Avey MT, Baker M, Browne WJ, Clark A, Cuthill IC, Dirnagl U, Emerson M, Garner P, Holgate ST, Howells DW, Karp NA, Lazic SE, Lidster K, MacCallum CJ, Macleod M, Pearl EJ, Petersen OH, Rawle F, Reynolds P, Rooney K, Sena ES, Silberberg SD, Steckler T, Würbel H. The ARRIVE guidelines 2.0: updated guidelines for reporting animal research. Plos Biol. 2020;18(7):e3000410.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references


We appreciate Heibei Baoding Qingyuanlvze and Shuangfeng dairy farms for providing the blood samples and pedigree data for the Chinese Holstein cows used in this study.


This work was financially supported by the Hebei Provincial Key Research Projects (21326358D), National Key Research and Development Program of China (2021YFF1000700), and Program for Changjiang Scholars and Innovative Research Team in University (IRT_15R62). All the funding bodies supported the design of this study. The Hebei Provincial Key Research Projects (21326358D) supported the collection of the experimental population, blood samples and phenotypic data. The National Key Research and Development Program of China (2021YFF1000700) and the Program for Changjiang Scholars and Innovative Research Team in University (IRT_15R62) supported the analysis and interpretation of data. All the funding bodies supported the final manuscript.

Author information

Authors and Affiliations



DS conceived and designed the experiments. YL collected blood samples and conducted experiments, association analysis and data analysis. WZ and PP participated in sample collection and data analysis. BH participated in experimental design and data analysis. JN, CY, GJ, YM and JL collected the experimental population, blood samples and analyzed phenotypic and pedigree data. DS, YL and BH prepared the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Junqing Ni or Dongxiao Sun.

Ethics declarations

Ethics approval and consent to participate

All methods were performed in accordance with the relevant guidelines and regulations. All protocols for collection of the blood samples of Chinese Holstein cows were reviewed and approved by the Institutional Animal Care and Use Committee (IACUC) at China Agricultural University (Beijing, China; permit number: DK996). Sample collection specifically for this study followed the standard procedures with the full agreement of the Heibei Baoding Qingyuanlvze and Shuangfeng dairy farms who owned the animals. This study was carried out in compliance with the ARRIVE guidelines 2.0 [61].

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, Y., Han, B., Zheng, W. et al. Identification of genetic associations and functional SNPs of bovine KLF6 gene on milk production traits in Chinese holstein. BMC Genom Data 24, 72 (2023).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: