Abstract
Cortisol is a steroid hormone with important roles in regulating immune and metabolic functions and organismal responses to external stimuli are mediated by the glucocorticoid system. Dysregulation of the afferent and efferent axis of glucocorticoid signaling have adverse effects on growth, health status, and well-being. Glucocorticoid secretion and signaling show large interindividual variation that has a considerable genetic component; however, little is known about the underlying genetic variants. Here, we used trait-correlated expression analysis, screening for expression quantitative trait loci (eQTL), genome-wide association (GWA) studies, and causality modeling to identify candidate genes in porcine liver and muscle that affect or respond to plasma cortisol levels. Through trait-correlated expression, we characterized transcript activities in many biological functions in liver and muscle. Candidates from the list of trait-correlated expressed genes were narrowed using only those genes with an eQTL, and these were further prioritized by determining whether their expression was predicted to be related to variation in plasma cortisol levels. Using network edge orienting (NEO), a causality modeling algorithm, 26 of 990 candidates in liver were predicted to affect and 70 to respond to plasma cortisol levels. Of 593 candidates in muscle that were correlated with cortisol levels and were regulated by eQTL, 2 and 25 were predicted as effective and responsive, respectively, to plasma cortisol levels. Comprehensive data integration has helped to elucidate the complex molecular networks contributing to cortisol levels and thus its subsequent metabolic effects. The discrimination of up- and downstream effects of transcripts affecting or responding to plasma cortisol concentrations improves the understanding of the biology of complex traits related to growth, health, and well-being.
- expression quantitative traits loci (eQTL)
- genome-wide association (GWA)
- network edge orienting (NEO)
- quantitative trait networks
- cortisol concentration
CORTISOL, a steroid hormone released in response to stress and a low level of blood glucocorticoids, acts to increase blood sugar through gluconeogenesis; suppress the immune system; and aid in fat, protein, and carbohydrate metabolism. Produced in the adrenal cortex under the control of the hypothalamic–pituitary–adrenal axis and the sympatho-adrenomedullar (SAM) system, cortisol is most important during times of stress, when it modulates many of the acute responses to illness and “resets” metabolism in favor of providing substrates for oxidative metabolism (Andrews and Walker 1999). Stress factors can unbalance metabolic homeostasis, directing the demand for nutrients and energy for essential physiological processes away from growth (Heetkamp et al. 1995). Individual cortisol concentration is variable and depends on sensitivity to stress, which is highly genetically determined (e.g., heritability estimates of 0.40–0.70 in pigs) (Geverink et al. 2002; Kadarmideen and Janss 2009). According to the function of cortisol in the maintenance of basal and stress-related homeostasis the adrenocortical steroidogenesis is sensitive to exogenous biotic and abiotic stressors but also to endogenous signals released from peripheral tissues and organs different from the hypothalamic–pituitary–adrenal (HPA) axis or SAM system. Such effective signals are associated with the plasma cortisol concentration that depends on the rate of adrenal steroidogenesis regulated either directly at the level of genes encoding the corresponding enzymes or indirectly by triggering the HPA or SAM systems. Glucocorticoids mainly exert their effects via members of the nuclear receptor superfamily of ligand-dependent transcription factors, primarily the glucocorticoid receptor (GR) and to a lesser extent through the mineralocorticoid receptor (MR). The steroid-receptor complex binds to hormone-responsive elements on the chromatin and regulates gene transcription; i.e., genes addressed by this effective mechanism are responsive to cortisol plasma concentration.
Cortisol concentrations in plasma have been shown to be both directly and indirectly related to behavioral, physiological, and metabolic disease in mammals (de Groot et al. 2001; Murani et al. 2010, 2011). However, the association patterns of gene expression with cortisol concentrations remain little understood.
Pigs, valued as agricultural commodities, are typically raised in intensive production units that likely introduce various stresses to the animals (Black et al. 2001). Stress can negatively affect growth, meat quality, and the immune system. Importantly, pigs share both physiological and genomic similarities with humans and therefore provide a more tractable model in which to study genetic determination of behavioral, physiological, and metabolic traits (Bode et al. 2010). As stress can also negatively affect human health, the identification of genes related to cortisol concentrations in pigs may enhance the discovery of genes associated with metabolic, immunological, and behavioral traits in humans.
Processes occurring in both the liver and muscle can be influenced by cortisol concentrations. Indeed, the liver is of particular interest given its vital roles in maintaining homeostasis and health as well as regulating nutrient utilization. On the other hand, muscle is a major consumer of energy and is responsive to cortisol in terms of its metabolic and motor activities. Thus, we used expression analyses in liver and muscle to identify genes related to plasma cortisol concentrations to better understand the biological processes and underlying genes and functional networks influenced by this hormone in target tissues.
Functional genomics provide insight into essential cellular functions whose regulation determines a significant proportion of the phenotypic variance. Transcript profiling is a proximal endophenotype affected by genetic variation (Jansen and Nap 2001) and has been used to dissect the genetic architecture of expression regulatory variation in a number of systems (Morley et al. 2004; Cheung et al. 2005; DeCook et al. 2006; Schadt et al. 2008; Swanson-Wagner et al. 2009; Ponsuksili et al. 2010, 2011) and has aided in the identification of quantitative trait loci (QTL) (Bystrykh et al. 2005; Hubner et al. 2005). Pairwise associations between gene expression traits and other organismal traits can be used to define an undirected traits network. Recently, causality modeling algorithms have been developed by using genetic markers as causal anchors for orienting the edges of a trait network. These algorithms proved effective in establishing causality by integration of gene expression, QTL, and association (Farber et al. 2009a,b; Plaisier et al. 2009). Using this approach, we determined the biological function of genes in liver and muscle correlated with plasma cortisol concentrations. We then identified candidate genes in these target tissues that affect or respond to plasma cortisol concentrations by integrating trait-correlated expression analysis, detection of expression quantitative trait loci (eQTL), genome-wide association (GWA) study, and causality modeling. Finally, causal links between SNP, gene expression, and cortisol concentration were generated.
Materials and Methods
Animals, tissue collection, and phenotyping
Animal care and tissue collection procedures followed the guidelines of the German Law of Animal Protection and the experimental protocol was approved by the Animal Care Committee of the Farm Animal Biology (FBN). This study was based on trait measurements, genotyping records, expression profiles, and genome-wide association analyses done with performance-tested pigs from commercial herds of the crossbreed Pietrain × [German Large White × German Landrace (female and male castrates (53%)]. Total cortisol secretion was measured in the morning by collection of 50 ml of trunk blood from each pig during exsanguinations at an age of 180 ± 2 days. Blood was collected in a plastic tube containing 1 ml of 0.5 M EDTA, stored on ice until plasma preparation, and then stored long-term at −80°. Plasma cortisol concentrations (total) were determined using commercially available enzyme-linked immunosorbent assay (DRG, Marburg, Germany) in duplicate according to the manufacturer’s protocol. The intra- and interassay coefficients of variation were ≤7.0% and 9.8%, respectively.
Whole-genome scan and quality control
Illumina bead array technology was used to carry out all genotyping reactions in accordance with the manufacturer’s protocol for the SNP Infinium HD assay (http://www.illumina.com). Genotyping was performed using the PorcineSNP60 BeadChip (Illumina, San Diego). In brief, 200 ng of DNA was used for genome-wide amplification and subsequent fragmentation. DNA was hybridized to the 62,163 locus-specific 50mers covalently linked to the beads distributed on the surface of the microarray. Single-base extension of oligos on the BeadChip, using captured DNA as a template, was performed, incorporating detectable labels on the BeadChip. Signals of each wavelength were determined using an Illumina iScan that converted the images to intensity data. Intensity data for each SNP were normalized and assigned a cluster position and a resulting genotype with the GenomeStudio software (Illumina); a quality score was generated for each genotype. An overlapping set of 150 (liver) and 207 (muscle) samples was genotyped for 62,163 SNPs, and quality of the data was evaluated. Samples with call rates <95% were removed. Markers were excluded if they had low minor-allele frequency (MAF) <5%. Deviation from Hardy–Weinberg equilibrium was not considered because a three-way crossbreed pig population was used; deviation from Hardy–Weinberg equilibrium can be expected due to discordant allele frequencies in the parental breeds. The average call rate for all samples was 99.8% ± 0.2%.
Whole-genome expression profiling
Gene expression profiling of the liver and M. longissimus dorsi samples of pigs was conducted with the same 150 and 207 animals. In brief, total RNA was isolated using TRI Reagent (Sigma, Taufkirchen, Germany) and used for target preparation for microarray hybridization. According to Affymetrix protocols, 500 ng of total RNA was reverse transcribed into cDNA, transcribed into cRNA, and labeled using the Affymetrix One cycle synthesis and labeling kit (Affymetrix, High Wycombe, UK) to prepare antisense biotinylated RNA targets. Liver expression patterns were produced using 150 GeneChip Porcine Genome Arrays (Affymetrix). Quality of hybridization was assessed in all samples, following the manufacturer’s recommendations. Data were analyzed with the Affymetrix GCOS 1.1.1 software, using global scaling to a target signal of 500. Data were processed with MAS5.0 to generate cell intensity files (present or absent). Quantitative expression levels of the present transcripts were estimated using the probe logarithmic intensity error (PLIER) algorithm for normalization that was implemented in the Expression Console software (Affymetrix). Normalizations were performed separately per tissue. The microarray data related to all samples were deposited in the Gene Expression Omnibus (GEO) public repository [GEO accession nos. GSE25445 (liver) and GSE32112 (muscle)]. The overlapping set of successfully genotyped samples with available expression data amounted to 150 (liver) and 207 (muscle) individuals.
Correlation between traits and expression levels
Phenotypes and expression levels were adjusted for systematic effects by analysis of variance performed with the procedure “Mixed” of the SAS software package (SAS version 9.1; SAS Institute, Cary, NC) before analyzing their correlation. “Sex” was used as a fixed effect, “sire” and “slaughter day” as random effects, and “carcass weight” as a covariate. In muscle an additional fixed effect of RYR genotype was used. Subsequently, Pearson correlation coefficients were calculated between the residuals of log2-transformed expression intensities and cortisol concentrations in plasma. Genes that showed correlation at P ≤ 0.05 were analyzed further.
Pathway analysis and network generation
Based on BLAST comparison of the Affymetrix porcine target sequences with the porcine genome sequence (Ensembl_Sscrofa_10), 20,689 of the 24,123 probe sets on the Affymetrix Porcine GeneChip were localized and annotated. The list of significant trait-correlated transcripts was analyzed, referring to predefined pathways and functional categories of the Ingenuity Pathways Knowledge Base, using Ingenuity Pathways Analysis. Also, canonical pathways were identified from the Ingenuity Pathways Analysis library, which were most significant to the input data. The significance of the association between the data set and the predefined pathways and functional categories was measured by Fisher’s exact test and adjusted using the Benjamini–Hochberg correction (Benjamini and Hochberg 1995).
Population structure analysis
We first computed the genetic similarity matrix between the individuals as identity by descent (IBD) of each pair and visualized the result as a heatmap in Figure 2. A complex multileveled population structure and genetic relatedness are observed in the genetic similarity matrix. The QK mixed model was used to test for association between traits and SNP genotypes from a single SNP at a time while adjusting simultaneously for population structure and family relatedness (Yu et al. 2005). A regression testing for linear trend of SNP alleles was performed.
Genome-wide association
For the association analysis, cortisol concentration and log2-transformed expression levels of individuals were used as traits for whole-genome association analysis with SNPs genotyped by a mixed-model analysis of variance, using JMP Genomics (SAS Institute). “Genotype” and “gender” were used as fixed effects, “relationship matrix” and slaughter day as random effects, and carcass weight as a covariate. In muscle an additional fixed effect of RYR genotype was used. The sequences flanking the SNPs represented on the Illumina SNP chip, which were significantly associated with the level of transcription of any probe set of the Affymetrix expression microarrays, were assigned to the porcine genome sequences (Ensembl_Sscrofa_10). Annotation and localization of SNP sites and probe sets allowed discrimination of cis- and trans-regulation. We defined an eQTL as cis if an associated SNP was located within an area <10 Mb from the probe set/gene. All other eQTL were considered as trans. Results are reported at thresholds of P < 10−5, corresponding to a genome-wide false discovery rate (FDR) of 11% in both tissues. For the FDR calculation we used the q-value package in R (Storey and Tibshirani 2003). In total 80,725,037 P-values in liver and 45,625,920 P-values from muscle were used to estimate FDR. A significance threshold was used, representing a trade-off between avoidance of false-positive associations while taking into account the likely higher dependence among the expression traits (correlated expression) and among the SNP genotypes (linkage disequilibrium).
Causality modeling using network edge orienting
To infer causal flow between genetic markers, gene expression traits, and cortisol concentration, we utilized the network edge orienting (NEO) package (Aten et al. 2008). We used the residuals of gene expression and cortisol concentrations after correction for systematic effects as previously described. Because we were investigating a large amount of expression traits, to avoid the high complexity of the unknown quantitative traits network, we looked for prediction of pairwise causal relation only between each expression trait and cortisol concentration. On the other hand, concerning the genetic marker selection, we adopted a hybrid method combining manual selection and automatic selection embedded in NEO. For each pair of expression trait and cortisol concentration, a pool of SNPs is constructed using the SNPs found to be significantly associated with the expression level (eQTL) and cortisol concentration. Markers were first divided into two larger groups for each trait with respect to marker assignment consistency to ensure clear separation of two nonoverlapping groups of markers, and the best sets of markers were then picked up using multivariate regression model variable selection with respect to the Akaike Information Criterion (AIC). Before fitting the local structural equation models (SEM) model, the local SEM-based edge-orienting next-best (LEO.NB) orientation score was first evaluated for both directions and proceeded only with the edge direction having a LEO.NB orthogonal causal anchor (LEO.NB.OCA) score >0.3 and a significant P-value <0.05 for the Wald test of path coefficient. The recommended threshold of 0.3 for LEO.NB.OCA is substantially more conservative than the 5% false-positive rate. A local SEM model was therefore fitted for the selected direction of the edge, using the best genetic markers selected, and model probability was then evaluated. The causal model was dropped if a model probability <0.05 was detected. Further evaluation of causality could proceed by inspecting the model-fitting index root mean-square error of approximation (RMSEA) to see whether a smaller value was in agreement. As recommended by the package author, a robustness analysis was also applied, for maximum anchoring SNP number ranging from 1 to 5.
Results
We first determined the distribution of plasma cortisol in our population (N = 475) of performance-tested pigs from commercial herds of the crossbreed Pietrain × (German Large White × German Landrace); mean (± SD) cortisol was 93.9 (±34.6) ng/ml (Figure 1). To characterize the architecture of hepatic gene expression influenced by cortisol, a subset of 150 performance-tested pigs was analyzed by Affymetrix Porcine Genome Arrays containing 24,123 probe sets, of which 20,689 probe sets assigned to known genes were applied for transcription profiles. Similarly, gene expression was quantified in M. longissimus dorsi, another peripheral target for cortisol, in 207 performance-tested pigs from the same herd (132 of these overlapped with the set used for hepatic expression profiles). MAS5 analysis revealed consistent “present calls” for 10,945 and 11,191 probe sets in liver and muscle, respectively; these probe sets represent genes being “active”. Probe sets with present calls were further analyzed using the PLIER algorithm for normalization that was implemented in the Expression Console software (Affymetrix) and subsequently logarithmized PLIER values were subjected to association analysis.
The distribution of plasma cortisol concentrations (nanograms per milliliter) of all 475 individuals with mean ± SD of 93.9 ± 34.6 ng/ml.
The study population (N = 475) was also genotyped using the PorcineSNP60K BeadChip (Illumina). SNPs with genotyping call rates >95% and minor-allele frequencies >5% were selected for GWA analysis. The population structures in performance-tested pigs from commercial herds were estimated. One approach to using genetic markers to elucidate the relationships between individuals is genome-wide relatedness as calculated from IBD (Oliehoek et al. 2006). High IBD was found for paternal half-sibs (Figure 2) and was considered accordingly as a random factor before further analysis.
IBD relationship matrix heat map between individual pigs (n = 475), reflecting the population structure among performance-tested pigs from commercial herds. The dendrograms on the bottom and on the side of the heat map are based on the hierarchical clustering of identities by descent (IBDs) assessed from genome-wide genotypes. The genome-wide relatedness was calculated from IBD.
Identification of gene expression traits correlated with cortisol concentrations
In liver, 10,945 probe sets represented active genes. Analysis of correlation between residuals derived from the mixed-model analysis of cortisol concentrations and gene expression levels revealed 1531 probe sets having significant correlation of their expression with cortisol concentration. The correlation coefficients of expression and cortisol concentrations—adjusted for systematic effects and significant at P < 0.05—ranged between |0.16–0.42|. Cortisol concentrations were positively correlated with 624 probe sets and negatively correlated with 907 probe sets (Supporting Information, Table S1). Similarly, in muscle, of 11,191 probe sets representing active genes, 960 had expression levels correlated with plasma cortisol concentration (Table S2). Of these, 664 probe sets were positively correlated and 296 were negatively correlated, with r ranging between |0.31–0.14|. Cortisol concentrations were most strongly correlated with ARL4A in liver (r = 0.42, P = 1.12 × 10−7, FDR = 0.0009) and NFKBIA in muscle (r = 0.31, P = 6.99 × 10−6, FDR = 0.05). Further, 135 probe sets showed correlation with cortisol concentrations and were expressed in both liver and muscle. Indeed, 78 of 135 transcripts had the same directions of correlation between cortisol concentration and expression in these tissues (Table S3).
Functional annotation of liver transcripts correlated to plasma cortisol concentrations
We tested the significant positively and negatively correlated functional candidate genes for enrichment in functional annotation groups, as defined in the Ingenuity Pathways Knowledge Base. In liver, genes represented by probe sets with expression levels that correlated positively with cortisol concentrations were often involved in metabolism of lipids and carbohydrates as well as various macromolecules; genes with negative correlation were generally involved in dynamic cellular processes (Table S4). To further refine the functional annotation of gene sets with cortisol concentrations correlated to expression in liver, their assignment to canonical pathways was explored using the Ingenuity Pathway Knowledge Base. The positively and negatively correlated expression levels of functional candidate genes were uploaded as gene data sets for analysis. Significantly enriched canonical pathways for positively correlated transcripts included cytokine signaling [acute phase response signaling, interleukin (IL)-10, and IL-6], glucocorticoid receptor signaling, amino acid metabolism (glutamate, β-alanine, tryptophan, arginine and proline, and methionine), carbohydrate metabolism (propionate and butanoate), and fatty acid metabolism. Most of the negatively correlated transcripts were categorized in cellular immune response (antigen presentation pathway, interferon signaling, cytotoxic T-lymphocyte–mediated apoptosis of target cells, and IL-4 signaling), and phospholipase C signaling pathways.
Functional annotation of muscle transcripts correlated to plasma cortisol concentrations
Transcripts in muscle positively correlated with cortisol concentrations were categorized as genes involved in protein synthesis and functions of dynamic cellular processes. Transcripts negatively correlated with cortisol concentrations were less often assigned to biological functions, with just 30 genes involved in cell cycle and 44 in cell death (Table S5). Three canonical pathways (NRF2-mediated oxidative stress response, HMGB1 signaling, and IL-17A signaling in fibroblasts) were significantly enriched for genes positively correlated with plasma cortisol. Most of the negatively correlated transcripts belonged to PDGF signaling or an antiproliferative role of TOB in T-cell signaling.
Detection of cis and trans-eQTL in liver
To characterize the genetic architecture of hepatic gene expression relevant to cortisol concentrations, only correlated expression levels were used for further analysis. eQTL represent variations perturbing the expression of a gene (Rockman and Kruglyak 2006). eQTL of genes for which expression levels are correlated to traits are strong positional and functional candidates. To identify positional and functional candidates regulating cortisol concentrations, we performed an eQTL analysis, using microarray expression profiles of liver and muscle tissue from our population of pigs. Two types of eQTL were distinguished: those that mapped near (<10 Mb from) the gene encoding the transcript (cis) and those that mapped elsewhere in the genome (trans) (Table 1). eQTL were identified for the 1531 probe sets with trait-correlated expression in liver, with 52,727 SNP markers using the mixed model with systematic adjusting (see Materials and Methods) (Figure 3). This analysis identified 6608 eQTL P < 10−5 (11% FDR): 1019 in cis and 4873 in trans. The mean negative log values of the P-values of cis and trans-eQTL were 8.18 and 6.34, respectively. The 4873 trans-eQTL and 1019 cis-eQTL corresponded to 927 and 116 genes, respectively (Table 1). Thus, ∼64% of genes with expression correlated with plasma cortisol were associated with eQTL; 12% of these eQTL act in cis compared to 94% regulated in trans. The 6608 SNPs explain 12–64% of the variation in transcript levels, with, on average, 23% for cis-eQTL and 17% for trans-eQTL.
Manhattan plots displaying the significance of associations of SNP markers and probe sets with trait-associated expression in liver (top) and muscle (bottom).
Detection of cis and trans-eQTL in muscle
In muscle, eQTL analysis was performed for 960 probe sets representing genes whose expression correlated with plasma cortisol concentration. We used the same mixed model for analysis of liver data but with the genotype at the SNP c.1843C > T of the porcine ryanodine receptor gene (RYR1) as an additional fixed effect because of its known impact on muscle metabolism (Figure 3). A total of 4115 eQTL were identified at P < 10−5 (11% FDR), with a mean negative log of the P-value of 7.44. We detected 2001 cis and 1663 trans-eQTL with mean negative log of the P-values of 8.62 and 6.02, respectively. These eQTL corresponded to 593 genes: 139 in cis and 538 in trans. Similar to liver, in muscle 62% of transcripts correlated with cortisol had associated cis- or trans-eQTL. However, in muscle a higher proportion of cis-eQTL, 23%, was found; 91% of the eQTL had trans effects. The eQTL associated with the 4115 SNPs explained 9–39% of the variation in transcript levels, with, on average, 15% for cis-eQTL and 11% for trans-eQTL.
Predicting causal links between genes with eQTL and cortisol concentrations
eQTL of genes showing expression levels correlated with plasma cortisol concentration potentially point either to key drivers of changes in cortisol concentrations or to transcripts responding to cortisol concentrations. However, eQTL and correlation alone are not enough to indicate the direction of the relationship and to support causality. To predict causal relationships between gene expression traits and cortisol concentrations, we used the NEO R software package (Aten et al. 2008), which takes SNP genotypes as input into structural equation models that compute local-structure edge-orienting scores (LEO.NB). Here, we used SNPs associated with the expression levels of probe sets (eQTL markers) and SNPs associated with plasma cortisol concentrations (GWA markers) from the same population to predict causal interactions between the two traits, i.e., gene expression and cortisol concentrations.
In total, 6608 markers from eQTL in liver and 4115 markers from eQTL in muscle, at LOD scores >5, as well as 11 significant GWA markers (Table 2), at P < 10−4, were used to predict the direction of the relationship of 990 and 593 genes with cortisol concentration-correlated expression that showed eQTL in liver or muscle, respectively. NEO was used to study causal relationships between expression levels of genes and cortisol concentrations, i.e., either gene expression → cortisol (gene expression effects change in cortisol concentration) or cortisol → gene expression (cortisol concentration changes gene expression), to identify the up- and downstream molecular pathways of cortisol. In liver, 26 of 990 genes (Table S6) were predicted as causal (SNP → expression→ cortisol ← SNP). Twenty genes predicted as causal together with their function are presented in Table 3. The most significant prediction of a causal link was for CD99 (Ssc 3205.1.S1.at), with a LEO.NB.OCA score of 1.07. The three genes Pla2g10, LSR, and muCL had cis-eQTL (Table 3). Another 70 genes (Table S7) were predicted to be regulated by cortisol (SNP → cortisol → expression ← SNP). Among these genes being responsive to plasma cortisol concentrations, those annotated were involved in lipid metabolism (SGK1, ABCA3, LIPG, ACSL1, APOC3, PCTP, TLR4, NFKBIA, ITGB3, RDH11, and INSIG1), transcription regulation (CCRN4L, TAF1A, CEBPB, EPAS1, TEF, AKIRIN2, CEBPD, YY1, SREBF2, and NFIB), and cell death (ASAP3, TP53INP1, SLC30A10, GAS1, DBNL, GNB1, MCTS1, CFLAR, POLB, CYR61, and p53RFP) (Table 4). The most significant responsive gene was CCRN4L (Ssc.9136.1.S1.at), with a LEO.NB.OCA score of 2.38. Additionally, SGK1, CEBPB, and NFKBIA—all known to be regulated by cortisol—were predicted as downstream of cortisol. Four genes (ATP6V1H, PCTP, DERL1, and TAF1A) of 70 had cis-eQTL (Table 4). In contrast, in muscle only 2 of 593 genes were predicted as upstream (RPS11 and ARRB2), and 25 transcripts were predicted as responsive to cortisol (Table S8). Both genes predicted to affect cortisol concentration (RPS11 and ARRB2) had trans-eQTL, whereas four genes being responsive to cortisol (MPDU1, TAX1BP3, CYBASC3, and POLR2G) revealed cis-eQTL (Table 5). The highest LEO.NB.OCA scores (0.854 and 0.804) were for genes NFKBIA (Ssc.4759.1.S1.at), CXCL3 (Ssc.19692.1.S1.at), and POLR2G (Ssc_1925_1_S1_at), involved in glucocorticoid receptor signaling. Other responsive genes to plasma cortisol concentrations were involved in inflammatory response (RGS2, F13A1, CAPG, and HLA-A) transcription regulation (SERTAD1, ZNF235, BHLHE41, and RUNX1), and cellular growth and proliferation (PSMB2, TAX1BP3, and CAMK2N1) (Table 5).
After screening 6608 markers from eQTL analysis in liver, together with 11 markers from a GWA study, the best set of markers serving as causal anchors for orienting the edges of the network related to cortisol concentrations were DIAS0000681, MARC0052121, and MARC0000128. The best set of markers serving as causal anchors from eQTL analysis in muscle, together with markers from the GWA study, were MARC0001638, DIAS0000681, ASGA0027973, and MARC0091789. We found linkage disequilibrium between MARC0000128 and MARC0001638 with ProbChi = 2.93 × 10−27. A robustness analysis was also applied, for maximum anchoring SNP number ranging from one to five. According to the robustness evaluation for three or more SNPs anchored to each trait, the LEO.NB.OCA scores are generally higher and more stable (Table S9 and Table S10). An outline of the systems biology approach to dissecting the effects of cortisol concentrations in target organs of liver and muscle is shown in Figure 4.
Outline of the systems biology approach to dissecting effects of cortisol concentrations in target organs, liver and muscle. Expression levels in liver and muscle that correlated with cortisol concentrations revealed biological functions. These transcripts were integrated with SNP genotypes to map expression QTL. Significant markers from genome-wide association together with eQTL were merged to prioritize genes using network edge orienting. This approach identified gene networks in muscle or liver that affect or respond to cortisol concentrations.
Discussion
Cortisol is released from the adrenal glands and exerts diverse actions on the immune system, glucose metabolism, and the central nervous system. Few studies have investigated genetic variants in relation to both cortisol secretion and expression pattern in the target organs. In the current study we integrated eQTL and causality modeling to identify candidate genes that affect or respond to cortisol concentrations in the target organs of liver and muscle. The identified candidate genes associated with cortisol concentrations may enhance discovery of genes associated with immune, metabolic, and behavioral traits in humans. Similar studies have identified candidate genes associated with cortisol concentrations that are also causally associated with depression (Gottesman and Gould 2003; Zobel et al. 2010; Velders et al. 2011).
Biological function and molecular pathways correlated with cortisol concentrations
The application of global gene expression analyses to studies of cells and tissues has provided a wealth of data relevant to complex traits. Here, ∼14% of liver transcripts and 9% of muscle transcripts were found to correlate with plasma cortisol concentrations. Finding statistically significant correlations between a trait and particular genes suggests a biological relationship between them (Ponsuksili et al. 2008, 2009, 2011). The potential biological functions and molecular pathways of cortisol-correlated expression differed between liver and muscle. In liver, more genes than expected by chance were identified as involved in lipid and carbohydrate metabolism; in contrast, in muscle, genes encoding various macromolecules related to protein synthesis, cell cycle, and cell death were represented.
Canonical pathway analysis highlighted pathways involving cytokine signaling, including acute phase response signaling, IL-10, and IL-6. The acute phase response (APR) is the body’s innate immune defense against infection, inflammation, stress, physical trauma, and pathological conditions such as cardiovascular disease (Ramji et al. 1993; Baumann and Gauldie 1994). The liver is the major site for production of acute-phase protein (APP), whose expression is regulated at the transcriptional level by cytokines (Baumann and Gauldie 1994; Moshage 1997). We identified 18 molecules in liver and 16 molecules in muscle that were positively correlated to cortisol and are involved in APR signaling. Inflammatory cytokines including IL-1, IL-6, tumor necrosis factor-α (TNF-α), and interferon-γ (IFN-γ) are considered to be the mediators of immunological and pathological responses to stress and infection (Kelley et al. 1994; Warren et al. 1997). In the present study, FKBP5, which encodes a member of the immunophilin family, was significantly correlated with cortisol concentrations (r = 0.23, P = 0.004, FDR = 0.14) but expression levels were not regulated by an eQTL. FKBP5 is a cochaperone of HSP90, which is part of a receptor complex that regulates the sensitivity of the GR. Interestingly, a recent report found that an FKBP5 SNP is associated with a decrease in cortisol and the likelihood of depressive symptoms (Velders et al. 2011). Indeed, several studies have shown the association between high cortisol concentrations and depression (Vreeburg et al. 2009; Velders et al. 2011). Associations between inflammatory markers, the immune system, and depression have also been described (Miller et al. 2009; Miller 2010). In animals, the elevation of APPs by transportation and weaning stress has been reported (Richards et al. 2000; Arthington et al. 2003). Weaning stress may induce an APR through increased cortisol production and modulation of inflammatory cytokines (Kim et al. 2011). The relationship between metabolic syndrome, oxidative stress, and inflammation has been described (Bastard et al. 2006; Holvoet 2008). Similarly, cortisol is a potent glucocorticoid in energy homeostasis in Syrian hamsters, but cortisol alone does not mediate stress-induced increases in food intake or body mass (Solomon et al. 2011). Stress hormones released in response to activation of the HPA axis have been shown to affect aspects of the immune system, the metabolic system, and behavior, including depression. Here, we show for the first time the relationship between systemic plasma cortisol concentration and the transcriptome of two tissues/organs, muscle and liver, that both have implications for organismal metabolism. This study points to complex molecular pathways of target molecules of cortisol in liver and muscle.
Whole-genome association analyses for abundance of transcripts with cortisol concentration-correlated expression (eQTL)
We present the first eQTL study of genes influencing or influenced by plasma cortisol concentrations, using whole-genome association analysis in a commercial crossbred population. In the first step we analyzed the correlation between gene expression in liver and muscle and plasma cortisol concentrations, which revealed biologically meaningful relationships. In the second step eQTL were identified for transcripts that showed trait-correlated expression, which informed us about the genomic location of putative regulatory loci. This strategy eliminated several thousand eQTL that are not associated with cortisol concentrations. For the identification of sequence variants that have cis- and trans-regulatory effects on expression traits, SNPs located within a 10-Mb window at the location of the probe set were defined as cis-eQTL. The trans-acting eQTL represent transcripts whose abundance is regulated by loci remote from the genomic locus of each of these genes. Our goal was to identify a small number of differentially expressed candidates from the thousands of transcripts correlated with cortisol concentration. The eQTL analysis provides evidence for the presence of genomic variation that affects the abundance of transcripts in relation to plasma cortisol concentrations. In fact, for 65% and 62% of the transcripts whose expression in liver and muscle, respectively, correlated with cortisol concentration, significant eQTL were found. Thus, evidence exists for heritable gene regulation.
Candidate gene screening from eQTL and causality modeling
Recently, causality modeling was successfully applied in a variety of different settings (Farber et al. 2009a,b; Plaisier et al. 2009; Park et al. 2011). With the causality modeling algorithms, we integrated GWA and eQTL markers to identify up- and downstream molecules related to cortisol concentration. In liver, plasma cortisol concentration was anchored by three SNPs (DIAS0000681, SSC18, 41.8 Mb; MARC005212, unknown; and MARC0000128, SSC9, 20.5 Mb) that explained ∼3–4% of the variation in cortisol concentrations (n = 475). MARC0000128 is located with the gene LOC100627939 (disks large homolog 2-like, DLG2). Recently, a genome-wide association study in human revealed significant association of DLG2 with circulating phospholipid levels that play a crucial function in several diseases, with diverse neurological, psychiatric, and metabolic consequences (Demirkan et al. 2012). In the neighborhood (±10,000 bp) of the other markers there are no currently annotated genes in the porcine genome. Additionally, 63 known genes were regulated by cortisol and anchored by markers, most of which corresponded to eQTL markers. These eQTL markers explained, on average, >15% of the variation in expression levels. Many genes are regulated by corticosterone (e.g., SGK1, CEBPB, and NFKBIA), some of which were also predicted as downstream of cortisol in this study. For instance, SGK1, whose expression is directly affected by corticosterone with an emphasis on metabolism, is involved in cellular functions such as glucose transport (Jeyaraj et al. 2007), regulation of voltage-gated channels (Boehmer et al. 2008), and apoptotic processes (Schoenebeck et al. 2005). Other molecules in our study that were regulated by cortisol with an emphasis effect on lipid or transport metabolism were ACSL1, APOC3, PCTP, SREBF2, INSIG1, LIPG, NFKBIA, and ITGB3. Ten of 63 known genes (CEBPB, CCRN4L, TAF1A, EPAS1, CEBPD, YY1, TEF, AKIRIN2, NFIB, and SREBF2) that were reactive to cortisol have transcription factor activity. CEBPB and CEBPD are transcription factors belonging to the C/EBP family. These contain a conserved basic leucine zipper motif at the carboxy terminus that is involved in dimerization and DNA binding (Ramji and Foka 2002). Binding sites for the C/EBPs have been identified in the regulatory regions of a large number of APP genes. Further, the main stimulators of APP production are inflammation-associated cytokines such as IL-6, IL-1, tumor necrosis factor-α, and IFN-γ. Three genes (CEBPB, TLR4, and AKIRIN2) downstream of cortisol are involved in IL-6 regulation. Toll-like receptor 4 (TLR4), a key member of the TLR family, has been well characterized as inducing inflammatory products of innate immunity and stress response (Eicher et al. 2004; O’Loughlin et al. 2011; Wang et al. 2011).
While most studies have covered genes downstream of cortisol, we identified 24 known genes that are upstream of cortisol concentrations. These 24 genes can be categorized into Gene Ontology (GO) terms of transcription regulator activity (TCF21, IKBKAP, NFRKB, and ETV3), intracellular signaling cascade (IKBKAP, MAP3K3, PREX1, STK4, and PMEPA1), or KEGG_PATHWAY of MAPK signaling pathway (PLA2G10, MAP3K3, and STK4). The most significant transcript influencing cortisol concentrations was CD99. CD99 is a leukocyte surface glycoprotein and plays a key role in several biological processes, including the regulation of T-cell activation (Pata et al. 2011), as well as in pig-to-human xenotransplantation (Schneider et al. 2009). Another gene that was identified as contributory for cortisol concentrations is phospholipase A2, group X (PLA2G10). PLA2G10 has been implicated in diverse biological functions. One study showed that lungs from mice deficient in GX sPLA2 have significantly reduced ovalbumin-induced inflammatory responses, suggesting a role for GX sPLA2 in asthma (Henderson et al. 2007). Additionally, macrophage secretory phospholipase A2 group X enhances anti-inflammatory responses, promotes lipid accumulation, and contributes to aberrant lung pathology (Curfs et al. 2008). Recently, a study showed that C57BL/6 mice with targeted deletion of group X secretory phospholipase A2 (GX KO) have ∼80% higher plasma corticosterone concentrations compared with wild-type (WT) mice under both basal and adrenocorticotropic hormone (ACTH)-induced stress conditions (Shridas et al. 2010). Further, GX sPLA2 negatively regulates corticosterone production by altering the expression of steroidogenic acute regulatory protein (StAR), likely by suppressing LXR-mediated StAR promoter activation. Here, we confirmed that expression levels of PLA2G10 negatively correlate with cortisol concentrations (r = −0.243, P = 0.003, FDR = 0.12). Indeed, we found that it is regulated by a cis-eQTL (NegLog10 of P = 5.52) and that it affects cortisol concentration by causality (LEO.NB.OCA = 0.74).
In muscle, cortisol concentrations were anchored by two significant SNPs (MARC0001638, SSC9, 20.8 Mb; DIAS0000681, SSC9, 20.5 Mb). The results correspond well with those obtained for liver: Marker DIAS0000681 was also anchored to cortisol concentration in liver, and MARC0001638 is in linkage disequilibrium to MARC0000128, which anchored cortisol in liver. As expected, most of the causal relationships were “cortisol modulates gene expression.” One of the two transcripts identified in muscle as upstream of cortisol concentrations was ARRB2. Recently, a study showed that disruption of β-arrestins (Arrb1−/−, Arrb2−/−) blocks GR in mice, suggesting that GR is an important downstream effector of the β-arrestin signaling pathway involved in regulation of lung and liver development (Zhang et al. 2011). Another study reported that peptidoglycan-induced immune protection in stress-induced changes of cytokine levels appears to require β-arrestin 2, a multifunctional adaptor and signal transducer (Li et al. 2011). Here, we show for the first time the up- and downstream effects on cortisol concentrations in muscle. As in liver, many molecules are involved in inflammatory response (ARRB2, CAPG, CXCL3, HLA-A, NFKBIA, RGS2, and F13A1) and glucocorticoid receptor signaling (ARRB2, NFKBIA, POLR2G, and CXCL3). Other molecules we identified were noted in the literature as related to stress and therefore indirectly related to cortisol. Polymorphisms in one candidate, RGS2 (regulator of G-protein signaling 2), are associated with anxious behavior in mice (Leygraf et al. 2006) and anxiety in humans (Smoller et al. 2008).
Causality modeling has been successfully applied in our study to understand the role of cortisol as either a downstream target or an effector molecule. This approach allowed us to identify two different types of transcripts and associated markers: those influencing expression of individual genes that, in turn, influence cortisol concentrations and, alternately, those influencing cortisol concentrations that, in turn, influence gene expression. We thus identified a list of genes with upstream relationships to plasma cortisol concentrations, as well as a list of genes that are regulated downstream by cortisol concentrations in liver and muscle. Interestingly, the list of our candidate genes in this study did not overlap with that of a previous study by using a comparative, translational systems genetics approach addressing cortisol in mice and pigs (Kadarmideen and Janss (2007). Whereas the previous study of Kadarmideen and Janss (2007) focused on gene expression from mouse INIA brain mRNA, our study aimed to discriminate up- and downstream effects of transcripts affecting or responding to plasma cortisol concentrations in the target tissue of liver and muscle.
In summary, looking at trait-correlated expression patterns in genes in two target tissues of cortisol—liver, the central metabolic organ and source of bioactive molecules, and muscle, the major energy consumer and storage tissue—has helped to elucidate the complex molecular networks contributing to cortisol concentrations in plasma, which, in turn, relates to behavioral, physiological, and metabolic disease in mammals. Many of the candidate genes previously identified in humans and mice have been confirmed here in the porcine model. Moreover, novel candidate genes associated with cortisol concentrations, identified for the first time in a pig model, are valuable for human research due to the genomic and physiological similarities between pigs and humans. The application of this information can facilitate the detection of trait-associated SNPs related to behavioral, physiological, and metabolic traits. The discrimination of up- and downstream effects of transcripts affecting or responding to plasma cortisol concentrations further supports gaining a better understanding of the biology of complex traits.
Acknowledgments
The authors thank Annette Jugert, Joana Bittner, and Hannelore Tychsen for excellent technical help. This research was partly supported by the German Research Foundation (Deutsche Forschungsgemeinschaft grants FOR 753, PO 753/7-1, and WI 1754/10-14) and by matched funding from the FBN.
Footnotes
Communicating editor: F. F. Pardo Manuel de Villena
- Received June 19, 2012.
- Accepted August 6, 2012.
- Copyright © 2012 by the Genetics Society of America