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

Two groups of samples were collected: the first group comprised blood samples from 43 unrelated purebred animals (Angus, n = 5; Brahman, n = 5; Creole, n = 5; Hereford, n = 5; Holstein, n = 5; Limousin, n = 4; Nellore, n = 4; Shorthorn, n = 5; Wagyu, n = 5), which were used to identify polymorphisms in the bovine PPARG and CEBPA genes. The second group comprised 260 steers (15–29 month-old), born between 2006 and 2010, which were used to perform population analyses, including SNP validation and association tests. This group of animals had been used in previous studies to evaluate crossbreeding systems under pasture grazing with strategic supplementation at the Experimental Station of the National Institute of Agricultural Technology (INTA, Balcarce, Argentina; Coordinates: 37°49S 58°15 W). Steers included: purebred Angus -A- (n = 44) and Hereford -H- (n = 26) steers, their crossbreeds F1 and F2 -½AH- (n = 95), reciprocal backcrosses -¾A and ¾ H- (n = 54), and steers produced by mating Limousin -L- sires with F1 crossbred cows -LX- (n = 41) (Additional file 1: Table S1). Fifty-four sires were used, including 17 A sires (1-16 steers), 18 H sires (1-11 steers), 8 AH sires (1-7 steers), 8 HA sires (1-7 steers) and 4 L sires (1-34 steers). L sires were mated only with AH and HA cows and, every year, some of the A and H sires were mated with more than one genetic group.

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].

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.

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].

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.

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²) 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].

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 RNA-binding 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.

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.

Fig. 1. Haplotypes (upper part) and linkage disequilibrium (lower part) among SNPs in the PPARG gene estimated in a mixed sample panel (n = 43). Blocks were estimated using the four gamete rule (a) and solid spine of LD (b). In both cases, r2 values are indicated inside the boxes and blocks are indicated in thick lines

Download Original Figure

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, our variation values translate to one SNP every 762 bp for Bos taurus and one SNP every 544 bp for B. indicus. As expected, variability was higher in the Zebuine group. A few years ago, the Bovine HapMap Consortium [26] obtained one SNP every 714 bp for Angus or Holstein, and one SNP every 285 bp for Brahman. Therefore, the variability obtained in this work was similar in the case of Taurine breeds but lower for Zebuine breeds. On the other hand, the variability observed here was lower than that reported lately by our group for the LIPE gene in this same sample panel, where a SNP was detected every 123 bp [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.

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-variant 3.

Regarding the SNPs in PPARG, rs207671117 showed low variability in subpopulations A, ¾ H, ½AH and LX (Minimum Allele Frequency [MAF] ≥ 0.02), and showed no variability in H and ¾A. On the other hand, SNPs rs41610552 and rs42016945 showed moderate allele frequencies among the subpopulations (MAF ≥ 0.12). Surprisingly, SNP rs109613657 showed no variability at all (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.

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 marker-assisted selection programs.

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 PhysBinder. According to this tool, rs133517803 was located over an ESRRB (estrogen-related receptor beta) binding site (threshold = 308), whose role is still not clear.