STARNET is a cohort-based study of 600 individuals undergoing coronary artery bypass grafting (CABG) for coronary artery disease (CAD) and was used as the primary discovery cohort in this study. These individuals underwent blood genotyping preoperatively for 951,117 genomic markers, and during surgery, seven different tissue samples were obtained and underwent RNA-sequencing (RNA-seq): liver, skeletal muscle, atherosclerotic aortic root, internal mammary artery, visceral abdominal fat, subcutaneous fat, and whole blood. STARNET data are available through a database of Genotypes and Phenotypes (dbGaP) application (accession no. phs001203.v2.p1). A detailed description of data processing can be found in the Supplemental Material of this article (section S1.1).

The Stockholm Atherosclerosis Gene Expression (STAGE) study (n = 114) (20) and the Metabolic Syndrome in Man (METSIM) study (n = 982) (21) were used in the replication of causal gene networks identified using STARNET. Gene expression data for the METSIM and STAGE studies are available publicly at Gene expression omnibus (GEO) (accession no. GSE70353 and GSE40231, respectively). Microarray data for the liver, subcutaneous fat, and visceral abdominal fat were used from the STAGE study, and gene expression data from subcutaneous fat were measured in the METSIM study using RNA-seq.

A list of SNPs associated with plasma cortisol was obtained from the summary statistics of the 2021 GWAMA conducted by the CORNET consortium (available at https://datashare.ed.ac.uk/handle/10283/3836) (8). We filtered this list to obtain SNPs that were found to be associated with plasma cortisol at a level of genome-wide significance (p< 5 × 10^-8) that were taken forward 68 and tested against all genes across STARNET tissues.

The secondary linkage test (P2) is a likelihood ratio test in the Findr package (22) (version 1.0.8) that was used to identify associations between a given SNP (E) and a gene (B) using categorical regression. P2 proposes a null hypothesis where E and B are independent and alternative hypotheses where E is causal for B (E →B). Maximum likelihood estimators are then used to obtain a log likelihood ratio (LLR) between the alternative and null hypotheses. The LLR is then converted to the posterior probability of the alternative hypothesis ℋ a l t ( P 2 ) being true with empirical estimation of the local false discovery rate (FDR) as a value from 0 to 1 (Equation 1).

Multiple datasets were used to identify genes that had prior evidence of putative regulation by GR (23–27). These datasets have been filtered to include targets for NR3C1, the gene that encodes GR.

Trans-genes were categorized according to evidence of GR regulation from datasets shown in Supplementary Table S1 . Genes were scored against these criteria: 1) appearing in a transcription factor database (ENCODE, TRANSFAC, CHEA); 2) identified as a GR target from chromatin immunoprecipitation sequencing (ChIP-seq) experiment in adipocytes from Yu et al. (23); 3) differentially expressed in response to dexamethasone treatment in adipocytes from Yu et al. (23); and 4) murine homolog of human gene differentially expressed in response to dexamethasone treatment using adrenalectomized mice (FC >1; p-value<0.05) (24). Genes were then ranked according to how well they met the criteria for GR regulation (+1 for each item matched from criteria 1–4).

Pairwise causal inference was used for the reconstruction of cortisol-responsive transcriptional networks across STARNET tissues using cis-eQTL genotypes as genetic instruments with gene expression data from STARNET, as implemented by the Findr software (22). A detailed description of these methods can be found in the Supplementary Material of this article (Section S1.2).

Lists of known transcription factor targets for both NR3C1 and IRF2 were obtained from ENCODE and TRANSFAC datasets, respectively. These datasets were used to test for an enrichment of known transcription factor targets within novel gene sets derived from gene network targets. This was performed using Fisher’s exact test from the Python module Scipy Stats (28) and involved the creation of a 2 × 2 contingency table based on a tissue-specific background consisting of all genes available in the corresponding tissue.

Correlations between gene network targets were calculated using gene expression data from STARNET, STAGE, and METSIM. Gene expression matrices were filtered to only include the target genes under investigation. Correlation matrices of corresponding Pearson correlation coefficients as absolute values were constructed in Python.

A background gene set was constructed from the overlapping genes between the STARNET gene expression set that was used for network discovery and the corresponding gene expression set that was being used for replication. The previously described correlation analysis was then repeated using a random set of genes (the same size as the target set) selected from the background gene set. The Kruskal–Wallis test was implemented in Python using Scipy Stats (28) to test if the targeted and randomly sampled correlations follow the same distribution. Both the targeted and random correlations were then plotted as a boxplot using the Python plotting package Seaborn (29).

Hierarchical clustering was performed on correlation values between network targets using the discovery (STARNET) gene expression data and hierarchical clustering from Scipy Stats (28) in Python. The leaves list that resulted from the clustering of the discovery dataset was then extracted and applied to the correlations between target genes from the corresponding replication dataset. Both sets of clustered correlation values were then plotted as opposing correlation heatmaps with Seaborn (29).

SNPs associated with plasma cortisol at the SERPINA6/SERPINA1 locus have previously been linked as expression single-nucleotide polymorphisms (eSNPs) for SERPINA6 in the liver (8). Using genotype and tissue-specific RNA-seq data from the STARNET cohort, we explored the hepatic and extrahepatic consequences of genetic variation for plasma cortisol using 73 cortisol-associated SNPs at genome-wide significance (p< 5 × 10^-8) identified from the CORNET GWAMA (8). We identified 704 eQTL associations in cis and trans between plasma cortisol-associated SNPs and genes measured across all STARNET tissues, composed of 262 unique genes and 72 SNPs at a 15% FDR threshold ( Supplementary Tables S2 , S3 ).

The tissues with the greatest number of trans-genes were the liver, subcutaneous fat, and visceral abdominal fat, with a combined total of 157 trans-genes and 422 total SNP–gene associations (FDR = 15%) ( Figure 1A ). The vast majority of trans-eQTL associations were specific to a single tissue. A single trans-gene, the glycosyltransferase-encoding gene OGT, was identified in both the liver and visceral abdominal fat. However, as this was the only cross-tissue trans-gene identified, suggesting that the transcriptional impact of genetic variation at the SERPINA6/SERPINA1 locus is highly tissue-specific. The CORNET GWAMA describes four blocks of SNPs in linkage disequilibrium (LD), which represent the cortisol-associated variation at the SERPINA6/SERPINA1 locus (8). We observed that LD blocks 2 and 4 represent the majority of the variation across all tissues in the trans-gene sets ( Figures 1B, C ).

As the GR is the primary mechanism by which cortisol influences transcription, we sought to identify a subset of cortisol-associated trans-genes that were also regulated by the GR. The cortisol-associated trans-genes identified in this study were compared to sets of known GR targets identified from different sources as described in Supplementary Table S1 . This included large projects such as ENCODE, TRANSFAC, and CHEA that predict transcription factor-binding targets from high-throughput transcription factor-binding assays. We also included predicted GR targets from perturbation-based experiments in specific tissues. ChIP-seq and microarray analysis has been used to identify 274 glucocorticoid-regulated genes in 3TS-L1 adipocytes, a murine-derived cell line (23). In addition, RNA-seq data in subcutaneous fat from adrenalectomized mice treated with dexamethasone, a GR agonist, have been used to identify genes that are differentially expressed (24).

The greatest number of unique cortisol-associated trans-genes was identified in the liver (n = 43), subcutaneous fat (n = 54), and visceral abdominal fat (n = 59) at a 15% FDR threshold. The involvement of these tissues in glucocorticoid signaling and physiological effects has been well documented in the literature (31–34); therefore, the identification of GR-regulated trans-genes was restricted to these tissues. Comparisons of genes identified as glucocorticoid-regulated in 3T3-L1 adipocytes were only made with subcutaneous and visceral adipose trans-genes. Likewise, as the murine RNA-seq experiments were restricted to subcutaneous adipose, only subcutaneous adipose trans-genes were compared to these differentially expressed genes.

In the liver trans-gene set, 19/43 genes were identified that were present in either the ENCODE, TRANSFAC, or CHEA datasets (FDR = 15%) ( Figure 1D ; Supplementary Table S4 ). This includes SERPINA6 that is cis-associated with genetic variation for plasma cortisol, as described previously (8). One gene, CPEB2, was identified in more than one dataset and was present in both ENCODE and CHEA. CPEB2 (posterior probability = 0.89) is a regulator of translation, splice variants of which have been linked to cancer metastasis (35).

Visceral adipose tissue had the largest number of cortisol-associated trans-genes. Here, 21/59 of these genes had some evidence of being targets of GR ( Figure 1E ; Supplementary Table S5 ). There were five genes that had been identified as GR targets from both high-throughput transcription factor-binding assays and adipose-specific experiments. These include CD163 and LUC7L3. CD163 is a hemoglobin scavenger protein that is expressed in macrophages and involved in the clearance of hemoglobin/haptoglobin complexes that may play a role in the protection from oxidative damage. It also plays a role in activating macrophages as part of the inflammatory response (36). LUC7L3, also known as CROP, encodes a protein that is involved in alternative splicing and is associated with human heart failure (37). It has also been shown to play a role in the inhibition of hepatitis B replication (38).

Of the cortisol-associated trans-genes identified in subcutaneous adipose (FDR = 15%), 28/54 genes were either present in a transcription factor dataset or identified from the adipose-specific perturbation datasets ( Figure 1F ; Supplementary Table S6 ). There were 13 genes that had been identified as GR targets from both high-throughput transcription factor-binding assays and adipose-specific experiments. These include RNF13 that encodes IRE1α-interacting protein that plays an important role in the endoplasmic reticulum (ER) stress response through regulation of IRE1α, a critical sensor of unfolded proteins (39). Also IRF2, encoding the transcription factor Interferon Regulatory Factor 2 that plays an important role as a repressor of IRF1 that in turn is involved in the interferon-mediated immune response (40). Furthermore, IRF1 has previously been identified as a marker for glucocorticoid sensitivity in peripheral blood (41).

Having identified cortisol-associated trans-genes that are regulated by GR, causal estimates were obtained for pairwise relationships between GR-regulated trans-genes and all other genes within the given tissue. This was carried out for all GR-regulated trans-genes in the liver, subcutaneous fat, and visceral abdominal fat with a valid cis-eQTL instrument (12, 19, and 7 genes, respectively) ( Supplementary Table S7 ). A 10% global FDR threshold was then imposed for each gene set ( Table 1 ). Primary networks were obtained by filtering to include only GR trans-genes with a minimum of four target genes at the global FDR threshold.

In the liver, we identified a single gene network driven by CPEB2, which was found to be trans-associated with the cortisol-associated SNP rs4905194 ( Figure 2A ). This network contained 48 causal interactions driven by CPEB2 at a 10% FDR threshold ( Figure 2D ; Supplementary Table S9 ). It is notable that CPEB2 appears as the only network regulator in the liver considering it was also the cortisol-associated trans-gene with the strongest links to GR regulation from the liver trans-gene set. A detailed description of the CPEB2 network and all other networks identified can be found in the Supplementary Information (Section S2.1).

In subcutaneous fat, two major subnetworks were identified under the regulation of the genes RNF13 and IRF2. This includes a total of 343 causal relationships across both subnetworks, including two genes shared by both subnetworks. RNF13 was found to be trans-associated with the cortisol-associated SNP rs11622665 ( Figure 2B ) and represents the largest subcutaneous fat subnetwork with 215 gene targets at a 10% FDR threshold ( Figure 2E ; Supplementary Table S10 ).

The transcription factor IRF2, which was associated with the cortisol-linked SNP rs8022616 ( Figure 2C ), was found to putatively regulate a network of 128 genes (FDR = 10%) ( Figure 2F ). Some notable targets of IRF2 include LDB2 (posterior probability = 0.94) and LIPA (posterior probability = 0.91). GWAS suggests functions for LIPA related to CAD and ischemic cardiomyopathy (42), while LDB2 has been demonstrated to be involved in the development of atherosclerosis (43). Additionally, cortisol has been shown to induce a 5-fold reduction in LDB2 expression in adipocytes (44).

Predicted IRF2 transcription factor targets have been previously described as part of the TRANSFAC dataset. We examined the overlap between predicted IRF2 targets in TRANSFAC, and gene targets within the IRF2 causal networks were identified in subcutaneous fat. A true network of IRF2 targets would be expected to show an enrichment of predicted IRF2. Using Fisher’s exact test on data from subcutaneous fat, at a 10% FDR threshold, the IRF2 network had 128 target genes, 35 of which were also predicted IRF2 targets (p = 0.08); at a 15% FDR threshold, 104/247 causal targets were also predicted targets of IRF2 in TRANSFAC (p = 0.005). Decreasing the global FDR beyond this threshold increased the number of TRANSFAC targets within the pool of causal targets, however at a lower enrichment (p = 0.046) ( Supplementary Table S12 ).

In addition to examining the prevalence of IRF2 targets within the IRF2 causal network, we investigated the overlap between network genes that are also regulated by GR. We observed an enrichment of ENCODE GR targets at 15% and 20% FDR thresholds (p< 0.05) including 68 and 138 GR targets, respectively. No GR enrichment was observed in either CHEA or TRANSFAC datasets for IRF2 networks.

Causal gene networks represent coordinated changes in gene expression in response to changes in the expression of network regulators. Therefore, it is possible to examine if these changes in gene expression are present in independent datasets using gene expression data alone. We used RNA-seq and microarray data from the METSIM and STAGE datasets, respectively, to compare patterns in gene expression within causal networks predicted from STARNET. As METSIM only contains gene expression data for subcutaneous fat, analysis was restricted to the causal networks identified in STARNET subcutaneous fat.

Absolute correlation coefficients between the targets of the previously described network regulators were calculated, and their distributions were compared to distributions of random sets of genes selected from the replication gene expression data, the same size as the corresponding target gene set. The difference between targeted and random distributions was formalized using the Kruskal–Wallis test for each subnetwork ( Table 2 ).

In the liver, correlations between network targets of the single subnetwork under the regulation of CPEB2 were observed in STARNET and STAGE. Hierarchical clustering within the STARNET liver also revealed clustering of correlated genes that were retained when the clustered gene order was then applied to the STAGE liver ( Figure 3A ). Correlations between the 44 CPEB2 target genes in the STAGE liver were stronger than their random counterparts (p = 8.2 × 10^-32), with this shift also being observed in the STARNET liver (p = 2.32 × 10^-197) ( Figure 3D ).

In subcutaneous fat, correlations were observed between the network targets of RNF13 and IRF2, and hierarchical clustering patterns from STARNET were applied to the replication datasets of STAGE and METSIM ( Figures 3B, C ). For RNF13, similar patterns of co-expression were observed in the STAGE subcutaneous fat following clustering; however, this was not the case in the METSIM dataset ( Figure 3B ). Despite this, RNF13 targets appeared more highly correlated than their randomly selected counterparts in STARNET (p< 1.0 × 10^-300), STAGE (p< 1.0 × 10^-300) and to a lesser extent in METSIM (p = 2.3 × 10^-7) ( Figure 3E ).

In subcutaneous fat, patterns of co-expression between IRF2 targets were conserved most prominently in METSIM; however, co-expression was less strongly correlated compared with RNF13 targets ( Figure 3C ). IRF2 subcutaneous fat subnetwork targets were more strongly correlated than their random counterparts in STARNET (p< 1.0 × 10^-300), STAGE (p = 8.35 × 10^-86), and METSIM (p< 1.0 × 10^-300) ( Figure 3F ).