eQTLs identify regulatory networks and drivers of variation in the individual response to sepsis

Abstract

Sepsis is a clinical syndrome of life-threatening organ dysfunction caused by a dysregulated response to infection, for which disease heterogeneity is a major obstacle to developing targeted treatments. We have previously identified gene-expression-based patient subgroups (sepsis response signatures [SRS]) informative for outcome and underlying pathophysiology. Here, we aimed to investigate the role of genetic variation in determining the host transcriptomic response and to delineate regulatory networks underlying SRS. Using genotyping and RNA-sequencing data on 638 adult sepsis patients, we report 16,049 independent expression (eQTLs) and 32 co-expression module (modQTLs) quantitative trait loci in this disease context. We identified significant interactions between SRS and genotype for 1,578 SNP-gene pairs and combined transcription factor (TF) binding site information (SNP2TFBS) and predicted regulon activity (DoRothEA) to identify candidate upstream regulators. Overall, these approaches identified putative mechanistic links between host genetic variation, cell subtypes, and the individual transcriptomic response to infection.

Keywords: co-expression; eQTL; genomics; infection; sepsis; transcriptomics.

Graphical abstract

Figure 1

Genetic variants associated with gene expression in the context of sepsis (A) Schematic of cohort design for SRS1ever genome-wide association study (GWAS) using all patients with genotyping data and at least one gene-expression time point for SRS assignment. (B) Common SNPs (MAF ≥ 1%) were tested for association with the SRS1ever vs. never phenotype. Manhattan plot showing −log10(p value) for each variant plotted against its genomic position. The most significant SNP in each locus is highlighted in orange. (C) Schematic of cohort design for eQTL (cis-eQTL and co-expression module QTL) analysis using all samples from patients with genotyping data and RNA-seq data. Co-expression modules were defined using the full RNA-seq dataset. (D) Histogram showing the distribution of the numbers of independent signals detected through conditional analysis for each eGene. (E) eQTL interactions with source of sepsis (CAP or FP). Each point represents an independent eSNP-eGene pair, with the interaction effect size plotted against the genotype effect. eQTLs with bigger effects in FP compared to CAP are therefore found in the top right and bottom left quadrants. Red indicates a significant interaction between genotype and source of sepsis (FDR < 0.05), with the most significant results labeled with the eGene name. (F) Sepsis-dependent eQTL effects identified with mashr. Each point represents a lead SNP-eGene pair from the first-pass eQTL mapping in sepsis patients that was also tested for whole-blood eQTL in the European subset of GTEx. Posterior effect sizes estimated by mashr are plotted for GTEx against sepsis, and eQTLs are categorized based on the difference between these estimates. eQTLs significant in sepsis are "shared" if the mashr posterior effect size is in the same direction as and within a factor of 0.5 of the GTEx effect size. Those with bigger effects in the same direction in sepsis or GTEx are “sepsis-magnified” and “sepsis-dampened,” respectively. Those significant in both GAinS and GTEx but with opposite directions of effects are “opposite direction of effect” (n = 53). Those significant only in sepsis are also classed as “sepsis-magnified,” and those significant in neither cohort are “not significant.” Please see also Figures S1–S9 and Tables S2, S3, S4, S5, S6, S7, S8, S9, and S10.

Figure 2

Genotype-by-environment interactions find widespread variation in eQTL effects across sepsis patients (A) eQTL interactions with SRS1 status. Each point represents an independent eSNP-eGene pair, with SRS interaction effect size plotted against the genotype effect. eQTLs with bigger effects in SRS1 compared to non-SRS1 are therefore found in the top right and bottom left quadrants. Red indicates a significant interaction between genotype and SRS1 status (FDR < 0.05), with the most significant results labeled with the eGene name. (B) An exemplar sepsis-magnified eQTL that also has a significant positive interaction with SRS1 status. Gene expression residuals were modeled with an SRSq-by-genotype interaction to illustrate the continuous relationship of SRS status with the genotype effect, with point color indicating number of copies of the minor allele of rs4378192. (C) UpSet plot showing sharing of magnifying (red) and dampening (blue) eQTL interaction effects between environmental variables tested. (D) Enrichment of eSNPs for eQTLs with an SRS interaction in different chromatin states from the Roadmap Epigenomics core 15-state genome segmentation annotations across relevant cell types, compared to eSNPs for eQTLs without a significant SRS interaction (only significant enrichments shown [FDR < 0.05]). Please see also Figures S10 and S11 and Tables S11 and S12.

Figure 3

Identification of putative driver transcription factors for SRS from eQTL interactions (A) Schematic of strategy to identify transcription factor (TF) binding sites (TFBSs) enriched in SRS interaction eQTLs. For each of 12,959 eQTL signals tested for an SRS interaction effect, query SNPs were defined as those in LD (r² ≥ 0.8) with the lead SNP. We then identified instances where binding motifs for 124 human TFs were interrupted or introduced by these query SNPs using SNP2TFBS. Each independent eQTL signal was scored for each motif as having ≥1 or 0 binding sites altered by the signal SNP or its LD proxies. We then tested for enrichment of each TF motif among eQTLs with a significant interaction effect compared to eQTLs with no significant interaction using a one-tailed Fisher’s exact test. (B) TF binding sites identified with SNP2TFBS enriched in eQTLs with an SRS interaction vs. no interaction, with point color indicating significance. (C) Volcano plot showing comparison of inferred TF activity between SRS1 and non-SRS1 samples. Each point represents one TF, with adjusted p value plotted against effect size estimated from a linear mixed model. Red indicates significance. (D) Adjusted p values from (B) are plotted against those from (C) to highlight TFs with both enriched binding sites in SRS interaction eQTLs and differential activity between SRSs. (E) Heatmap showing Spearman correlation between estimated cell proportions and inferred TF activity for the first sample available per patient. White indicates a non-significant correlation (FDR > 0.05), red a positive correlation, and blue a negative correlation. Please see also Figures S12 and S13 and Tables S13, S14, and S15.

Figure 4

Co-expression modules pinpoint several trans-regulatory networks relevant to sepsis outcomes (A) Heatmap showing enrichment of sepsis cell markers in module member genes. Modules were tested for enrichment of leukocyte marker genes identified in a sepsis cohort. Modules shown had significant enrichment for at least one signature. (B) Heatmap showing significant associations between module eigengenes (MEs) with a modQTL and clinical phenotypes. modQTLs passing stringent sensitivity analysis are marked with bold font. MEs were tested for differential expression with measured cell proportions, SRS1 status, diagnosis (CAP or FP), and time point (day 1, 3, or 5) using a linear mixed model. Association of each ME with survival up to 28 days was tested using a Cox proportional hazards model. (C) Module 92 eigengene (ME_92) plotted by rs821470 genotype, with partial residuals calculated from the linear-mixed-model fit. (D) Circos plot showing the chromosomal locations of the genes contained in module 92 and the lead eSNPs associated with the ME. Member genes that are cis-eGenes for these eSNPs are highlighted in orange. (E) ME_47 plotted by rs16843927 genotype, with partial residuals calculated using the linear-mixed-model fit. (F) Circos plot showing the chromosomal locations of the genes contained in module 47 and the lead eSNPs associated with the ME. Genes that are cis-eGenes for these eSNPs are highlighted in orange. Please see also Figures S14–S20 and Tables S16, S17, S18, S19, S20, S21, S22, S23, and S24.