Variants of the bovine retinoic acid receptor-related orphan receptor C gene are in linkage disequilibrium with QTL for milk production traits on chromosome 3 in a beef × dairy crossbreed population

The bovine retinoic acid receptor-related orphan receptor C gene (RORC) is located on bovine chromosome 3 in the vicinity of several QTL for milk production traits, and an association of RORC genetic variants with carcass fatness traits had been confirmed by several studies. In a F2 resource population, a chromosome-wise QTL and association study with 26 genetic markers (microsatellites and SNPs from genes involved in fat metabolism) on BTA3 was performed. In the analysis, 183 first and 152 second lactation records of 183 cows were included. The QTL study revealed five QTL affecting milk yield (MY), fat percentage (FP) and protein percentage (PP) as well as the milk fat/milk protein ratio (FPR), respectively, in the chromosomal region harbouring the RORC gene. The association study displayed a significant association between RORC c.1138+65A>G and MY and also associations between RORC c.934-262T>G and MY, FP and FPR, respectively. A combined QTL association study showed that these SNPs included as fixed effect in the corresponding model resulted in a prominent drop of the QTL test statistic. For the RORC c.934-262T allele, which had been confirmed to increase intramuscular fat deposition, an increasing effect on milk fat content traits was detected. The RORC c.1138+65G allele that had been shown to positively affect rump fat thickness, was associated with increased MY in our study. In conclusion, we found indications that the RORC SNPs, which previously had been highlighted due to their effects on carcass fat deposition, are also associated with milk production traits.


Introduction
Most traits of economic interest in livestock are determined by a combination of environmental factors and a number of unknown genes.Linkage analysis and association studies enabled the dissection of the unknown genetic factors into quantitative trait loci (QTL), which vary in number depending on the specific target trait (Hayes et al. 2010).
In addition to the chromosomes harbouring major milk production QTL (BTA6, BTA14, BTA20), there were reports from several dairy breeds on a QTL for milk fat production on BTA3 (Heyen et al. 1999, Viitala et al. 2003).In addition, there were confirmed results regarding a QTL affecting body fat deposition in the same chromosomal region in beef cattle populations (Casas et al. 2003).This QTL was demonstrated to be associated with genetic variants within the Retinoic acid receptor-related orphan receptor C (RORC) gene.RORC, also known as RORγ, belongs to the orphan nuclear receptor superfamily.Meissburger et al. (2011) found that RORC affects the size of adipocytes, modulates insulin sensitivity in obesity, and is involved in the control of adipogenesis.When Barendse et al. (2007) analysed genetic variation in the RORC gene, the strongest association was found between the RORC c.934-262T>G mutation and marbling score, whereas the RORC c.1138+65A>G mutation was highly associated with rump fat in Angus and Shorthorn cattle.RORC haplotypes were shown to explain most of the QTL effect on body fat deposition on BTA3 (Barendse et al. 2007).The RORC c.934-262T allele that increased marbling score was also the favourable allele that increased marbling score and intramuscular fat content in a confirmation study in other cattle breeds (Barendse et al. 2010).
While the RORC mutations have been extensively investigated for their effects on body fat deposition, there are no reports regarding corresponding effects on milk fat production in beef cattle populations.To address this issue, we took advantage of a crossbreed cow population of a dairy (German Holstein) and a beef breed (Charolais), which was raised and kept under dairy production conditions (Kühn et al. 2002), but predominantly showed a beef phenotype (Hammon et al. 2010).In this population, Hammon et al. (2010) and Hamada et al. (2012) had observed a negative correlation between milk yield (MY) on one hand and body fat deposition and carcass weight on the other hand.

Animals and data collection
Our study comprised 183 cows of a F 2 resource population established from a cross between Charolais and German Holstein (SEGFAM resource population, Kühn et al. 2002), for which 183 first full lactation records and 152 second lactation records until slaughter at day 30 of lactation were available.All cows finishing the 30th day of the second lactation were slaughtered for experimental reasons.All animals were housed in a loose stall barn under identical conditions at the Leibniz Institute for Farm Animal Biology (FBN) in Dummerstorf, Germany and were milked twice a day in a commercial tandem milking parlour.The total mixed ration (TMR) diet fed ad libitum had an energy content of 6.6 MJ of NEL/kg.The diet was supplemented by 1 kg concentrates per 2 kg milk, when milk yield exceeded 17 kg/day (Hammon et al. 2010).
Milk yield was recorded daily during the normal milking process.Milk fat, milk protein and lactose content was measured weekly using accredited procedures with the MilkoScan FT+/Fossomatic FC-Combination (Foss GmbH, Rellingen, Germany) at the laboratory of the Landeskontrollverband für Leistungs-und Qualitätsprüfung Mecklenburg-Vorpommern, in Güstrow, Germany.The SEGFAM cow population showed a production level found in beef type cattle and a strong deviation from the typical lactation curve of dairy cows, because milk yield strongly decreased in the interval day 30 to day 100 (Hammon et al. 2007).Thus, for our study, the lactation period was dissected into specific intervals (Table 1) for the yield traits (milk volume: MY, fat yield (FY), protein yield (PY) and lactose yield (LY).In addition, our analysis specifically considered milk fat (FP) and milk protein content (PP) as well as the milk fat to milk protein ratio (FPR) as indicator of the energy balance of the individuals (Buttchereit et al. 2011).

DNA extraction and genotyping
DNA was isolated from blood using the NucleoSpin Blood L kit (Macherey-Nagel GmbH & Co. KG, Düren, Germany) according to the manufacturer's instructions.
The initial QTL study included seven microsatellite markers distributed across BTA3 and 19 SNPs in genes involved in fat metabolism and located in the region of previously described QTL on BTA 3 (Table 2).
Genotyping of the microsatellite markers was performed according to standard procedures as described by Kühn & Weikard (2007) on a MegaBace capillary sequenzer (GE Healthcare, Little Chalfont, UK).The genotyping of the investigated SNPs was performed on an Illumina Bead station (Illumina Inc., San Diego, CA, USA) (Oliphant et al. 2002) as part of a targeted 384-SNP GoldenGate assay (Eberlein et al. 2009).

QTL and association analyses
For recalculating a genetic map of BTA3, CRI-MAP v2.50 (Green et al. 1990) was used, incorporating modifications by Ian Evans and Jill Maddox.If marker pairs were sharing identical positions on the resulting genetic map, a distance of 0.01 cM was assumed for providing a minimal distance, which is required for QTL linkage mapping (Table 2).The calculated genetic map was utilised for a variance component QTL model for initial linkage analysis using Qxpak v5.02 (Perez-Enciso & Misztal 2004): where y is a vector of phenotypes, b represents a vector of fixed effects (age at calving, calving year, calving season), u symbolises the vector of the individual infinitesimal polygenic effects, g is a vector of additive QTL effects; F, Z and Q symbolise the incidence matrices for the fixed, polygenic and the additive QTL effect, respectively and e represents the vector of random residuals.An association linkage disequilibrium (LD) analysis was performed to evaluate whether the used SNPs were associated with additive effects on milk performance and composition traits by using Qxpak v5.02 (Pérez-Enciso & Misztal 2004).The following model was applied: where y i is the record of individual i, f e represents the fixed effects (age at calving, calving year, calving season), λ i is an indicator variable that depends on the individual genotype at the analysed loci and is defined as 1, 0, or −1, δ k is the additive allelic effect of allele k, u i is the infinitesimal polygenetic effect of individual i, and e dik is the residual.A chromosome-wise Bonferroni correction was made to avoid false positive associations.SNPs significant in the association analysis were included in a combined QTL model as a fixed effect to test their impact on the QTL test statistic essentially as described by Eberlein et al. (2009).
The statistical significance of the QTL analyses was tested using a likelihood-ratio test (LRT).The variance in the model explained by the QTL was calculated as the relative reduction of the residual variance due to including the QTL in the model (Knott et al. 1996).

QTL analyses
The initial QTL analysis identified five significant QTL with a consistent peak of the test statistic in the proximal region of BTA 3: four QTL for first lactation traits [MY from lactation day (d) 1 to 5 (P<0.01) and from d 6 to 30 (P<0.1),FP on d 30 (P<0.01) and FPR on d 30 (P<0.05)] and one QTL for a second lactation trait [PP on d 14 (P<0.01)](Figures 1-3).

Association study
Ten of the 19 SNPs showed significant associations with milk performance and composition traits after Bonferroni correction (Table 3).For several traits, more than one SNP showed a significant association, most prominently for MY from d 1 to d 5.For this trait, five SNPs in three genes showed a significant association for first or second lactation data: RORC c.934-262T>G, RORC c.1138+65A>G, ANXA9 c.251T>C, PHGDH c.*116A>G and PHGDH c.945+1014T>C.SNP effects were calculated as difference between alternative alleles at a SNP locus.The results showed that the allele RORC c.1138+65G is associated with an increased MY from d 1 to d 5 and from d 6 to d 30, while the RORC c.934-262T allele is associated with higher FP and a higher FPR on d 30 (Table 3).Aligning the data from the association study to the results from the linkage analysis, significant associations were revealed between RORC c.934-262T>G and FP and FPR on d 30 corresponding to the QTL data (Table 3 and Figure 2).

Combined QTL and association study
The SNPs with significant associations to MY from d 1 to d 5 and FP as well as FPR on d 30 were used as fixed effects in a QTL study to analyse their relevance for the detected QTL.Furthermore, nominally significant SNP associations for the traits MY from d 6 to 30 (first lactation) and PP on d 14 (second lactation) (data not shown) were also included in the corresponding QTL model.The results showed that RORC c.934-262T>G and RORC c.1138+65A>G, respectively, had the most prominent effect on the test statistic of the identified QTL.Although RORC c.1138+65A>G is located 1 cM adjacent to the 96.8 % QTL confidence interval for MY from d 1 to d 5, it had a large effect on the QTL test statistic (Figure 1).Effects analogous to this SNP were observed for the QTL affecting MY from d 6 to d 30.RORC c.934-262T>G had a prominent effect on the QTL test statistic for FP and FPR on d 30.
When including this SNP in the QTL model, the QTL test statistic for FPR on d 30 showed an almost complete drop (Figure 2).The SNP RORC c.934-262T>G had also a large effect on the test statistic for the QTL affecting PP on the d 14 in the second lactation (Figure 3).In all combined QTL and association analyses, which included the RORC SNPs as a fixed effect, the peaks of the QTL test statistics shifted towards the telomere of the chromosome compared to the QTL model without SNP effect included.

Discussion
Previous analyses within our resource population had indicated that higher body fat accretion was accompanied with a lower MY (Hammon et al. 2010), which suggested that genes modulating body fat deposition would also have an effect on milk production traits.We focused our study on BTA3 due to several QTL for milk production and milk composition traits reported for this chromosome (Cattle QTL database [http://www.animalgenome.org/cgi-bin/QTLdb/BT/index], Hu et al. 2007; Combined QTL Map of Dairy Cattle Traits [http:// firefly.vetsci.usyd.edu.au/reprogen/QTL_Map],Khatkar et al. 2004), and because the RORC gene that has been confirmed as a candidate gene for body fat deposition (Barendse et al. 2007(Barendse et al. , 2010) ) in several cattle populations is located in the vicinity of those QTL.Our data represent the QTL for milk production traits in a beef type cattle population.The results support the hypothesis that the RORC gene known for its association with body fat deposition may also be associated with milk production traits including milk fat percentage at the beginning of the lactation.We confirmed significant QTL for milk production traits on BTA3 and showed that the RORC variants explained a substantial proportion of the QTL variance as demonstrated by the drop of the QTL test statistic when adding the respective SNPs into the model.Three of the five detected QTL (MY from d 1 to 5 as well as d 6 to 30 in the first lactation and PP in the second lactation) had their highest test statistic peak in the region between INRA006 and BMS2522, where Viitala et al. (2003) found a QTL for FP and PP and Ashwell et al. (2004) for FY and PP.Furthermore, Kolbehdari et al. (2009) found an association between MY and LOC511740 (rs43709850), which is 1.2 Mb telomeric to RORC, and an association between the trait PP and the PDZK1 gene (rs41587408) (2.7 Mb telomeric to RORC).Those results are in positional disagreement with a recently published study combining fine mapping and resequencing for the detection of QTL for milk production traits on BTA3 in a Holstein dairy cattle family design by Cohen-Zinder et al. (2011).They highlighted four candidate genes (RAP1A, ADORA3, OVGP1 and C3H1orf88) in the region 27.1-36.2Mb between the microsatellites BL41 and TGLA263 as potentially underlying QTLs found for MY, FY and PY.However, it has to be considered, that the QTL identified in our study focussed on the very beginning of the lactation and might have another genetic background than QTL for phenotypes covering the entire lactation period.
The analysis of the allelic effects in this study revealed that the RORC c.934-262T allele is associated with a higher FP and a higher FPR on d 30 and a decreased MY from d 1 to d 5 in the first lactation.A similar significant effect on MY was observed for the RORC c.1138+65A allele, which was associated with decreased MY from d 1 to d 5 and d 6 to d 30.Comparing our results for milk production traits with the results described for RORC effects on carcass traits (Barendse et al. 2007(Barendse et al. , 2010)), it seemed that the wildtype allele RORC c.934-262T, which increased marbling score and intramuscular fat content but decreased rump fat thickness, may have positive effects on milk fat percentage.The RORC c.1138+65A allele related to decreasing marbling score and increasing rump fat thickness is mainly associated with effects on MY in early lactation.In summary, our study showed for the first time that genetic variants in a gene (RORC) that initially had been associated with fatness in cattle are also in linkage disequilibrium with loci affecting milk production and composition traits.The regulatory functions of retinoidrelated orphan receptor γ (RORγ) encoded by the RORC gene, for modulating lactation are not yet fully understood.While the RORγ effects during adipogenesis are well described (Meissburger et al. 2011), it remains unclear whether the putative RORC effects detected in our study are direct results of RORγ modulating milk synthesis, or if the effects are indirect consequences of a primary effect on body fat deposition or if they are due to another, yet undetected mutation in the vicinity of the RORC gene.Barendse et al. (2010) concluded from their data that neither of the detected RORC variants was causal for the effects on body fat deposition.In this study, an association was found between RORC alleles and MY in the first days of lactation.At this time, typical dairy cows are in a negative energy balance due to the fact that the amount of energy via feed intake is not sufficient enough to cope with the increasing energy demand.Hammon et al. (2010) showed remarkable differences in energy partitioning between cows towards milk production or body fat deposition.Thus, it will have to be established in further confirmation studies whether variation in the RORC gene is involved in modulating this process at the beginning of lactation.

Figure 1
Figure 1 QTL test statistic for milk yield (MY) on lactation day 1 to 5 (•) as well as 6 to 30 (■) in the first lactation.RORC c.1138+65A>G is located at position 13.41 cM.Solid symbols: RORC c.1138+65A>G effect not included in the QTL model, open symbols: RORC c.1138+65A>G effect included in the QTL model as fixed effect.

Figure 2
Figure 2 QTL for fat percentage (FP •) and fat protein ratio (FPR ■) on lactation day 30 in the first lactation.RORC c.934-262T>G is located at position 13.40 cM.Solid symbols: RORC c.934-262T>G effect not included in the QTL model, open symbols: RORC c.934-262T>G effect included in the QTL model as fixed effect.

Figure 3
Figure 3 QTL for protein percentage (PP) on lactation day 14 in the second lactation.RORC c.934-262T>G is located at position 13.40 cM.Solid symbols: RORC c.934-262T>G effect not included in the QTL model, open symbol: RORC c.934-262T>G effect included in the QTL model as fixed effect.

Table 2
Overview of the investigated loci and their genomic position on BTA3 a genomic marker position on the reference sequence Bos_taurus_UMD_3.1 (AC_000160.1),b Genetic marker position in the linkage map calculated by CRI-MAP v2.50.