Genetic characterisation of PPARG, CEBPA and RXRA, and their influence on meat quality traits in cattle

Background Peroxisome proliferator-activated receptor gamma (PPARG), CCAAT/enhancer binding protein alpha (CEBPA) and retinoid X receptor alpha (RXRA) are nuclear transcription factors that play important roles in regulation of adipogenesis and fat deposition. The objectives of this study were to characterise the variability of these three candidate genes in a mixed sample panel composed of several cattle breeds with different meat quality, validate single nucleotide polymorphisms (SNPs) in a local crossbred population (Angus - Hereford - Limousin) and evaluate their effects on meat quality traits (backfat thickness, intramuscular fat content and fatty acid composition), supporting the association tests with bioinformatic predictive studies. Results Globally, nine SNPs were detected in the PPARG and CEBPA genes within our mixed panel, including a novel SNP in the latter. Three of these nine, along with seven other SNPs selected from the Single Nucleotide Polymorphism database (SNPdb), including SNPs in the RXRA gene, were validated in the crossbred population (N = 260). After validation, five of these SNPs were evaluated for genotype effects on fatty acid content and composition. Significant effects were observed on backfat thickness and different fatty acid contents (P < 0.05). Some of these SNPs caused slight differences in mRNA structure stability and/or putative binding sites for proteins. Conclusions PPARG and CEBPA showed low to moderate variability in our sample panel. Variations in these genes, along with RXRA, may explain part of the genetic variation in fat content and composition. Our results may contribute to knowledge about genetic variation in meat quality traits in cattle and should be evaluated in larger independent populations. Electronic supplementary material The online version of this article (doi:10.1186/s40781-016-0095-3) contains supplementary material, which is available to authorized users.


Background
Fat content and composition are considered major economically important traits in livestock, since variations in these two factors affect several meat properties [1]. These traits are the result of several biological processes, such as adipogenesis, lipolysis and fatty acid-transfer. Therefore, a part of the variability produced by these processes may be attributed to the genetic variants of the pathway members. Peroxisome proliferator-activated receptor gamma (PPARG), CCAAT/enhancer binding protein alpha (CEBPA) and retinoid X receptor alpha (RXRA) are important nuclear transcription factors involved in numerous cellular processes [2], and are considered key molecules in regulation of adipogenesis. PPARG and CEBPA are induced early in the signaling pathway, they work together to trigger the process and regulate each other [3]. PPARG acts as heterodimer with RXRA, which belongs to a family of nuclear receptors that act as homodimers and heterodimers. In view of their roles, the genetic control of adipogenesis by PPARG, CEBPA and RXRA may be important and helpful for animal improvement.
During the last years, SNPs in PPARG and CEBPA have been associated with a group of meat quality traits in Chinese and Korean cattle, including tenderness, backfat thickness, water holding capacity, fatty acid composition, weight at slaughter and marbling, among others [4][5][6][7][8][9]. However, those works have been performed almost entirely using Asian cattle under feedlot conditions and the results might not be necessarily comparable with researches performed with other breeds under pasture-based feeding. For instance, these two conditions may activate specific metabolic pathways governed by different genes. Nowadays, most of the exported beef in the world is produced on pasture-based systems [10], as in countries like Argentina, Brazil, New Zealand, Paraguay and Uruguay, among others.
In this context, we searched for gene variants in PPARG and CEBPA in a sample set composed of nine cattle breeds with different meat quality. Then, we validated some of these SNPs, along with SNPs in the RXRA gene, in a local Angus-Hereford-Limousin crossbred population (N = 260) fed on pasture-based conditions. We used this population to evaluate the association of the SNPs with intramuscular fat content (IF), backfat thickness (BT) and fatty acid composition. Finally, we analysed the molecular effects of these SNPs through bioinformatic predictive tools.
All animals grazed a sown pasture (predominantly Lolium multiflorum, Dactylis glomerata, Bromus catarthicus, Trifolium repens and Trifolium pratense). They were slaughtered in eight groups and meat blocks were taken from the 13th rib to perform association studies (Additional file 2: Table S1). The decision to sample this experimental population instead of other commercial cattle populations was based on the availability of reliable information in terms of phenotypic data, management and genetic background of the animals.
DNA was isolated from blood lymphocytes using Wizard® Genomic DNA purification kit (Promega, Madison, WI, USA) following the instructions of the supplier, and from meat samples as previously described in [11].
Re-sequencing study of the bovine PPARG and CEBPA genes To amplify the coding regions of the PPARG and CEBPA genes, ten pairs of primers were designed according to the DNA sequences available in GenBank [Gene IDs: 281677, 281993] (Additional file 2: Table S2). PCR reactions were performed, verified, purified and sequenced as described in [12]. Sequences were aligned using CLUSTAL-X 2.1 [13] and variants were defined by direct comparison with the bovine reference sequences.

SNP selection and genotyping
As the crossbred population used to validate SNPs and perform association tests included animals from Angus, Hereford and Limousin (pure or crossbred), only some of the SNPs detected in the re-sequencing stage were considered for further validation. In other words, many of the SNPs detected by re-sequencing showed no variation in the Taurine breeds and were not considered, since they would probably show no variation in the crossbred population. For this reason, additional SNPs were selected from dbSNP [14] to have a better covering of the length of the genes using markers with proven variation in Taurine breeds. At this stage, the addition of another candidate gene was decided, and SNPs in the RXRA gene were also selected from dbSNP to be validated and tested for associations. Genotyping was performed by the Neogen genotyping service (USA) using the Sequenom platform [15].

Meat quality measurement
Fatty acid content and composition measurements were gathered from blocks of meat obtained from the 260 animals of the crossbred population to perform association studies. These blocks, corresponding to the Longissimus dorsi muscle (13th rib) were extracted from the carcass 24 h after slaughter. Backfat Thickness (BT) was measured perpendicular to the outer surface at a point three quarters of the length of the Longissimus dorsi muscle from the end of the loin bone, and expressed in millimetres. The Intramuscular Fat (IF) and fatty acid composition were then measured as described in [12]. The measured fatty acids were: myristic acid (C14:0); myristoleic acid (C14:1); palmitic acid (C16:0); palmitoleic acid (C16:1); stearic acid (C18:0); oleic acid (C18:1 cis-9); linoleic acid (C18:2 cis-9,12); γ-linolenic acid (C18:3 cis-6,9,12); αlinolenic acid (C18:3 cis-9,12,15); total saturated fatty acids (SFA); total monounsaturated fatty acids (MUFA); and proportion between omega-6 and omega-3 fatty acids (Ω6/Ω3). The fatty acid contents were expressed as percentage of total fatty acids per sample. C20:0 and other long-chain fatty acids were not included in the analysis since their percentages were lower than 0.5 %. The means, standard deviations, minimum and maximum values of all these measurements were included in Additional file 3: Table S3.

Statistical Analysis and association with meat quality
Haplotypes and linkage disequilibrium (LD) among SNPs were estimated and visualized on HAPLOVIEW v4.2 [16] using the four gamete rule and the solid spine of LD. Allele frequencies and Hardy-Weinberg equilibrium (HWE) were analysed using GENEPOP software [17]. The 95 % confidence intervals for allele frequencies were computed using the binomial distribution implemented in R with the binom.confint function (http:// cran.r-project.org/web/packages/binom/). Values for unbiased expected (h e ) and observed (h o ) heterozygosity were calculated using ARLEQUIN v3.5 [18].
The association of genotypes with BT, IF and fatty acid composition was evaluated using mixed models. BT and IF were analysed according to the following model: Where Y ijkl = observed value of the phenotypic variable, μ = intercept, GG j = fixed effect of the j th genetic group, SG k = fixed effect of the k th slaughter group, a = additive effect for the SNP, Z1 i is the incidence variable for the additive effect (that is 0 for one of the homozygous genotypes, 1 for the heterozygous genotype and 2 for the alternative homozygous one), d = dominant effect for the SNP, Z2 i is the incidence variable for the dominance effect (that is 0 for both of the homozygous genotypes and 1 for the heterozygous genotype), βW i = animal weight covariate for the i th animal, S l = random effect of the l th sire, e ijkl = random error. The same single trait model was used for fatty acid composition variables, but using ether extract instead of animal weight as covariate.
All statistical analyses were performed using the MIXED procedure of SAS software [19]. When the additive or dominance effects of the SNP were statistically significant (P < 0.05), the substitution effect (α) was calculated considering the frequencies of the major and minor alleles (p and q, respectively) and using the following equation [20]: The variance explained by the SNP (σ SNP 2 ) was also estimated for each SNP-trait test as follows: where RMS is the residual of the reduced model (SNP effect excluded), and FMS is the residual of the full model (SNP effect included). After all trait-SNP tests were performed, the false discovery rate (FDR) for multiple comparisons was controlled through the Benjamini & Hochberg method [21].

Bioinformatic analyses
The SNPs were also analysed through different bioinformatic prediction tools. For synonymous SNPs, changes in codon frequency usage were analysed through the Codon Usage Database [22]. For SNPs located in 5' UTR regions, the complete 5' UTR fragments of the RNA sequences were run on the Mfold Web Server [23] to compare stability among variants. These fragments were also analysed using RBPDB [24], the database of RNAbinding protein (RBP) specificities, considering a threshold of 0.8 to identify putative RBP binding sites. Finally, variations located in promoter regions were analysed through PhysBinder [25], considering all human and murine models available and the "average" threshold to predict putative transcription factor binding sites.

Results and discussion
Re-sequencing study A total of seven SNPs were identified in the PPARG gene using the panel of nine breeds. All of them had been previously reported and most showed very low frequencies. In fact, no homozygous genotypes were detected for any of the alternative alleles. Three of the SNPs were located in UTR regions: rs207671117 (5' UTR), rs211388309 (3' UTR) and rs207724742 (3' UTR); the first in Angus and Hereford, and the other two in Brahman and Nellore. The remaining four SNPs (rs207739706, rs41610552, rs110194439 and rs42661651) were detected in non-coding regions in different breeds ( Table 1).
The haplotype and LD analysis, estimated with the four gamete rule, showed two blocks: a small one (3' end), composed by two completely linked SNPs (rs42661651 and rs110194439), and a big one that included the other five SNPs and consisted of five haplotypes, with three of them in very low frequencies (Fig. 1). The same study, but estimated with solid spine of LD, showed one big block constituted by all the SNPs, with three haplotypes in very low frequencies.
Only two SNPs were detected in the CEBPA gene. These SNPs were located in the coding region of the gene, which comprised only one exon, and one of them had no previous reports in dbSNP. This novel SNP (ss1751108604) was detected in both Zebuine breeds (Brahman and Nellore) and the Japanese breed Wagyu. This SNP caused an amino acid change involving two neutral and polar residues (Ser139Asn). The other one was a synonymous SNP (rs110793792), widely distributed among the breeds, with the exception of Nellore and Shorthorn, which showed different homozygous genotypes for the mutation and no variability within the samples (Table 1).
Nine mutations were detected in total in this first stage. In other terms, we detected one SNP every 423 bp over 3809 bp analysed. Considering subspecies distribution,  [12]. This is consistent with the roles these genes play in lipid metabolism, since PPARG and CEBPA are key regulators in the first stages of fat deposition, among other processes, and LIPE codifies for an enzyme with very specific functions in lipid hydrolysis. In this scenario, LIPE may be subject to a less selective pressure than PPARG and CEBPA. In the case of PPARG, mutations were generally detected in low frequencies, which generated large linkage blocks with few haplotypes carrying most of the variation.

Allelic and genotypic frequencies
A group of SNPs was selected either from the results of the re-sequencing study or dbSNP to perform validation.
In the case of PPARG, two SNPs were selected from the re-sequencing study (rs207671117 and rs41610552) and two other SNPs were selected from dbSNP. These two SNPs from the dbSNP had been previously associated with a group of meat quality traits: rs42016945, located upstream of PPARG-2 and also part of the 5' UTR of splice-variant 1 (PPARG-1), and rs109613657, which caused an amino acid change in the seventh exon (Glu448His). For CEBPA, we selected rs110793792 from the re-sequencing stage, since the novel SNP (ss1751108604) showed no variability in the Taurine breeds despite its novelty, and rs210446561 from the dbSNP. We also selected four SNPs from the RXRA gene. Since this gene was included in the study after the re-sequencing stage and there were no previous reports of associations with meat quality traits to our knowledge, the SNPs were chosen directly from dbSNP. These SNPs were rs209839910 (Pro108Ser), located in the second exon, rs136289117 (synonymous), located in the ninth exon, and rs133517803 and rs207774429, which may be located in the first intron or a putative promoter region for the splice-  Table 2). The HWE test showed no significant deviations from the theoretical proportions, with the exception of rs42016945 in subpopulation LX (P < 0.05) ( Table 3). The unbiased expected heterozygosity (h e ) of the two balanced SNPs (rs41610552 and rs42016945) varied between 0.22 (¾ H) and 0.49 (LX). Observed heterozygosity (h o ) varied between 0.25 (¾ H) and 0.55 (LX). When linkage disequilibrium was analysed, rs207671117 and rs42016945 showed a small block with three haplotypes, two of them with more than 95 % of the haplotype frequencies (Additional file 4: Figure S1). SNP rs110793792, located in CEBPA, could not be genotyped by the Sequenom platform, reason why it was discarded from the analysis. The other SNP from this gene, rs210446561, was genotyped efficiently in only 143 samples. The analysis showed higher frequencies for allele C, with MAF ≥ 0.09 in the subpopulations ( Table 2). The HWE test showed no deviations from the theoretical proportions. Values for h e ranged from 0.16 (LX) to 0.38 (A), and h o ranged from 0.17 (LX) to 0.50 (A) ( Table 3).
Two of the SNPs from RXRA, rs207774429 and rs133517803, showed relatively balanced allele frequencies (MAF ≥ 0.06). It is worth mentioning that rs207774429, as happened with rs210446561 (CEBPA), was genotyped efficiently in only 153 samples. The other two SNPs, rs136289117 and rs209839910, showed very low variation and were not considered for association (Table 2). According to the HWE test, significant deviations were observed for rs207774429 in the whole population (P < 0.01) and rs133517803 in subpopulation ¾A (P < 0.01). The remaining SNPs showed no significant deviations from the theoretical proportions. Three of the SNPs in this gene were part of a linkage block constituted by three haplotypes, with two of them accounting for over 97 % of the haplotype frequencies (Additional file 4: Figure S1).
As we already mentioned, some of the SNPs selected for validation were not efficiently genotyped by Sequenom, but the rest showed different allele frequencies among the genetic groups. In general, the highest or lowest frequencies were observed in the LX subpopulation, which is historically and productively the most different breed, since Limousin is a European continental breed, and Angus and Hereford are Scottish. Deviations from the HWE proportions were observed mainly in crossbreeds, as expected. These deviations would also suggest the violation of some of Hardy-Weinberg assumptions and the existence of phenomena like selective mating, small population size, endogamy, and especially migration, considering the original purpose of the population, i.e., the evaluation of crossbreeding systems.

Association study
When individual tests were performed, several genotype effects were observed on the evaluated traits. The least square means for each genotype, percentages of phenotypic variance explained by the SNPs, and substitution effects were estimated for all the evaluated traits and presented in Table 4. Regarding SNPs in PPARG, rs207671117 showed a significant effect on Ω6/Ω3 (P < 0.05), but small genetic variation in general. An additive effect on BT was detected (P < 0.05) for rs42016945, while dominance effects were significant (P < 0.05) for rs41610552 on C18:1 cis-9 and rs42016945 on C18:0. For these SNPs, the percentages of phenotypic variance explained ranged from 0.31 to 2.09 %. SNP rs133517803 (RXRA) showed significant effects on several measures: additive effects were detected on C18:1 cis-9; C18:3 cis-6,9,12; C18:3 cis-9,12,15; MUFA and BT; while a dominance effect was detected on C18:1 cis-9; C18:3 cis-9,12,15 and MUFA. For these traits, the SNP explained from 0.35 to 1.77 % of the phenotypic variance (Table 4). No significant effects were observed on the other traits. It is worth mentioning that none of these effects reached the threshold for statistical significance in the FDR control by means of the Benjamini-Hochberg method considering 70 tests.
The explained percentages of variance were interestingly high for oleic acid and MUFA compared with other traits such as BT and stearic acid. In particular, the percentage of phenotypic variance explained by SNP rs133517803 in RXRA for oleic acid, and subsequently for MUFA, was high, which may suggest an important role of this gene in the oleic acid deposition in muscle.
If we consider the first approach, our results were partially consistent with previous reports from other authors. Sevane et al. [9] reported associations for SNP rs42016945 with several Ω-3 fatty acids. Here, rs42016945 showed a significant effect on fatty acid composition, but on C18:0 instead. We also found a possible effect on BT, which was interesting and expectable given the nature of PPARG. None of the mutations in PPARG reported in Chinese and Korean breeds [5,7] was detected in our panel, probably due to different breed origins, although it is important to mention the finding of significant effects on BT in Chinese cattle, which is concordant with our findings in the European breeds. SNP rs110793792 (CEBPA) had been associated with BT and marbling in Chinese breeds [8], but these effects could not be tested in our population due to a genotyping problem. To our knowledge, there are no previous works reporting associations for RXRA, so this work may provide evidence that SNPs in this gene may be helpful for animal improvement by means of markerassisted selection programs.

Bioinformatic predictions
Since SNP rs207671117 was located in the 5' UTR of PPARG, we analysed the possible effects on mRNA stability and putative RBP binding sites. Two highly similar structures were obtained running the UTR sequences of the alternative variants in The Mfold Web Server, but the structure for allele A (ΔG = -42.60 kcal/mol) seemed slightly more stable than that for allele G (ΔG = -42.20 kcal/mol). When RBP binding sites were analysed (threshold = 0.8), we found that this SNP was located one base away from a putative binding site for protein FUS (SCORE = 7.36), which is important in maintaining genomic integrity. These same studies were performed for rs42016945 and we found that the structures provided by Mfold were quite different despite the similarity of energy values. The structure for allele G (ΔG = -27.70 kcal/mol) seemed more stable than structure for allele A (ΔG = -26.90 kcal/mol). According to the analysis on RBPDB, no sites for RBP were identified promptly at the mutated site, but the SNP was immediately next to a NONO (non-POU domain-containing octamer-binding protein) binding site (SCORE = 8.95). NONO is a protein involved in numerous nuclear processes like unwinding, recombination, DNA binding and regulation of splicing. SNPs rs133517803 (RXRA) was located in a possible alternative promoter region. Therefore, we searched for putative transcription factor binding sites through Phys-Binder. According to this tool, rs133517803 was located over an ESRRB (estrogen-related receptor beta) binding site (threshold = 308), whose role is still not clear.

Conclusions
PPARG and CEBPA showed low to moderate variability in our mixed sample panel. Variations in these genes, along with RXRA, may explain part of the phenotypic variation in fat content and composition of meat, especially SNPs in RXRA, which explained an important part of the variation in the highly heritable oleic acid percentage and MUFA. The molecular bases of the phenotypic differences may be partially explained by changes in RNA structures, RBP binding sites, codon usage frequencies and TF binding sites. The SNPs we analysed should be evaluated in independent populations with in-vitro and in-vivo analyses to explain the mechanisms by which these polymorphisms may be involved in the traits.