Determination of Genetic Effects of LIPK and LIPJ Genes on Milk Fatty Acids in Dairy Cattle

In our previous genome-wide association study (GWAS) on milk fatty acids (FAs) in Chinese Holstein, we discovered 83 genome-wide significant single nucleotide polymorphisms (SNPs) associated with milk FAs. Two of them were close to lipase family member K (LIPK) and lipase family member J (LIPJ), respectively. Hence, this study is a follow-up to verify whether the LIPK and LIPJ have significant genetic effects on milk FAs in dairy cattle. By re-sequencing the entire exons, and 3 kb of 5′ and 3′ flanking regions, two and seven SNPs were identified in LIPK and LIPJ, respectively, including a novel SNP, ss158213049726. With the Haploview 4.1 software, we found that five of the SNPs in LIPJ formed a haplotype block (D′ = 0.96 ~ 1.00). Single-locus association analyses revealed that each SNP in LIPK and LIPJ was significantly associated with at least one milk FA (p = < 1.00 × 10−4 ~ 4.88 × 10−2), and the haplotype-based association analyses showed significant genetic effects on nine milk FAs (p = < 1.00 × 10−4 ~ 3.98 × 10−2). Out of these SNPs, the missense mutation in LIPK gene, rs42774527, could change the protein secondary structure and function predicted by SOPMA, SIFT, and PROVEAN softwares. With the Genomatix software, we predicted that two SNPs, rs110322221 in LIPK and rs211373799 in LIPJ, altered the transcription factors binding sites (TFBSs), indicating their potential regulation on promoter activity of the genes. Furthermore, we found that both LIPK and LIPJ had relatively high expressions in the mammary gland. In conclusion, our research is the first to demonstrate that LIPK and LIPJ genes have significant associations with milk FAs, and the identified SNPs might be served as genetic markers to optimize breeding programs for milk FAs in dairy cattle. This research deserves in-depth verification.


Introduction
In the dairy industry, milk production traits are the most important factors when considering breeding goals, especially milk composition, milk protein and fat, which account for main focus of genetic improvement programs designed to increase the productivity of dairy cattle. Fatty acid,

Animal Population and Phenotypic Data
In this study, we used a total of 1065 Chinese Holstein cows from 44 sires, which is a different population from those used in our previous GWAS. The cows were from 23 dairy farms belonging to the Sanyuanlvhe Dairy Farming Center (Beijing, China), where the Dairy Herd Improvement system (DHI) has been carried out since 1999. In the dairy industry, DHI means a monthly routine standard milk production measurement system that is used to test the phenotypic values of milk protein percentage, fat percentage, and somatic cell count etc. using individual milk samples of dairy cows. During November to December of 2014, 50 mL milk samples were collected from the 1065 cows to be used for measuring DHI in the Beijing Dairy Cattle Center (www.bdcc.com.cn). After DHI measuring, we used 2 ml milk samples to measure the milk FAs by gas chromatography.

SNPs Identification and Genotyping
Based on the genomic sequences of LIPK (GenBank accession number: AC_000183.1) and LIPJ (GenBank accession number: AC_000183.1), we respectively designed 21 and 23 pairs of primers using the Primer 3 (http://bioinfo.ut.ee/primer3-0.4.0/) in the entire exons with their partial adjacent introns, and 3 kb of 5 and 3 flanking regions of the two genes (Supplementary Table S2). These primers were synthesized at the Beijing Genomics Institute (Beijing, China). DNA was extracted from each semen sample of 44 sires using a salting-out procedure, and the quantity and quality were measured by a NanoDrop 2000 spectrophotometer (Thermo Scientific, Hudson, DE, USA) and gel electrophoresis, respectively. We randomly mixed them into two pools used for all the polymerase chain reaction (PCR) amplifications with equal concentrations (50 ng/µL) for each DNA, and each pool included 22 sires. PCR amplifications for the two pooled DNA were performed in a final reaction volume of 25 µL, which consisted of 2 µL genomic DNA, 1.25 µL of each primer (10 mM), 12.5 µL Premix Taq TM (Takara, Dalian, China), and 8 µL DNase/RNase-Free Deionized Water (Tiangen, Beijing, China). The amplification conditions were as follows: 5 min at 94 • C for initial denaturing followed by 35 cycles at 94 • C for 30 s, 60 • C for 30 s, 72 • C for 30 s and a final extension at 72 • C for 7 min. Then each PCR product was sequenced by ABI3730XL DNA analyzer (Applied Biosystems, Foster, CA, USA), and the sequences were aligned with the reference genome (UMD3.1.1) using the BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi) to search the potential SNPs.
DNAs of the whole blood samples from 1065 Chinese Holstein cows were extracted using TIANamp Blood DNA Kit (Tiangen, Beijing, China) according to the manufacturer s instructions. We measured the quantity and quality of the blood DNA using a NanoDrop 2000 spectrophotometer (Thermo Scientific, Hudson, DE, USA) and gel electrophoresis, respectively. The identified SNP were genotyped using the blood DNAs of these cows by the matrix-assisted laser desorption/ionization time of flight mass spectrometry assay (MALDI-TOF MS, Sequenom MassARRAY, Agena, San Diego, USA).

Estimation of linkage disequilibrium (LD)
We performed the extent of LD among the nine SNPs of LIPK and LIPJ identified in this study using the Haploview 4.1 (Broad Institute, Cambridge, MA, USA), and the haplotype blocks were estimated.
In addition, we also used the Haploview 4.1 to estimate the LD among the SNPs of LIPK, LIPJ, and SCD genes, including two in LIPK, seven in LIPJ, and 24 in SCD (Table S4). The SNPs (one in LIPK, one in LIPJ, and 24 in SCD) of the previous studies were not genotyped in this study. Therefore, we used the 1577 Holsteins with worldwide correlation from the 1000 Bull Genomes Project [22] to estimate the LD.

Association Analyses
For the association analyses between the identified SNPs or haplotype blocks and 24 milk FAs, 3335 individuals had their pedigree traced back to three generations. The SAS 9.2 mixed procedure was employed for the analyses based on the following model: in which, Y ijklm was the phenotypic value of each trait of per cow; µ was the overall mean; G i was the fixed effect corresponding to the genotype or haplotype combination; h j and l k were the fixed effects of the farm and stage of lactation, respectively; a l was the random polygenic effect; M m was the fixed effect of calving month; b was the regression coefficient of covariate M; and e ijklm was the random residual. In addition, the additive (a), dominance (d), and allele substitution (α) effects were estimated by the following equations [23]: , and α = α + d(q − p), where AA, AB, and BB represented the least square mean of milk fatty acids corresponding to the genotypes, and p and q were the allele frequencies of A and B, respectively.

Protein Structure and Function Prediction
We used the NPSA SOPMA SERVER (https://npsa-prabi.ibcp.fr/cgi-bin/npsa_automat.pl? page=/NPSA/npsa_sopma.html) to predict the change of protein secondary structure caused by the missense mutation in the coding regions of genes, and the parameters included window width (17), similarity threshold (8), and number of states (4). In addition, we assessed whether the missense mutation altered the protein function using the SIFT (http://sift.bii.a-star.edu.sg/) and PROVEAN (http://sift.jcvi.org/index.php) software. The score thresholds of the SIFT and PROVEAN were 0.05 and −2.5, respectively. If the score was below the thresholds, the SNP was predicted as a "deleterious" variant.

Predication of Transcription Factors Binding Sites (TFBSs)
We used the Genomatix software (http://www.genomatix.de/cgi-bin/sessions/login.pl?s= 77bfbe2f9849561b2b3e91c76f124365) to predict the binding sites of the transcription factors (TFs) caused by the SNPs in 5 flanking and untranslated region (UTR; matrix similarity threshold, MST > 0.80). During the prediction, the matrix exhibited a high conservation at this position (Ci-value >60%).

Gene Expressions Assay of LIPK and LIPJ
To further detect the potential function of LIPK and LIPJ, we detected expression levels in eight types of tissue, including mammary, uterus, ovary, liver, spleen, heart, kidney, and lung. Three lactating Chinese Holstein cows were selected from the Sanyuanlvhe Dairy Farming Center (Beijing, China). Eight types of tissue were collected from each cow and then stored at liquid nitrogen.
We isolated the total RNAs from the eight types of tissue from each cow using Trizol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer s instructions. Subsequently, we measured the quantity and quality of RNA with the NanoDrop 2000 spectrophotometer (Thermo Scientific, Hudson, DE, USA) and gel electrophoresis, respectively. We used PrimerScriptH RT reagent Kit (TaKaRa Biotechnology Co., Ltd., Dalian, China) for the reverse transcription. We conducted the quantitative real-time PCR (qRT-PCR) using SYBR green fluorescence (Roche, Penzberg, Germany) with a volume of 15 µL consisted of 2 µL template of cDNA, 0.375 µL of each primer (10 mM), 4.75 µL distilled water, and 7.5 µL SYBR Green Mixture. The PCR conditions were: Denaturation 95 • C for 10 min; amplification 45 cycles at 95 • C for 10 s, 58 • C for 10 s, and 72 • C for 10 s. The qRT-PCR primers of LIPK, LIPJ, and GAPDH are shown in Table S4. We performed all the measurements in triplicate and the relative gene expression was normalized by the GAPDH with 2 −∆∆Ct method [24].

SNPs Identification
We identified two SNPs for LIPK gene, namely, rs110322221 in the 5 UTR and rs42774527 in the exon 11 (Table 1). For LIPJ gene, seven SNPs were identified (Table 1), including rs41606812 in the 5 flanking region, rs211373799 in the 5 UTR, rs42107056 in the 3 UTR, and four SNPs (rs42107122, ss158213049726, rs209219656, and rs42107114) in the 3 flanking region. Among the total nine SNPs, rs42774527 was a missense mutation with the substitution of amino acid from threonine (ACA) to lysine (AAA), and ss158213049726 was first identified in this study.

Linkage disequilibrium among the SNPs of LIPK, LIPJ and SCD Genes
To determine whether the genetic associations of LIPK and LIPJ on FAs were caused by the effect of SCD, we indirectly estimated the LD among the SNPs of LIPK, LIPJ, and SCD genes using the 1000 Bull data [22], and three haplotype blocks were observed (D = 0.93~1.00; Figure 2). The haplotype block 1 consisted of eight SNPs, including rs41606812, rs211373799, rs42107056, rs42107122, rs209219656, rs42107114, and rs209033376 of LIPJ and rs110933619 of LIPK. The haplotype block 2 and 3 were formed by six (rs211483324, rs41255693, rs41255692, rs41255691, rs41255690, and rs41255688) and four SNPs (rs42086690, rs42087679, rs42088948, and rs42088972) of SCD gene, respectively. The results indirectly showed that there was no linkage (r 2 > 0.8) between LIPK/LIPJ and SCD, indicating that the genetic associations of LIPK/LIPJ with milk FAs were not due to the LD with SCD.

rs42774527 Caused the Changes of the LIPK Protein Structure and Function
We used the SOPMA software to predict the changes of the protein secondary structure for the missense mutation (rs42774527) in the exon 11 of LIPK gene, and found that the alpha helix was changed from 34.09% to 33.08%, extended strand from 19.95% to 19.44%, beta turn from 6.06% to 6.31%, and random coil from 39.90% to 41.16%. In addition, the SNP, rs42774527, was considered as a "deleterious" mutation by SIFT (score = 0.04) and PROVEAN (score = −3.135), that potentially altered the protein functions.

Transcription Factors Binding Sites (TFBSs) Changed by rs110322221 and rs211373799
For the SNPs in 5 flanking and UTR of LIPK and LIPJ genes, we predicted the changes of TFBSs caused by them. The allele A of rs110322221 in the 5 UTR of LIPK gene was predicted to create the binding sites for the TF FAC1 (fetal alz-50 clone 1; MST = 0.97) (Figure 4). For the LIPJ gene, the allele C in the rs211373799 could create the binding sites for the TFs AIRE (autoimmune regulator; MST = 0.81) and MTBF (muscle-specific Mt binding site; MST = 0.91), and the allele A in this SNP could invent the TFBS for FAST1 (FAST-1SMAD interacting protein; MST = 0.87) (Figure 3). . D′ is the value of D prime between the two loci. Haplotype block 1 included rs41606812, rs211373799, rs42107056, rs42107122, rs209219656, rs42107114, and rs209033376 of LIPJ and rs110933619 of LIPK. Haplotype block 2 and 3 included six (rs211483324, rs41255693, rs41255692, rs41255691, rs41255690, and rs41255688) and four SNPs (rs42086690, rs42087679, rs42088948, and rs42088972) of SCD gene, respectively.

Transcription Factors Binding Sites (TFBSs) Changed by rs110322221 and rs211373799
For the SNPs in 5′ flanking and UTR of LIPK and LIPJ genes, we predicted the changes of TFBSs caused by them. The allele A of rs110322221 in the 5′ UTR of LIPK gene was predicted to create the binding sites for the TF FAC1 (fetal alz-50 clone 1; MST = 0.97) (Figure 4). For the LIPJ gene, the allele C in the rs211373799 could create the binding sites for the TFs AIRE (autoimmune regulator; MST = 0.81) and MTBF (muscle-specific Mt binding site; MST = 0.91), and the allele A in this SNP could invent the TFBS for FAST1 (FAST-1SMAD interacting protein; MST = 0.87) (Figure 3).

Expressions of LIPK and LIPJ Genes
By using the qRT-PCR method, we observed that LIPK and LIPJ genes were expressed in eight types of tissue, including mammary, uterus, ovary, liver, spleen, heart, kidney, and lung. The relatively high expression levels of these two genes in the mammary gland implied that they might be involved in milk composition synthesis (Figure 4).

Expressions of LIPK and LIPJ Genes
By using the qRT-PCR method, we observed that LIPK and LIPJ genes were expressed in eight types of tissue, including mammary, uterus, ovary, liver, spleen, heart, kidney, and lung. The relatively high expression levels of these two genes in the mammary gland implied that they might be involved in milk composition synthesis ( Figure 4).

Discussion
Our previous GWAS identified that LIPK and LIPJ genes were both promising candidates for FAs in dairy cattle. In this study, we first demonstrated that both LIPK and LIPJ genes were mainly associated with saturated medium-and long-chain milk FAs.
SFAs are reported to increase the risks for diseases, while UFAs lower the risks [25,26]. C6:0 is related with the depressive symptoms phenotype [27]. C14:0 causes the nonalcoholic steatohepatitis associated with lipodystrophy [28]. C14:1 inhibits the formation of large multinucleated osteoclasts related with osteoporosis and other bone-lytic diseases [29]. C17:1 can be affected by the fatty acid two-hydroxylase to regulate intestinal homoeostasis in elegans [30]. C18:1cis-9 can reduce inflammation [31]. To date, no report has highlighted the correlation between LIPK and LIPJ and milk fatty acids in dairy cattle. Our association analysis results showed that genotypes AA of rs110322221 and AA of rs41606812 decreased milk contents of C8:0 and C20:0, respectively. Genotypes TT of rs42107122 deceased C14:0 compared with genotypes CC. Both the genotypes CC of rs42107056 and GG of ss158213049726 showed lower C14:0, and higher C17:1. Genotype TT of rs209219656, and haplotype H2 formed by five SNPs (rs211373799, rs42107056, rs42107122, ss158213049726 and rs209219656) increased C17:1. In the present study, we further verified that no LD (r 2 > 0.8) was found between LIPK/LIPJ and SCD, implying that the associations of LIPK and LIPJ with milk FAs were not caused by the linkage to the SCD. In addition, we observed the relative high mRNA level of LIPK and LIPJ in the mammary gland, suggesting that the two genes might be related to the composition of milk fat. Furthermore, we observed the expressions of the two genes in the uterus, ovary, liver, spleen, heart, kidney, and lung. There has been no report about the expression of LIPK and LIPJ in these tissues until now. A recent study showed the LIPK expression level in granular keratinocytes of human abdominal skin, and suggested that the gene might play a highly specific role in the last step of keratinocyte differentiation to regulate the lipid metabolism of the most differentiated epidermal layers [32]. We hypothesize that LIPK and LIPJ could play roles in some biological functions specific for these types of tissue. Considering the already identified significant genetic effects of LIPK and LIPJ on milk fatty acids, the six SNPs, including rs110322221 of LIPK, rs41606812, rs42107056, rs42107122, ss158213049726, and rs209219656 of LIPJ, could be used as molecular markers for genetic improvement programs of dairy cattle through marker-assisted selection to provide the market with better milk with high UFA and low SFA contents. SNP, rs42774527, as a missense mutation, was predicted to alter the protein secondary structure, and was a "deleterious" mutation changing protein functions. The association analysis showed that the cows with the AA genotype of rs42774527 had significantly higher C14:0, and lower C6:0, C17:0, C18:1cis-9, and total index than those with CC genotype. Hence, the changes in LIPK protein's secondary structure might impact the LIPK protein's ability to form milk FA traits.
Two SNPs, rs110322221 of LIPK and rs211373799 of LIPJ, in their 5 UTR created the TFBSs for FAC1, AIRE, MTBF, or FAST1. The association analysis showed that the cows with the genotype AA of rs110322221 had significant lower C8:0 than those with genotype GG, and the genotype CC of rs211373799 significantly lowered the C6:0 and C8:0, and increased the C20:0, than the genotype AA. We didn't find any reports about the regulation of TFs, FAC1, AIRE, MTBF, and FAST1 on lipid metabolism in the mammary gland. Previous studies showed that FAC1 represses transcription of the APP gene to influence neuronal development and neuro-degenerative diseases in human [33], AIRE transcriptionally activates the IL-6 in androgen-independent cells to modulate the prostate tumor microenvironment in mice [34], and FAST1 directly regulates mesodermal gene expression to mediate the form of mesodermal in the Xenopus embryo [35]. We herein deduced that FAC1, AIRE, MTBF, and FAST1 could up-/down-regulate the expression of LIPK or LIPJ to affect the milk FA traits.
There was a limitation in this study. We used DNA pooling to search SNPs, which may have resulted in rare SNPs with very low allele frequencies being missed. To identify as many possible variants in our population as possible, we randomly pooled the 44 sires into two pools in order to ensure the concentration (50 ng/µL) for each sire was high enough for the identification of potential SNPs. Moreover, we also have performed PCR product sequencing based on pooled DNA to identify SNPs in several previous studies [36][37][38][39][40][41], and we found that this method can detect the SNP with a low frequency of 0.0419. In addition, as for the variant with a very low frequency, the smaller number of individuals with the homozygous genotype of rare alleles could impact the statistical power of the association analysis result.

Conclusions
Our initial GWAS identified LIPK and LIPJ genes as two promising candidates for FAs in dairy cattle. In the study, we further demonstrated that both LIPK and LIPJ showed significant associations with milk FA traits in the Chinese Holstein population. These SNPs, including rs110322221 of LIPK, rs41606812, rs42107056, rs42107122, ss158213049726, and rs209219656 of LIPJ, and the haplotype H2 could be used as genetic markers to reduce the milk SFA or increase milk UFA.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4425/10/2/86/s1, Table S1: Summary of the phenotype values. Table S2: PCR primers of LIPK and LIPJ genes. Table S3: SNPs used for the LD analysis between LIPK, LIPJ, and SCD. Table S4: Primers used for the amplification of bovine LIPK, LIPJ, and GAPDH genes by quantitative real-time PCR. Table S5: Associations of nine SNPs of LIPK, and LIPJ genes with fatty acid traits in dairy cattle (Least square mean ± Standard error).