The Functional SNPs in the 5’ Regulatory Region of the Porcine PPARD Gene Have Significant Association with Fat Deposition Traits

Peroxisome proliferator-activated receptor delta (PPARD) is a key regulator of lipid metabolism, insulin sensitivity, cell proliferation and differentiation. In this study, we identified two Single Nucleotide Polymorphisms (SNPs, g.1015 A>G and g.1018 T>C) constituting four haplotypes (GT, GC, AC and AT) in the 5’ regulatory region of porcine PPARD gene. Functional analysis of the four haplotypes showed that the transcriptional activity of the PPARD promoter fragment carrying haplotype AC was significantly lower than that of the other haplotypes in 3T3-L1, C2C12 and PK-15 cells, and haplotype AC had the lowest binding capacities to the nuclear extracts. Transcription factor 7-like 2 (TCF7L2) enhanced the transcription activities of promoter fragments of PPARD gene carrying haplotypes GT, GC and AT in C2C12 and 3T3-L1 cells, and increased the protein expression of PPARD gene in C2C12 myoblasts. TCF7L2 differentially bound to the four haplotypes, and the binding capacity of TCF7L2 to haplotype AC was the lowest. There were significant associations between -655A/G and fat deposition traits in three pig populations including the Large White × Meishan F2 pigs, France and American Large White pigs. Pigs with genotype GG had significantly higher expression of PPARD at both mRNA and protein level than those with genotype AG. These results strongly suggested that the SNPs in 5’ regulatory region of PPARD genes had significant impact on pig fat deposition traits.

We are interested in studying pig fat deposition because it plays important roles in animal agriculture and can be used as biomedical model for human obesity. Pig carcass composition such as backfat thickness and lean meat percentage are controlled by polygenes with pleiotropic effects. Identification of these polygenes or linked markers is necessary for understanding the genetic basis of carcass traits and the application of marker assisted selection (MAS) in breeding programs [22]. A significant quantitative trait loci (QTL) for backfat thickness (BFT) has been consistently mapped to sus scrofa chromosomes 7 (SSC7) p1. 1-q1.4 in several pig populations [23][24][25][26][27][28][29][30][31]. Using the F 2 resource population derived from the intercross of Large White boars and Meishan dams, we have identified significant QTLs for carcass and meat quality traits on SSC7 [23,32]. Haplotypes of the porcine PPARD gene are associated with backfat thickness [33], and PPARD gene also functions in fatty acid metabolism and fat metabolism [34]. Therefore, PPARD gene was considered to be a promising positional candidate gene for the fat deposition traits. Previous studies also showed that one missense mutation within the coding sequence of PPARD gene was significantly associated with ear size [35,36]. However, whether and how the genetic variants within upstream regulatory region of PPARD gene affected fat deposition traits were still unknown in pigs. In this study, we identified two functional SNPs (GU565976.1: g.1015 A>G and g.1018 T>C) in the 5' regulatory region of PPARD which altered the binding capacity transcription factor TCF7L2 to the promoter region, and found that the SNP g.1015A>G was significantly associated with fat deposition traits in three pig populations.  [37]. The fat deposition and carcass traits were as following: backfat thickness at shoulder (BFT1, cm), backfat thickness at thorax-waist (BFT2, cm), backfat thickness at buttock (BFT3, cm), average backfat thickness at shoulder, thorax-waist and buttock (ABT, cm), leaf fat wight (LFW, kg), backfat thickness between 6 th and 7 th ribs (67RIBBF, cm), rib number (RN), carcass length from the first cervical vertebra to anterior border of pubic symphysis (CL1, cm), carcass length from the first thoracic vertebra to anterior border of pubic symphysis (CL2, cm), Internal fat rate (IFR). The live backfat thickness were measured by B-ultrasound machine (ESAOTE, Mylob Touch VET) between 85 kg and 115 kg body weight and then corrected to the backfat thickness at 100 kg body weight according to the national standard of swine performance testing of P. R. China (NY/T822-2004). Sequencing and haplotypes inference were conducted in 90 pigs from seven breeds, including 9 Chinese Bamei pigs, 6 Erhualian pigs, 4 Huainan pigs, 5 Duroc pigs, 20 American Large White pigs, 25 France Large White pigs, and 21 Large White × Meishan F 2 pigs.

SNP detection and haplotype inference
In order to detect SNPs in the PPARD 5' regulatory region, we sequenced the promoter region (GenBank accession no. GU565976.1) from five Meishan pigs and five Large White pigs, respectively. PCR products of 544 bp, 568 bp and 757 bp spanning nucleotides − 1415 bp to 464 bp (1880 bp length) were sequenced to detect the SNPs in the 5' regulatory region. The mix of PCR contained 1.5 mM MgCl 2 . The primers for DNA sequencing, amplification and annealing temperature are listed in S1 Table. The purified PCR products of five Meishan and five Large White pigs were cloned into the pMD18-T vector (TaKaRa, Japan), and sequenced commercially (Sangon, China). Then these sequences were aligned using Clustalw (http:// www.ebi.ac.uk/Tools/msa/clustalw2/) to detect the SNPs. Direct sequencing of 568 bp PCR products containing the SNPs (g.1015 A>G and g.1018 T>C) were further conducted to infer the haplotypes. Haplotypes were inferred using Haploview software (Version 4.2) [38].
After transfection for 24 h, cells were lysed in 100 μl of lysis buffer, and firefly and renilla luciferase activities were measured using the Dual-Luciferase Reporter Assay System (Promega, Madison, WI, USA). The luciferase activity was measured using PerkinElmer 2030 Multilabel Reader (PerkinElmer, USA). To normalize the transfection efficiency, the pRL-TK plasmid vector (Promega, USA) was co-transfected with the reporter construct as described above. The experiments were performed in three replicates for each construct and the data were the average of three replicates.

Quantitative real-time PCR
For tissue expression pattern analysis, nine tissues including heart, liver, spleen, lung, kidney, stomach, intestines, longissimus dorsi (LD), and backfat were sampled from three France Large White pigs. For differential gene expression analysis between different genotypes, muscle and adipose tissues were sampled from three France Large White pigs with genotypes AG and GG, respectively. RNA samples from different tissues of the France Large White pigs were isolated using the TRIzol reagent (Invitrogen, USA). Reverse transcription of RNA (1 μg) was performed using random primers (Invitrogen, USA) and Moloney Murine Leukemia Virus (M-MLV) Reverse Transcriptase (Invitrogen, USA) according to the manufacturer's instructions. The primers used to amplify the cDNA corresponding to PPARD, TCF7L2, and β-actin genes are listed in S3 Table. Quantitative RT-PCR was performed on a LightCycler 480 II Real-Time (Roche, Switzerland) using the FastStart DNA Master SYBR Green I reagent (Roche, Switzerland). The gene expression results were normalized with regard to the expression of the β-actin. The Ct (2 -ΔΔCt ) method was used to analyze the relative gene expression [39].

Electrophoretic mobility shift assay (EMSA)
Nuclear extracts from C2C12 myoblasts, 3T3-L1 cells and pig longissimus muscle were prepared by Nuclear Protein Extraction Kit (BestBio, China). The protein concentration was determined by Bicinchoninic Acid assay (Beyotime, China). The double-stranded DNA probes of four haplotypes were prepared by annealing the desired sense and anti-sense oligonucleotides which were 5' end-labeled with biotin. The DNA binding activity of protein was detected by chemiluminescent EMSA Kit (Beyotime, China). The protein-DNA complexes were analyzed by electrophoresis in 6.0% polyacrylamide gels at 100 V for 1 h in 0.5 TBE running buffer (44.58 mM Tris Base, 44.58 mM Boric Acid, 1.25 mM Na 2 EDTAH 2 O), and then were transferred to a nylon membrane. The dried nylon was visualized using the ECL (Bio-Rad, USA). The bands were detected using the the ImageQuant LAS 4000 mini (GE, USA) according to the manufacturer's instruction.

DNA pull-down
The probes for DNA pull-down assay were the same as EMSA assay. The probes were conjugated with M-280 Streptavidin Dynabeads (Invitrogen, USA) in binding buffer (10 mM Tris-HCl, pH 7.5, 50 mM KCl, 1 mM MgCl 2 , 1 mM EDTA, pH 8.0, 1 mM Na3VO4, 5 mM DTT, 5% glycerol, 0.3% NP-40) for 40 min at room temperature. Non-denaturing total proteins were extracted from 3T3-L1cells using Western and IP Cell lysis Buffer (Sangon Biotech, China). The proteins (400 μg) were incubated with unconjugated Dynabeads for 2 h at 4°C using Dynal MPC-S magnetic particle concentrator (Dynal Biotech, Norway), and the non-specific binding proteins were removed. The supernatant (supernatantI) was collected for binding reaction. The probes were incubated with the Dynabeads for 5 h at 4°C, and then the Dynabeads-probes complexes were washed three times with lysis buffer and incubated with supernatant I for 2 h at 4°C in the presence of 25 μg/ml poly (dI:dC) which could prevent the nonspecific binding of protein with DNA. The precipitates (Dynabeads-probes-protein) were washed with lysis buffer, and then were eluted in SDS sample buffer. The supernatantII was collected for β-actin detection. The precipitates and supernatant were assayed by western blotting.

SNP genotyping and statistical analysis
The SNP (g.1015A>G) was genotyped by PCR-RFLP using EcoRI restriction enzyme. 8.5 μl PCR products were digested with 5 U EcoRI restriction enzyme (Fermentas, Canada) for 6 h at 37°C, and then were separated by electrophoresis on a 1.5% agarose gel with ethidium bromide in 1× TAE buffer.
The association of SNP genotypes with fat deposition and carcass traits in Large White × Meishan F 2 pigs, France, and American Large White pigs were conducted with the general linear model (GLM) procedure of SAS version 8.0 (SAS Institute, Cary, NC, 2000). Both additive and dominance effects were estimated using the REG procedure of SAS version 8.0. The data are presented as mean ± S.D., and statistical significance level was set at p < 0.05.
The model of association analysis in 274 Large White × Meishan F 2 pigs was as following: Where, T ijkl is the observed values of a given trait; μ is the overall mean; S i is effect of sex (i = 1 for male or 2 for female); Y j is the effect of year (j = 1 for year 2000 or 2 for year 2003); G k is the effect of genotype (k = AA, AG and GG); F l is the effect of family (l = 37); b ijkl is the regression coefficient of the slaughter weight for carcass traits, X ijkl is the slaughter weight; e ijkl is the random residual.
The model of association analysis in 710 France Large White pigs and 409 American Large White pigs was as following: Where, T ijkl is the corrected backfat thickness at 100 kg live weight; G k is the effect of genotype (k = AA, AG and GG); S i is the effect of sex (i = 1 for male or 2 for female); F j is effect of family (j = 17 for France Large White and j = 20 for American Large White); B l is the effect of batch (l = 20 for France Large White and l = 8 for American Large White); e ijkl is the random residual.

Results
Identification of two SNPs in the upstream of the 5' regulatory region of PPARD gene Two SNPs, g.1015A>G and g.1018 T>C (GenBank accession no. GU565976.1) were detected at the upstream − 655 bp (A/G) and − 652 bp (T/C) (relative to the transcription start site + 1) of the 5'regulatory region of PPARD gene (Fig 1). The corresponding variation ID of two SNPs in the SNP database (http://www.ncbi.nlm.nih.gov/SNP) was rs80912557 and rs80964726, respectively. The haplotypes of two SNPs were inferred in 90 pigs randomly selected from seven breeds including Bamei, Erhualian, Huainan, Duroc, American Larger White, France Larger White, and Large White × Meishan F 2 pigs. In these pig populations, the two SNPs were in strong linkage disequilibrium (D = 0.925, R 2 = 0.633) and the frequency of haplotypes GT, GC, AC and AT was 52.9%, 8.9%, 36.7% and 1.6%, respectively. In American Large White pigs, the frequency of haplotypes GT, GC, AC and AT was 36.7%, 6.7%, 57.7%, and 0%, respectively. In France Large White pigs, the frequency of haplotypes GT, GC, AC and AT was 86.7%, 10.0%, 3.3%, and 0%, respectively. In Large White × Meishan F 2 pigs, the frequency of haplotypes GT, GC, AC and AT was 32.1%, 25.9%, 42.9%, and 0%, respectively. The results suggested that the allele A at rs80912557 of the PPARD was mostly linked with allele C at rs80964726, and allele G at rs80912557 was linked with allele T or C at rs80964726 in the pig population studied.
The SNPs in the promoter of PPARD gene significantly affected its transcriptional activity The bioinformatics prediction of transcription factor binding sites showed that five transcription factors could bind to the segment containing the two SNPs, and the prediction scores of the binding capacity of transcription factors with different haplotype promoter fragments are listed in S4 Table. According to the bioinformatics prediction of transcription factor binding sites by TESS and Promoter 2.0 software, we constructed a series of deletion fragments by fusing different fragments of the 5' flanking region to luciferase reporter vector pGL3-basic ( Fig  1). All constructs were transfected into C2C12 myoblasts and their promoter activities were examined. The luciferase activity of the promoter fragments PGL3-260 (+ 205 bp to + 464 bp) was the lowest in all the luciferase activities (Fig 2A), suggesting that this region between nucleotides − 100 bp to + 205 bp was the core promoter region. The transcription activity of the promoter fragments PGL3-1077 (− 612 bp to + 464 bp) was lower than that of the other two promoter fragments PGL3-1880GT (− 1415 bp to + 464 bp), and PGL3-1146GT (− 681 bp to + 464 bp). The results showed that the promoter fragments from position − 1415 bp to − 612 bp contained the potential positive regulatory elements. To determine the effects of the SNPs on the promoter transcriptional activities, we designed and studied two sets of promoter luciferase report gene vectors carrying the four haplotypes constructed by site-directed mutagenesis using pGL3-1146GT and pGL3-1880GT as templates. The promoter vectors were transfected into C2C12, 3T3-L1, and PK cells, respectively. In C2C12 myoblasts, the transcriptional activity driven by pGL3-1146AC decreased 8.5-fold, 12.9-fold, and 9.7-fold compared with that driven by pGL3-1146AT, pGL3-1146GT and pGL3-1146GC, respectively (Fig 2B). The same trends were found in 3T3-L1 cells and PK cells. The luciferase expression level driven by pGL3-1880AC was also the lowest among the four haplotype constructs (Fig 2C and 2D). The results indicated that the SNPs in the promoter of PPARD gene significantly affected its transcription activities. The transcription activities of PPARD gene carrying haplotype AC were the lowest and these of haplotype GT the highest.

Potential transcription factors binding to the promoter fragments carrying the four haplotypes
We used genomatrix (http://www.genomatrix.de/) to investigate the potential transcription factors in the promoter region harboring the SNPs identified. The results showed that there were potential ETS1, MEL1, Foxh1, POU2F1, and TCF7L2 binding sites in the promoter region, and the binding capacities of Foxh1, POU2F1 and TCF7L2 to the promoter changed when the SNPs were introduced to the DNA sequences (S4 Table). EMSA was performed to detect the binding capacities of different haplotypes to the nuclear extracts in C2C12 myoblasts, 3T3-L1 cells and pig LD muscle. The results of EMSA showed that the binding capacity of haplotype AC to the nuclear extracts was the weakest in C2C12 myoblasts, 3T3-L1 cells and pig LD muscle (Fig 3A, 3B and 3C), confirming that there are potential transcription factors differentially binding to the promoter fragments carrying the different haplotypes. EMSA was further conducted in C2C12 myoblasts with the probes of haplotypes AC and GC which showed the largest difference in the binding capacities (Fig 3D). The incubation of nuclear extracts of C2C12 myoblasts with probe GC formed a DNA-protein complex. The complex became weaker with 1 × cold probes in the mixture. The complex did not change in the mutation cold probe reaction. By contrast, the incubation of nuclear extracts of C2C12 myoblasts with probe AC did not form the DNA-protein complex. These results demonstrated that transcription factors can bind to the DNA sequence with haplotype AT, GT and GC, but not to haplotype AC.

The PPARD gene SNPs showed different interactions with transcription factor TCF7L2
The binding motif of transcription factor TCF7L2 contains AC/GA/TTCAAAG [40], which is similar to the reverse sequence of the region harboring the haplotypes. In the present study, a classic motif of TCF7L2 was not identified in PPARD promoter, but a variant TCF7L2 binding sequences ccgggATCAAACtcaca was located at − 660/− 644 bp. TCF7L2 could bind to the promoter of PPARD and increase its expression in human colorectal cancer (CRC) cells [41]. Therefore, we hypothesized that SNPs alters transcriptional activities of the promoter through affecting the binding capacities of TCF7L2 to its binding sites. We initially studied the tissue  expression patterns of TCF7L2 and PPARD genes in ten tissues including heart, liver, spleen, lung, kindey, brain, stomach, intestine, LD, and backfat of the four-month-old Large White pigs. The tissue expression profiles of two genes were very similar. They were all broadly expressed in the ten tissues, and highly expressed in liver and adipose which are active in fat synthesis (Fig 4A). To investigate the effects of TCF7L2 overexpression on the promoter transcriptional activities of PPARD gene, we constructed the eukaryotic overexpression vector of the porcine TCF7L2 (pCDNA3.1-TCF7L2) and co-transfected with the promoter fragments carrying the four haplotypes into C2C12 myoblasts and 3T3-L1 cells. TCF7L2 overexpression significantly improved the transcriptional activities of the promoter fragments with haplotypes AT, GT, and GC, whereas there were not significant effects with haplotype AC (Fig 4B). The protein expression levels of PPARD were also up-regulated by TCF7L2 overexpression in C2C12 myoblasts (Fig 4C and 4D). The binding capacities of TCF7L2 to the promoter fragments carrying four haplotypes were further analyzed by DNA pulldown in 3T3-L1 cells. The binding capacities of the fragments carrying haplotypes AT, GT and GC to TCF7L2 were significantly stronger than those of haplotype AC, and the binding capacity of haplotypes GT was the strongest. These results were in agreement with the results of the luciferase assay (Fig 5).
The SNP g.1015A>G was significantly associated with fat deposition traits in three pig populations We used EcoRIPCR-RFLP to genotype the polymorphism g.1015A>G. Digestion of the PCR fragments with EcoRIrestriction enzyme generated one fragment (genotype GG, 568bp), two fragments (genotype AA, 314 bp + 254 bp) or three fragments (genotype AG, 568 + 314 + 254 bp) (S1 Fig). Allele frequency distributions of EcoRI-RFLP polymorphism in thirteen pig populations are shown in Table 1. This polymorphic SNP was segregated in the tested pig populations, most of which had high frequencies of allele G at this site. The association analysis of the genotypes with the fat deposition and carcass traits were conducted in three pig populations including the Large White × Meishan F 2 pigs, France, and American Large White pigs. There were significant differences in the fat deposition and carcass traits among different genotypes in Large White × Meishan F 2 populations. Pigs carrying genotype GG had higher BFT, RIBBF, LFP and IFR, but lower CL than the AA or AG animals (Table 2). Furthermore, significant additive effects were also observed for these traits ( Table 2). The significant associations of g.1015 A>G with the corrected live backfat thickness were also found in France and American Large White pigs (Table 3). Taken together, we concluded that this functional SNP was significantly associated with fat deposition traits.

The SNPs could significantly affect the expression of PPARD gene in vivo
To estimate whether the SNP g.1015A>G could affect the gene expression in vivo, we selected France Large White pigs with genotypes AG and GG to compare the expression of PPARD gene in adipose and muscle tissues. The expression of PPARD gene between two genotypes was significantly different at both mRNA and protein level. Pigs with genotype GG had significant higher expression level than that of pigs with genotype AG (Fig 6).

Discussion
Functional mutations in the regulatory region have effects on the gene expression and phenotypic variation. PPARD is a promising candidate gene for pig carcass traits [27,33,42]. Currently, there is not report about any functional SNPs in the 5' regulatory region of the porcine PPARD gene. In this study, we identified the functional SNPs in the 5' regulatory region which directly affected the promoter transcriptional activities. The transcriptional activity of the PPARD promoter fragment carrying haplotype AC was the lowest, and haplotype AC had the lowest binding capacities with TCF7L2. TCF7L2 is a member of the T cell factor (TCF) family of transcription factors, and involved in the control of cell growth and signaling of winglesstype MMTV integration site family (Wnt) receptors [43]. TCF7L2 interacts with β-catenin in the nucleus and regulates the expression of Wnt target genes by acting as transcriptional activator or repressor [44][45][46][47]. TCF4-binding sites are located at long distances from transcription start sites, and defined as evolutionarily conserved A-C/G-A/T-T-C-A-A-A-G motifs [40]. In our study, three TCF7L2 binding sites were predicted in the promoter region between nucleotides −1636 bp and + 464 bp, and one of the three sites was located within the promoter sequence harboring the two SNPs. TCF7L2 can bind to the promoter of PPARD gene and alter gene expression in human CRC cells [41]. We also confirmed that the expression of PPARD gene was up-regulated by TCF7L2 overexpression in C2C12 cells. Besides, the porcine TCF7L2 gene improved the transcriptional activities of the PPARD promoter fragments  The SNPs in PPARD Gene Have Association with Fat Deposition Traits carrying the haplotypes AT, GT, and GC. However no significant effects on the haplotype AC were found, because there was no interaction between promoter fragments carrying haplotype AC and TCF7L2 protein.
PPARD gene has been implicated in lipid metabolism via several signaling pathways including Wnt/FZD and Wnt/β-catenin pathways. On the one hand, PPARD promote the accumulation of cholesterol and serum high density lipoprotein (HDL), and decrease the level of triglyceride in human macrophages by lipogenesis pathways such as Wnt/FZD and Wnt/βcatenin pathway [48,49]. On the other hand, it reduces body weight, lipid droplet number and size, and up-regulates gene expression related to fatty acid oxidation [34]. In human population, three SNPs in the PPARD affect (lifestyle intervention) LI-induced changes in overall adiposity hepatic fat storage [50]. However, the molecular mechanism by which genetic variants within upstream regulatory region of PPARD gene affect fat deposition traits still remain unclear. In this study, the allele A was mostly linked with allele C at two loci in all the pig breeds studied, and no haplotype AT was found in American Large White pigs, France Large White pigs and Large White × Meishan F2 pigs studied. We also found that the functional variant (g. 1015A>G) had significant association with fat deposition traits in three pig populations and allele A was associated with lower fat deposition traits as well as lower gene expression of PPARD in vitro and in vivo. Previous studies showed that PPARD and CDKN1A gene in preadipocytes had lower expression levels in heterozygote MSQTL7/LWQTL7 (the heterozygous  [51], which was in agreement with our results. Taken together, it can be inferred that haplotype AC inhibited the promoter transcriptional activity through decreasing the binding capacities of TCF7L2 to its binding sites, thereby repressed the PPARD expression and resulted in decreased fat deposition traits. Therefore, the expression change of PPARD gene may be an important cause of the major QTL effects for fat deposition.
In pig breeding schemes, the lean percentage in carcass can be indirectly selected by live backfat thickness because there is significant negative correlation between these two traits. In this study, there were significant differences in fat deposition traits between different genotypes. Pigs with genotypes AA or AG had significant lower fat deposition trait values compared with those of genotype GG in the pig populations studied. Considering the effects of allele A in this study, enhancing the allele frequency of allele A will lead to the reduction of backfat thickness in the pig populations of which the allele was not fixed at this locus. However, selection of allele A also has adverse effects on intramuscular fat and meat quality, as intramuscular fat and backfat deposition are negatively correlated.
In conclusion, our results provided a new evidence that the newly identified functional SNPs in the 5' regulatory region of PPARD genes have significant effects on the fat deposition and carcass traits. Further work will be necessary to confirm the effects of this functional SNP in more pig populations and investigate the molecular mechanism regarding the roles of SNPs in phenotype variation.