Adrenal cortex expression quantitative trait loci in a German Holstein × Charolais cross.

BACKGROUND: The importance of the adrenal gland in regard to lactation and reproduction in cattle has been recognized early. Caused by interest in animal welfare and the impact of stress on economically important traits in farm animals the adrenal gland and its function within the stress response is of increasing interest. However, the molecular mechanisms and pathways involved in stress-related effects on economically important traits in farm animals are not fully understood. Gene expression is an important mechanism underlying complex traits, and genetic variants affecting the transcript abundance are thought to influence the manifestation of an expressed phenotype. We therefore investigated the genetic background of adrenocortical gene expression by applying an adaptive linear rank test to identify genome-wide expression quantitative trait loci (eQTL) for adrenal cortex transcripts in cattle.
RESULTS: A total of 10,986 adrenal cortex transcripts and 37,204 single nucleotide polymorphisms (SNPs) were analysed in 145 F2 cows of a Charolais × German Holstein cross. We identified 505 SNPs that were associated with the abundance of 129 transcripts, comprising 482 cis effects and 17 trans effects. These SNPs were located on all chromosomes but X, 16, 24 and 28. Associated genes are mainly involved in molecular and cellular functions comprising free radical scavenging, cellular compromise, cell morphology and lipid metabolism, including genes such as CYP27A1 and LHCGR that have been shown to affect economically important traits in cattle.
CONCLUSIONS: In this study we showed that adrenocortical eQTL affect the expression of genes known to contribute to the phenotypic manifestation in cattle. Furthermore, some of the identified genes and related molecular pathways were previously shown to contribute to the phenotypic variation of behaviour, temperament and growth at the onset of puberty in the same population investigated here. We conclude that eQTL analysis appears to be a useful approach providing insight into the molecular and genetic background of complex traits in cattle and will help to understand molecular networks involved.

Background

eQTL analyses and the integration of proteomic and metabolomics data in genetic analyses contributed to the understanding of genes and molecular networks of complex traits in humans [1, 2] and farm animals [3-5]. In cattle most studies investigating the genetic background of complex traits such as milk performance or reproduction used quantitative trait loci (QTL) mapping and genome-wide association studies. Recent advances of high throughput technologies promoted the utilization of transcriptomic, proteomic or metabolomics data in genetic analyses [6].

Caused by interest in animal welfare and the impact of stress on economically important traits in farm animals the adrenal gland and its function within the stress response is of increasing interest. Stress-related effects on productivity [7-9], health [10, 11] and reproduction [12, 13] in cattle are not fully understood. However, insight into molecular networks affecting these traits in cattle could contribute to the understanding of the mutual relationship between stress and obesity or stress related effects on health in humans.

Adrenal cortex is a promising target to investigate these effects, due to its impact on stress response and regulation of essential metabolic pathways [14]. Human diseases such as Addison’s disease, the Cushing’s syndrome or Conn’s syndrome that are related to aberrations in glucocorticoid or mineralocorticoid levels further accentuate its relevance [15-17]. In cattle first observations indicating the importance of the adrenal gland for economically important traits were related to reproduction and lactation [18, 19]. Recent studies in cattle and pig showed that genetic variants in the glucocorticoid receptor and mineralocorticoid receptors are associated with cortisol concentrations and meat quality traits, for example [20, 21].

Therefore we investigated the genetic background of adrenocortical gene expression in a F2 cattle population (SEGFAM). The SEGFAM population is a cross of Holstein cows that were selected for high milk yield and Charolais bulls that were selected for meat production for a long time. By crossing the two breeds a high genetic and phenotypic variance in regard to milk and meat production traits was expected due to segregation of different QTL alleles in the F2 cross [22]. The population was established to investigate the genetic and physiological background of nutrient transformation types, and first studies showed a huge phenotypic variance in regard to growth at the onset of puberty [4, 23], milk yield [24] and temperament [25, 26]. By identifying eQTL in this population we assumed to gain insight into the genetic regulation of gene expression of different nutrition types that could provide new targets for further research of genes that have an effect on economically important traits such as milk production and meat quality traits in cattle.

The aims of this study were first to identify eQTL in the adrenal cortex in cattle and secondly, to evaluate whether associated genes and molecular pathways are known to affect economically important traits in cattle. The analysed data comprised genotypes of 37,204 SNPs and expression values for 10,986 transcripts. In order to account for systemic effects and to de-correlate the expression data a sire-dam model was used during pre-processing. To account for relatedness the complete pedigree data comprising the founder and F1 generation was used. Then we performed a genome-wide eQTL analysis using an adaptive linear rank test [27] on the residuals of the sire-dam model and the genotype data. Finally, we compared the results of this study with previous studies within the same population and performed Ingenuity® Pathway Analysis (IPA, QIAGEN Redwood City, www.qiagen.com/ingenuity) and literature reviews to identify the molecular pathways involved and evaluate potential targets for further research.

Results

Expression quantitative trait loci

In total, 505 SNPs were identified to be associated with the expression of 129 transcripts. A complete list of eQTL including the chromosomal position and annotation of the SNP and the gene encoding the transcript is provided in Additional file 1. In addition to Benjamini-Hochberg adjusted p-values of the adaptive linear rank test we also included the unadjusted p-values of adaptive linear rank test, Kruskal Wallis-test and ANOVA for completeness. The number of eQTL per chromosome is shown in Fig. 1. Identified associations and the corresponding genomic position for both SNP and the gene encoding the transcript are shown in Fig. 2. eQTL were identified on all chromosomes except on chromosomes X, 16, 24 and 28. The highest number of associations was observed on chromosome 23 comprising 55 SNPs associated with 11 transcripts. Apparently, in the case of most of these associations (482 of 505) both the SNP and the gene encoding the transcript were located on the same chromosome (cis associations). Only for 17 associations SNPs and genes were located on different chromosomes (trans associations). For the remaining six associations the position of the SNP was unknown.

Fig. 1

Distribution of eQTL across chromosomes. The number of cis (red) and trans eQTL (blue) per chromosome is shown as bars

Fig. 2

Position of eQTL in regard to the genomic location of the SNP and the gene encoding the transcript. The position of eQTL is plotted in regard to the genomic location of the SNP (x-axis) and the gene encoding the transcripts (y-axis). The size of the circles indicates the magnitude of the p-value with larger circles indicating smaller p-values. Red circles indicate cis and blue circles trans eQTL

Genomic distances between SNPs and genes for all 482 cis associations are shown in Fig. 3. Most of these SNPs (343 of 482) were located in close vicinity (less than 5 MB) to the gene encoding the transcript. On the other hand, the expression of LYPD6B (Chr. 2), GRIP1 (Chr. 5), SNCB (Chr. 7), COQ3 (Chr. 9) and ALDH5A1 (Chr.23), for example, was associated with 21, 15, 35, 28 and 20 SNPs, respectively. These SNPs were located up to 30 MB up- or downstream of the gene encoding the transcript. Considering that we investigated a F2 population that showed considerable linkage disequilibrium ([28] and Additional file 2) a clear discrimination between local and distant or cis and trans eQTL was hindered for these eQTL [29, 30]. Therefore, we assumed a cis or chromosome specific effect for all eQTL for which the SNP and the associated gene were located on the same chromosome. For the 17 eQTL for which the SNP and the gene encoding the transcript were located on different chromosomes a trans-regulatory effect was assumed.

Fig. 3

Distance between SNP and the gene encoding the affected transcript. The distance between the genomic location of the SNP and the start position of the gene encoding the transcript is shown for all 482 eQTL for which the SNP and the gene encoding the affected transcript were located on the same chromosome

Ingenuity® pathway analysis

In order to gain insight into molecular and cellular functions and regulatory processes, in which the genes encoding the 129 identified transcripts are involved in, Ingenuity® pathway analysis was performed. The analysis indicated that the genes are mainly involved in molecular and cellular function comprising free radical scavenging, cellular compromise, cell morphology, lipid metabolism and small molecule biochemistry. A complete list of the top five diseases and bio functions including the related genes is provided in Table 1. The results of the Ingenuity® Pathway upstream analyses are provided in Tables 2 and 3. Prolactin was the most important upstream regulator. The most important causal networks were related to ESRRG, GNRH1 and the two retinoic acid receptors RARG and RARB as master regulators (Table 3). Most of the affected diseases and bio functions as well as the most important causal networks and the most important upstream regulators are related to a relatively small number of genes that are involved in many of the identified pathways and causal networks, such as CYCS, LHCGR, TF, VIM, XDH, CYP27A1 and NQO1.

Table 1

Top five diseases and bio functions of the Ingenuity® pathway analyses

Category	p-value	Molecules
Molecular and Cellular Function
Free Radical Scavenging	2.69E-02-9.18E-02	TF, EPAS1, NQO1, CYCS, XDH
Cellular Compromise	2.69E-02-5.98E-02	VIM, NQO1, CYCS, XDH
Cell Morphology	2.69E-02-9.04E-02	VIM, MGAT5, LHCGR, NQO1, CYCS, GDPD2, XDH
Lipid Metabolism	2.69E-02-9.39E-02	VIM, EPAS1, NQO1, LHCGR, CYCS, XDH, CYP27A1, PHYH
Small Molecule Biochemistry	2.69E-02-9.39E-02	VIM, LYVE1, GSTO1, SLC2A9, EPAS1, NQO1, LHCGR, CYCS, XDH, CYP27A1, PHYH, EXTL2
Diseases and Disorders
Inflammatory Response	2.69E-02-9.39E-02	VIM, GSTO1, TF, MGAT5, EPAS1, NQO1, LHCGR, CYCS, XDH
Neurological Disease	2.69E-02-9.39E-02	VIM, GSTO1, TF, MGAT5, LHCGR, NQO1, XDH, CYP27A1, PHYH
Organismal Injury and Abnormalities	2.69E-02-9.48E-02	VIM, TF, PTP4A2, MGAT5, EPAS1, LHCGR, NQO1, XDH, PHYH
Psychological Disorders	2.69E-02-8.3E-02	VIM, GSTO1, TF, NQO1, CYP27A1, PHYH
Reproductive System Disease	2.69E-02-8.7E-02	TF, PTP4A2, LHCGR, PHYH
Physiological System Development and Function
Reproductive System Development and Function	2.69E-02-9.39E-02	TF, PTP4A2, EPAS1, LHCGR
Tissue Morphology	2.69E-02-9.39E-02	VIM, TF, EPAS1, LHCGR, NQO1, XDH, PHYH
Embryonic Development	2.69E-02-7.84E-02	VIM, LYVE1, TF, PTP4A2, EPAS1, LHCGR, NQO1, GDPD2
Cardiovascular System Development and Function	2.69E-02-8.7E-02	VIM, TF, MGAT5, EPAS1, CYCS, XDH
Hematological System Development and Function	2.69E-02-9.39E-02	TF, MGAT5, EPAS1, NQO1, LHCGR, CYCS, XDH

The top five molecular and cellular functions, diseases and disorders as well as physiological system development and functions are displayed, Benjamini-Hochberg adjusted p-values and the related genes

Table 2

Top five upstream regulators of the Ingenuity® upstream analyses

Upstream regulator	Molecule type	P-value of overlap	Target molecules in dataset
PRL	cytokine	6.07E-04	LHCGR, LYVE1, VIM, XDH
NRF1	transcription regulator	7.77E-04	CYCS, NQO1
MAPK7	kinase	1.13E-03	NQO1, VIM
NFE2L2	transcription regulator	2.57E-03	EPAS1, GSTO1, NQO1
CHD5	enzyme	2.80E-03	VIM

Table 3

Top four master regulators of the Ingenuity® upstream analyses

Master regulator	P-value of overlap	Target molecules in dataset
ESRRG	5.07E-05	CYCS, LHCGR, LYVE1, TF, VIM, XDH
GNRH1	1.44E-04	CYCS, CYP27A1, EPAS1, EXTL2, LHCGR, LYVE1, NQO1, P2RX5, PHYH, RPL13, SLC2A9, TF, VIM, XDH
RARG	1.73E-04	CYCS, CYP27A1, EPAS1, GSTO1, LHCGR, LYVE1, MGAT5, NQO1, TF, VIM, XDH
RARB	1.81E-04	CYCS, CYP27A1, EPAS1, GSTO1, LHCGR, LYVE1, MGAT5, NQO1, TF, VIM, XDH

Discussion

In order to identify molecular pathways that potentially influence the phenotypic variation of economically important traits in cattle we investigated the genetic background of adrenocortical gene expression in an experimental F2 population of a Charolais × German Holstein cross. We identified 505 SNPs associated with the expression of 129 genes. We also showed that some of these genes are involved in molecular pathways affecting economically important traits in cattle. In comparison with eQTL mapping studies in other species the proportion of cis eQTL (482) compared with trans eQTL (17) observed in our study is high. This is in line with the observation that the number of identified local or cis eQTL is related to sample size, extent of linkage disequilibrium, size of the window applied for cis eQTL and effect size of eQTL which is assumed to be stronger for cis eQTL in comparison with trans eQTL [29-31]. A reason for the high number of cis eQTL observed in our study might be that we assumed a cis or chromosome specific effect for all associations for which SNP and gene were located on the same chromosome. We decided to use this definition for cis associations for two reasons. Firstly, the extent of linkage disequilibrium observed in the studied population is vast ([28] and Additional file 2) and, secondly, the relative genomic position of SNP and associated gene does not necessarily provide information about the underlying regulatory variation [29]. However, even a window size of 1 MB which is often used to differentiate cis and trans effects would result in a higher number of cis eQTL (314) in comparison with trans eQTL (191). In addition, the extent of linkage disequilibrium and the high number of transcripts that are associated with more than one SNP in close vicinity to the gene encoding the transcript also indicated that the identified SNPs are not necessarily the causal variants but are in linkage disequilibrium to the causal polymorphisms and that the identified associations rather indicate genomic regions harbouring the causal variants affecting the transcript abundance.

Ingenuity® pathway analyses of genes encoding affected transcripts revealed that some of the genes are involved in “Diseases and Bio Functions” like lipid metabolism, inflammatory response and reproductive system development and function. This indicates that the identified eQTL potentially contribute to the phenotypic variation of economically important traits in cattle. For some of the genes there is evidence that they have an effect on reproduction, milk performance, growth or health traits, based on genetic and molecular studies in cattle and other farm animals. The following sections comprise details for selected genes. In addition, comparisons with studies within the SEGFAM population indicate that some of the genes affected by eQTL contribute to the phenotypic variation of behaviour, temperament and growth at the onset of puberty. Since the effects of adrenocortical gene expression on most of the economically important traits in cattle still have to be investigated, for other genes only the contribution to disease phenotypes in human or mice or the involvement in specific metabolic pathways like reproductive system development and function suggests that they might contribute to the phenotypic variation of economically important traits in cattle. In addition, the sheer existence of genomic regions affecting the adrenocortical transcript abundance of genes that are known to influence economically important traits in cattle is no evidence that the adrenocortical transcript abundance is causal, but similar effects might be present in other organs and tissues [29] providing insight into molecular mechanisms that contribute to the phenotypic variation in cattle. Subsequently we will discuss some of the genes that were highlighted by the functional analyses and that have been discussed to be involved in the phenotypic manifestation in cattle in other studies.

Genes involved in reproduction and other traits

Vimentin (VIM) is a type III intermediate filament that is expressed in cells of mesenchymal origin and has key functions in the formation of the cytoskeletal network, organelle positioning, cell migration and adhesions as well as in cell signaling [32]. VIM expression has been shown to be required for the development and establishment of the embryo in cattle. It was reported to be differentially expressed between nuclear transfer and by in vitro fertilization produced embryos [33, 34]. Nowadays in vitro production of embryos is common and approximately 15 % of all embryos are produced by in vitro technologies in cattle [35, 36].

The luteinizing hormone/chorigonadotropin receptor (LHCGR) is a G-protein-coupled receptor that is mainly known for binding luteinizing hormone and chorionic gonadotropin and mediating their cellular actions. In humans inactivating mutations in the LHCGR gene have been shown to cause infertility in both male and female [37]. In addition, LHCGR is discussed to be involved in the manifestation of the pregnancy induced Cushing’s Syndrome due to adrenocortical LHCGR overexpression that may lead to adrenocortical hyperfunction [38-40]. In cattle, polymorphisms within the LHCGR gene have been shown to be associated with calving interval, days to first service and milk yield [41] as well as with superovulation traits [42-44] that are important for generating a high number of eggs for the production of embryos [35, 36].

The cytochrome P450, family 27, subfamily A, polypeptide 1 (Cyp27A1) is a key-enzyme in Vitamin D3 metabolism [45] and is required for bile acid synthesis [46, 47]. In cattle, Vitamin D3 supplementation has been shown to affect dry-matter intake and beef tenderness [48]. In a previous study in the SEGFAM population, metabolites that belong to the sub class “Vitamine D3 and derivates” of Sterol Lipids were important for the classification of temperament types based on prefrontal cortex and serum metabolites [25]. In addition, likewise to LHCGR, Cyp27A1 has also been discussed to be important for reproduction [49] but has also been shown to affect milk yield and somatic cell score [50], an important production and an important health trait in dairy cattle, respectively.

Genes involved in milk production and other traits

The xanthine dehydrogenase (XDH) is a xanthine oxidoreductase that belongs to the group of molybdenum iron-sulfur flavin hydroxylases. XDH contributes to the detoxification of endogenous or xenobiotic compounds and is involved in the oxidative metabolism of purines [51, 52]. In cattle XDH has been shown to be upregulated during the lactation cycle due to its contribution to lipid droplet formation [53]. In addition, XDH protein abundance has been shown to increase in the first nine days of lactation in milk serum [54], and [55] suggested XDH to be a candidate gene for milk production and mastitis susceptibility based on studies in mouse and cattle.

Transferrin (TF) together with ferritin is important for acquisition, transport and storage of iron which is essential for many metabolic processes like the synthesis of hemoglobin for oxygen transport [56]. In cattle, TF has been shown to be highly polymorphic [57] and a study in Chinese native cattle indicated that TF polymorphisms are associated with protein yield, 305-day milk yield and mastitis susceptibility. In addition, they showed that TF mRNA expression was higher in mastitis affected in comparison with unaffected mammary tissue [58].

Genes potentially affecting behavior and temperament

The aldehyde dehydrogenase 5 family, member A1 (ALDH5A1) is a semialdehyde dehydrogenase that is important for the metabolism of γ-aminobutyric acid (GABA) the major inhibitory neurotransmitter in the brain [59]. A differential expression of ALDH5A1 has been reported in Angus cattle selected for high or low residual feed intake [60]. Increased levels of 2,4-dihidroxy-butanoic acid have been observed in the urine of patients suffering from semialdehyde dehydrogenase deficiency [61] and an accumulation of gamma-hydroxybutyric acid (GHB) and GABA in the central nervous system was reported in ALDH5A1−/− mice [62]. 2,4-dihidroxy-butanoic was priviously identified to be a potential biomarker in Alzheimer’s disease [63] and a temperament type related abundance of 2,4-dihidroxy-butanoic acid as well as of γ-aminobutyric acid acid has also been observed in temperament type specific metabolite profiles of the prefrontal cortex of animals with extrem temperament types deriving from the SEGFAM population [25]. Interestingly ALDH5A1 is located on chromosome 23 in the genomic region were bovine MHC genes are located. In this genomic region we identified a high number of SNPs that were not only associated with the expression of ALDH5A1 but also with the expression of some MHC genes like the MHC class I heavy chain BOLA and the major histocompatibility complex BOLA-DQA1. This genomic region has been reported to be associated with the antibody mediated immune response in Canadian Holstein cows [64].

Another study analyzed the impact of genetic variants on behaviour characteristics assessed in an open field and novel object test within the SEGFAM population [28]. We identified three SNPs (ARS-BFGL-NGS-2942/rs110027993, ARS-BFGL-NGS-98658/rs109674592, ARS-BFGL-NGS-27299/rs109313646) that affect the expression of NME6 and RAB32 in this study. The same SNPs were previously identified to be associated with the inactivity in an open field test and the explorative behaviour in a novel object test, respectively [28].

The NME/NM23 nucleoside diphosphate kinase 6 (NME6) is a nucleoside diphosphate kinase. Little is known about its function besides that it is assumed to play a role in cell growth and cell cycle progression [65, 66]. In addition, a recent study in pig indicated that NME6 is located in a QTL region associated with back fat [67].

RAB32, member RAS oncogene family (RAB32) belongs to the RAB family of small GTP-binding proteins. It is involved in the biogenesis of lysosome-related organelles like melanosomes [68]. Similar to NME6 we could not find any additional information other than a suggested role in lipid metabolism [69].

Ingenuity® upstream analysis

The upstream analysis revealed PRL, ESRRG and GNRH1 as upstream or master regulators, providing insight into regulatory mechanisms that are potentially affected by eQTL.

Prolactin (PRL) is an anterior pituitary protein hormone that is not only important for reproduction but is also involved in the control of behaviours [70]. In cattle especially the effects of PRL in regard to lactation and milk production have been investigated. In the last years several studies in cattle investigated the inhibition of PRL-release on milk yield and showed that the inhibition was associated with a decrease in milk yield additionally indicating an important role of PRL in the maintenance of lactation [71-73]. The SEGFAM population showed remarkable differences in their milk production [24].

Gonadotropin-releasing hormone 1 (GNRH1) encodes a protein that is secreted and cleaved to form the luteinizing hormone-releasing hormone and the prolactin release-inhibiting factor that regulate LH and FSH gene expression [74]. LH and FSH both have been early recognized to be involved in the onset of puberty in cattle and today Gonadotropin-releasing hormone is used to synchronize cows for timed artificial insemination especially in beef cattle [75, 76]. In addition, GnRH signaling was previously identified as a key modifier of differential growth at the onset of puberty in the SEGFAM population [4].

The Estrogen-Related Receptor Gamma (ESRRG) that belongs to the estrogen receptor-related receptor family has recently been shown to be important in the fibrate-mediated regulation of lipid metabolism genes in a human ApoA-I transgenic mice model [77]. Interestingly ESRRG was also identified to be associated with the onset of puberty in cattle [78].

Conclusions

In the present study we investigated the genetic background of adrenocortical gene expression and identified eQTL affecting the expression levels of genes that have been reported to influence economically important traits in cattle. In addition, some of the genes and molecular pathways affected by eQTL were previously shown to contribute to the phenotypic variation of behaviour, temperament and growth at the onset of puberty in the SEGFAM population. Therefore, eQTL analysis appears to be a useful approach to gain insight into the molecular and genetic background of complex traits in cattle. Additional studies are needed to identify the causal polymorphisms and underlying molecular mechanisms (cis or trans) and the effects of genetic variation on gene expression in other tissues.

Methods

Animals and tissue samples

This study comprised 145 F2 cows from the experimental resource cross population (SEGFAM) that was initiated at the Leibniz Institute for Farm Animal Biology (FBN) in the P0 generation by crossing five purebred Charolais sires to purebred German Holstein cows [22]. The population was established by multiple ovulation and embryo transfer, which resulted in five half-sib and full-sib families from which 49, 35, 52, 8 and 1 cows were included in this study, respectively. The cows from the fourth and fifth family were maternally related to animals in the first and third family, respectively. All animals were housed in a loose housing barn under identical environmental and feeding conditions at the Leibniz Institute for Farm Animal Biology (FBN) in Dummerstorf, Germany [24]. The animals were slaughtered at day 30 after parturition in their second lactation following a standardized protocol. Besides other tissues mammary gland tissue and the adrenal gland were immediately taken after slaughter for DNA and RNA extraction, respectively. The adrenal gland was further dissected to separate the adrenal cortex from the medulla. All tissue samples were cut in small pieces, snap-frozen in liquid nitrogen and stored at −80 °C or in liquid nitrogen.

DNA extraction and SNP genotyping

SNP genotyping was accomplished using Illumina® BovineSNP50 Beadchip v1.0 and v2.0. Genotype calling was performed with the Genotyping Module of the GenomeStudio V2011.1 Software (Illumina®). For this purpose genomic DNA was extracted from mammary gland tissue with the QIAmp DNA Mini Kit (QIAGEN, Hilden, Germany) and prepared for hybridization on the Illumina® BovineSNP50 Beadchip following manufacturer’s instructions. Genotype calling and quality control was performed according to [79] for v1.0 and v2.0 data separately. After quality control and filtering the v1.0 and v2.0 datasets were merged and all samples and SNPs with more than 10 % missing genotypes and a minor allele frequency of less than 5 % were excluded from further analyses. To identify inconsistencies between recorded genotypes and pedigree information, the software PedCheck [80] and a larger dataset including genotypes of P0 sires and all (male and female) F1 and F2 animals was used [81]. The final dataset comprised 37,204 SNPs of 145 animals. Chromosomal positions of SNPs are taken from [82] and are based on genome assembly UMD3.1.

RNA extraction and microarray hybridization

Total RNA was isolated from the adrenal cortex of 145 animals using TRI Reagent (Sigma, Taufkirchen, Germany). After DNase I treatment the RNA was further purified with the RNeasy Kit (QIAGEN, Hilden, Germany) following the manufacturers recommendations. RNA was quantified using the NanoDrop ND-1000 spectrophotometer (Peqlab, Erlangen, Germany). Integrity of RNA was checked by running 1 μg of RNA on a 1 % agarose gel. Absence of DNA contamination was verified by PCR amplifying fragments of GAPDH with RNA as a template. For hybridization, 500 ng of total RNA were amplified using Ambion WT Expression Kit (Affymetrix). Subsequently, 5.5 μg of the resulting cDNA was fragmented and labeled using the Affymetrix Terminal Labeling Kit. The fragmented cDNA was hybridized to the microarray using the Affymetrix Hyb-WashStain Kit and Affymetrix standard protocols. Fluidic station protocol was FS450_0001.

For expression profiling the custom GeneChip Bovine Gene v1 Array was used. Affymetrix designed the GeneChip Bovine Gene v1 Array based on Ensemble and RefSeq predictions for the Genome Bos Taurus Build 4.0. The design was targeted to develop a whole Genome Expression Array with approximately 25 probes per transcript distributed over the whole length of each transcript. In total 194,712 probe sets targeting at approximately 24,000 transcripts are implemented on the array. In addition standard Affymetrix controls for hybridization and labeling efficiency as well as for non-specific binding are included on the array. Gene-level analysis was performed using Affymetrix® Expression Console™ Software. Microarray raw data were preprocessed using the RMA algorithm [83]. Quality control was performed in accordance to [84]. In addition, the detection above background algorithm was used to identify expressed genes. Probe sets that were present in 75 % of all samples and transcripts with at least 50 % of all probe sets present were included in the analyses. The final dataset comprised expression values for 10,986 transcripts of 145 animals. For annotation the NetAffx [85] release 34 chromosomal positions of the transcripts were used that are based on genome assembly UMD3.1. In addition the Ingenuity® gene annotation is provided including information about transcript cluster ids that target the same gene. The microarray data are deposited at Gene Expression Omnibus database [86] (GEO: GSE75371).

Data pre-processing

After quality control and filtering the expression data was further pre-processed to account for systemic effects. Furthermore, the statistical test used for eQTL analysis required the data to be decorrelated. Since we expected gametic phase disequilibrium for some of the loci, we refrained from estimating kinship from SNP data but preferred the use of kinship estimates from pedigree information. Therefore the following sire-dam model was applied in ASReml [87] to account for fixed effects of season and year and age at slaughter as well as for relatedness:

where y is the expression value for transcript y of animal i (i = 1,…, 145), s is the fixed effect of year and season (j = 1,…, 25) at slaughter, x is the age in days of animal i at slaughter, a is the additive-genetic effect of the sire of animal i, a is the additive-genetic effect of the dam of animal i and ε is the random residual effect. To estimate the additive-genetic effects the complete pedigree data comprising 38 P0 and 52 F1 animals was used. The residuals comprise the Mendelian sampling effect as an uncorrelated genetic component as well as a random deviation. For further analysis the residuals were considered as de-correlated expression data.

eQTL analysis

For eQTL analysis an adaptive linear rank test [27] was applied. The test was optimized to meet the challenge of non-normally distributed expression data in eQTL studies. This two stage procedure first calculates selector statistics for skewness and tail length of the de-correlated expression data. Then an appropriate linear rank test depending on the selector statistics is applied to identify differences in the expression levels per genotype group. Thus, for each pair of SNP and transcript a linear rank test is selected that is supposed to have large statistical power for the respective actual distribution. The available tests are Kruskal-Wallis test [88], median test [89], long tails test [90], short tails test [91], right skewness and left skewness test [92]. To account for multiple testing we applied a global significance level of 0.1 for Benjamini-Hochberg adjusted p-values [93] (considering all 374 million performed tests), which is approximately equivalent to a significance level of 1.3x10−7 for unadjusted p-values. This global false discovery rate (FDR) threshold of 0.1 is less conservative for trans associations than usual, but it is very conservative for cis associations, since it corresponds to a FDR threshold of approximately 0.0035 for an eQTL analysis in which only cis associations (approximately 13 million) are considered for the calculation of the Benjamini-Hochberg adjusted p-values. This was done because the majority of associations was expected between SNPs and transcripts for which the encoding gene of the transcript was located on the same chromosome in close vicinity to the SNP.

The functional and upstream analyses were performed through the use of QIAGEN’s Ingenuity® Pathway Analysis (IPA, QIAGEN Redwood City, www.qiagen.com/ingenuity). The purpose of the functional analysis was to determine if the 129 identified genes are statistically significant over-represented (more genes than expected by chance) for a specific molecular or cellular function in comparison to a reference-set of genes. Similar to the functional analysis, the Ingenuity® upstream analysis was performed to measure whether there is a statistically significant overlap between the 129 identified genes and the genes that are regulated by a specific transcription regulator. For both tests IPA uses a Fisher’s Exact Test to determine the significance of the over-representation or overlap in comparison to a reference set. As reference set for the functional analysis the Ingenuity® Knowledge Base was chosen, because the reference network that is used to identify the upstream and master (upstream and downstream) regulators in the upstream analysis is based on the Ingenuity® Knowledge Base. Reported p-values are Benjamini-Hochberg adjusted.