Genetic Factors of the Disease Course after Sepsis: A Genome-Wide Study for 28Day Mortality.

Sepsis is the dysregulated host response to an infection which leads to life-threatening organ dysfunction that varies by host genomic factors. We conducted a genome-wide association study (GWAS) in 740 adult septic patients and focused on 28day mortality as outcome. Variants with suggestive evidence for an association (p≤10-5) were validated in two additional GWA studies (n=3470) and gene coding regions related to the variants were assessed in an independent exome sequencing study (n=74). In the discovery GWAS, we identified 243 autosomal variants which clustered in 14 loci (p≤10-5). The best association signal (rs117983287; p=8.16×10-8) was observed for a missense variant located at chromosome 9q21.2 in the VPS13A gene. VPS13A was further supported by additional GWAS (p=0.03) and sequencing data (p=0.04). Furthermore, CRISPLD2 (p=5.99×10-6) and a region on chromosome 13q21.33 (p=3.34×10-7) were supported by both our data and external biological evidence. We found 14 loci with suggestive evidence for an association with 28day mortality and found supportive, converging evidence for three of them in independent data sets. Elucidating the underlying biological mechanisms of VPS13A, CRISPLD2, and the chromosome 13 locus should be a focus of future research activities.

Introduction

Sepsis is the dysregulated host response to an infection which leads to life-threatening organ dysfunction according to the new Sepsis-3 definition (Singer et al., 2016, Seymour et al., 2016). It can result in 28 day mortalities of up to 60% (Engel et al., 2007, Angus and Wax, 2001). Consequently, there is an urgent need for new therapies but results from recent large scale phase III randomized controlled intervention trials (e.g Food and Drug Administration, 2011) have been disappointing. It has been proposed to go “back to the drawing board” (Angus, 2011) taking a fresh look at the biology that drives the sepsis processes (Cohen et al., 2015).

As part of this discussion, there is new interest in host genomic factors that are rooted in the landmark publication by Sørensen et al. (1988). These authors reported that if one biological parent died of an infection, the risk to die of an infection in the offspring was strongly increased (relative risk 4.52). This work stimulated the conduct of many candidate gene association studies for sepsis susceptibility with inconsistent and essentially weak results (e.g. reviewed in Clark and Baudouin, 2006). Moreover, focusing on sepsis susceptibility might be too challenging given that recent evidence strongly supported a stronger impact of the host genome to account for the variability during the clinical disease course after sepsis onset (Petersen et al., 2010). Thus, this and an accompanying report by Taudien et al. 2016 focus on host genomic factors related to differential clinical disease course after sepsis onset applying the new Sepsis-3 definition. While Taudien et al. 2016 report on deleterious single nucleotide variants and pathways, we describe a genome-wide association (GWA) study (GWAS) which by design is limited to common variants.

Of the two GWA studies related to sepsis reported so far (Man et al., 2013, Rautanen et al., 2015) the former focused on treatment response in 1446 patients with (severe) sepsis while the latter was aiming on 28 day mortality in 1533 patients with sepsis due to pneumonia. Both GWAS used the consensus definition of sepsis from 2001 which did not require the presence of an organ dysfunction (Levy et al., 2003) and only Rautanen et al. (2015) consider host genomic factors related to differential clinical disease course after sepsis onset. As their main finding, Rautanen et al. (2015) report that a common genetic variation in the FER (FER tyrosine kinase) gene is associated with a reduced 28 day mortality from sepsis due to pneumonia. They estimate an age-adjusted odds ratio (OR) of 0.56 (95% confidence interval (CI) [0.45–0.69]; p = 5.6 × 10− 8) for each C allele at rs4957796 in a joint analysis of discovery and replication samples (total 2078 patients).

Here we report results derived under a similar study design focusing on 28 day mortality in a discovery GWAS of 740 septic patients. We follow-up our best GWAS loci with single nucleotide polymorphism (SNP) allelic association signals below the significance level (p < 10− 5), i.e. suggestive evidence for an association, in the discovery meta-analysis by Rautanen et al. (2015) with 2534 patients with sepsis due to pneumonia or abdominal infections combined and in another independent GWAS of the PROGRESS consortium with 936 patients with confirmed community acquired pneumonia (CAP) – both with mortality outcome data. Next, we elucidate the potential differential organ impact of these variants by analyzing organ dysfunction scores after sepsis onset. Finally, we follow-up the loci with the most significant results previously identified by Rautanen et al. (2015) and all 21 candidate genes at or around our best GWAS loci in an independent exome sequencing study (Taudien et al. 2016 that included 74 patients with treated sepsis and 28 day mortality outcome data.

Material & Methods

Study Design and Patients

Discovery GWAS

Our discovery GWAS included patients that participated in two randomized controlled trials (RCTs) VISEP and MAXSEP of the SepNet Study group (Brunkhorst et al., 2008, Brunkhorst et al., 2012). Both RCTs ascertained patients of European ancestry who were admitted to German intensive care units (ICUs) with a diagnosis of sepsis (see Appendix for definitions). For VISEP, patients were recruited at 18 academic tertiary hospitals in Germany between 04/2003 and 06/2005 (n = 537). For MAXSEP, patients were recruited at 44 ICUs in Germany between 10/2007 and 03/2010 (n = 600). Here we analyzed a subgroup of patients from the two RCTs who gave additional written consent to participate in a genetic study and who met patient-wise quality control criteria (nVISEP = 410; nMAXSEP = 330). We included all 740 patients irrespective of treatment group but performed sensitivity analyses to address potential effects of study arm. Supplementary Fig. 16 shows the amount of organ dysfunction among (28 day) survivors and non-survivors based on SOFA (sub-)scores.

Validation GWA Studies

(1) We contacted Rautanen et al. (2015) who looked-up our best 14 GWAS hits in their meta-analysis of three discovery GWAS cohorts (GenOSept/GAinS; VASST; PROWESS) that included up to 2534 patients with sepsis and information on the 28 day mortality outcome. For details on the cohort descriptions and the quality control we refer to the original report (Rautanen et al., 2015). (2) In addition, we looked-up our best 14 GWAS hits in a GWAS of patients from the PROGRESS study. PROGRESS is a prospective multi-centric longitudinal observational study on patients hospitalized due to confirmed CAP. Patients were investigated for five consecutive days after enrolment including comprehensive clinical and laboratory assessments. Vital status was assessed at days 28, 180, and 360 after enrolment. PROGRESS is registered at ClinicalTrials.gov (registration number: NCT02782013).

Exome Sequencing Study

To further follow-up our findings, we performed a moderate-size whole-exome sequencing study in an independent cohort of 74 patients with treated sepsis again with European background which were recruited at two University hospitals (n = 15 at the Jena University Hospital, Germany and n = 59 at the University Hospital Athens, Greece). Sepsis patients for this study were selected for extremely different clinical disease courses - patients with co-morbidities who survived despite an inappropriate empirically administered antimicrobial treatment until the antibiogram became known (n = 37) vs. younger patients with a lack of comorbidities who had a bad disease course (as documented by SOFA trajectories) or died early in the presence of appropriate initial treatment (n = 37). A detailed characterization of all patients is provided in Taudien et al. 2016.

Ethics approval was granted for the individual centers and the study was conducted according to the ethical standards laid down in the Declaration of Helsinki. Written, informed consent was obtained from all patients or from a legal representative in case of critical illness. Table 1 shows patient characteristics of the analyzed patients in the discovery GWAS, the validation GWAS (PROGRESS) and the exome sequencing studies. Details on the validation GWA studies (GenOSept/GAinS; VASST; PROWESS) are provided in Rautanen et al. (2015).

Table 1

Characteristics of the patients with treated severe sepsis/septic shock in the discovery GWAS, the independent (unpublished) validation GWAS in the PROGRESS study and the exome sequencing study (Taudien et al., in press). For the cohort descriptions of the other validation GWA studies, we refer to the original report by Rautanen et al. (2015).

	Discovery GWAS (n = 740)	Validation GWAS PROGRESS (n = 936)	Exome sequencing study (n = 74)
Deaths (or qualified intensive carea) within 28 days (%)	149 (20)	95 (10)a	12 (15)
Females (%)	284 (38)	399 (43)	23 (32)
Median age (Q1; Q3)b	67.0 (56.0; 75.0)	61.0 (44.0; 73.0)	59.0 (47.0; 77.8)
Patients with pneumonia (%)	298 (40)	936 (100)	9 (12)
Median APACHE II scorec (Q1; Q3)b	20.0 (16.0; 24.0)	–	18.0 (14.0; 23.8)
Median SOFA scored (Q1;Q3)b	6.79 (4.94; 9.50)e	3.0 (2.0; 4.0)f	7.5 (5.0; 10.0)g
With microbiology (%)	603 (81)	612 (65)	74 (100)
Any pathogen identified (%)	534 (89)h	208 (34)h	69 (93)h
Gram-positive or gram-negative bacterial infection (%)	496 (82)h	177 (30)h	60 (81)h
Gram-positive infection only (%)	358 (59)h	–	8 (11)h
Gram-negative infection only (%)	324 (54)h	–	52 (70)h
Fungal infection (%)	172 (29)h	62 (10)h	2 (3)h

In PROGRESS the outcome is defined as death within 28 days or qualified intensive-care requiring ventilation, treatment with catecholamines, oxygenation or dialysis.

First and third quartile.

Acute Physiology and Chronic Health Evaluation II score.

Sequential Organ Failure Assessment score, in the discovery GWAS.

14 day mean SOFA data from 714 of 740 patients.

Worst SOFA score within 5 days.

SOFA score at baseline.

Percentages relative to the 603 or 612 patients with microbiology.

Procedures

For the GWAS data (HumanOmniExpressExome arrays) we applied stringent measures of quality control (QC) to remove unreliably genotyped patients or SNPs, population outliers as determined by performing a principal component analysis of the genome-wide data, and samples for which there were sex discrepancies (details see Appendix). The number of autosomal SNPs remaining for imputation were 644,699 which were subsequently imputed using IMPUTE2 (version 2.3.0) and with 1000 Genomes Project data (phase 1, version 3) as a reference panel. After additional QC of the imputed data, 7,993,459 SNPs were finally available for the genome-wide analysis (details see Appendix).

(1) Genotyping of the patients in GenOSept/GAinS was performed on Affymetrix 5.0 SNP arrays and Illumina Human OmniExpressBeadChip SNP arrays. VASST and PROWESS were both genotyped by Illumina Human 1 M-Duo BeadChip SNP array. All datasets were also imputed separately with IMPUTE2 and with 1000 Genomes Project data as a reference panel. For details see Rautanen et al. (2015). (2) Genotyping of the patients in PROGRESS was performed using the Affymetrix Axiom-CAP2 microarray. The CAP2 array is a genome-wide custom microarray. It contains Axiom-CEU content as backbone but is enriched with candidate SNPs and regions. Genotype calling was performed with Affymetrix power tools (version 1.15.1). Sample filters comprise dish-qc < 0.82, call-rate < 97%, implausible dish-qc vs. call-rate, implausible relatedness, sex-mismatches and ethnic outliers identified by the 6SD outlier criterion of SMARTPCA. SNP filtering comprise the cluster plot quality metrices proposed by Affymetrix (HetSO, HomRO, FLD), call-rate ≤ 97%, p-value of exact test for Hardy-Weinberg equilibrium ≤ 10− 6, p-value of plate association ≤ 10− 7, exclusion of monomorphic SNPs and non-autosomal SNPs. A total of 589,205 SNPs fulfilled these criteria and were used for imputation with the reference panel 1000 Genomes Project data (phase 1, version 3). SHAPEIT v2.r790 was used for pre-phasing and IMPUTE2 v2.3.1 was used for final imputation.

For the exome sequencing study, 2–3 μg DNA per sample was fragmented on a Covaris M220 focused ultra-sonicator. Exomes were enriched using Agilent SureSelect XT Human All Exon V5 + UTRs kits targeting 74,856,280 bp in the coding sequence and untranslated regions of 20,791 genes. The mean depth of sequence coverage was 91-fold (range: 52- to 159-fold). Relative to the human reference genome (GRCh37/hg19) we called 313,279 single nucleotide variants using the Genome Analysis Toolkit (GATK 2.5) (DePristo et al., 2011) and we again refer to Taudien et al. 2016 for details.

Statistical Analyses

We analyzed all autosomal GWAS variants for the dichotomous outcome 28 day mortality by logistic regression (log-additive genetic model as implemented in SNPTEST, version 2.5) with age (linear), sex, and the first three principal components as covariates (model 1). Age and sex are known to be strong determinants of mortality in patients with sepsis (Martin et al., 2003) and principal components were used to avoid confounding due to population structure. In addition, we added APACHE II scores (linear) as covariates as a summary measure of baseline morbidity (model 2). Finally, we performed sensitivity analyses by also including indicator variables for the treatment arm of the VISEP/MAXSEP trials for selected variants (p ≤ 10− 5 in the primary GWAS analysis). Association signals at SNP variants were summarized as GWAS loci if more than one SNP signal had a p-value ≤ 10− 5 within a region of ± 500 kb around the lead SNP. Details on the analysis of the X-chromosome and the corresponding results, which were not the focus of this report, are provided in the Appendix.

In the Appendix, we also provide details on the comparison-wise statistical power of our discovery GWAS with n = 740 patients to detect an association with the dichotomous 28 day mortality outcome after treated sepsis. Overall, our GWAS had a power > 80% to detect strong genetic effects on the 28 day mortality outcome.

(1) In short, the statistical analyses in the GWA studies (GenOSept/GAinS; VASST; PROWESS) were similar to those of the discovery GWAS; for details we again refer to Rautanen et al. (2015). (2) In PROGRESS we considered two outcome measures: first a combined binary endpoint of death within 28 days or necessity of intensive care (defined by a first-time requirement of ventilation, oxygenation, dialysis or catecholamines), and second, the worst SOFA score within five days after enrolment. Association analysis was performed using SNPTEST version 2.5 assuming a log-additive genetic model. We adjusted for age (linear), sex, and the first three principal components. A total of 936 individuals had complete genetic, phenotypic and covariate data.

We analyzed the exome sequencing data using the adjusted SKAT-O method as implemented in the package “SKAT” (Lee et al., 2012) in R (version 3.1.1). SKAT-O is a method to assess the cumulative effects of all variants in a genomic region (in our case gene coding regions ± 10 kb according to UCSC Genome Browser as reference (GRCh37/hg19)). Depending on the expected power for each region, SKAT-O runs as a “burden” or “nonburden” test. We searched for genomic associations with the dichotomous outcome 28 day mortality adjusting for age (linear), sex, and the first two principal components for population structure as covariates. In addition, we performed also sensitivity analyses by including the first five principal components and center as covariates. Note that the principal components were calculated within the exome study. We limited the SKAT-O result presentation to our 21 notable genes in proximity to the best 14 discovery GWAS SNPs and to the 14 reported genes in proximity to the best 35 GWAS SNPs from Rautanen et al. (2015). In particular, given that SKAT-O results summarize the evidence at the gene level, we decided to report all gene-level statistics up- and downstream of our discovery GWAS SNPs for cases in which the discovery GWAS SNP signal was not located within the coding region.

Bioinformatic Annotations and In Silico Analyses

We created regional plots of the best GWAS loci using LocusZoom (Pruim et al., 2010). As part of the functional assessment of the detected variants we assessed their potential deleteriousness as described in Kircher et al. (2014). These authors developed the C-score which is based on Combined Annotation Dependent Depletion. It estimates the deleteriousness of a SNP or indel across a number of bioinformatic tools (for details we refer Kircher et al., 2014) with larger values indicating higher pathogenicity potential. Next, we searched the NHGRI-EBI GWAS catalog (https://www.ebi.ac.uk/gwas/assessed at 07/25/2016; Welter et al., 2014) for all notable genes whereas pathways (Table 4) were checked using http://www.genecards.org/assessed at (07/25/2016). For loci with support from more than one data source, we applied a previously described approach (Geisel et al., 2016) using GTEx (release version 6) and ENCODE and the UCSC Genome Browser (GTEx Consortium, 2015, Rosenbloom et al., 2013; see Appendix) in order to further annotated the loci.

Table 4

Validation of the notable genes from the discovery GWAS in the independent exome sequencing study for the outcome 28 day mortality among patients with treated sepsis. Note that for 5 of the 21 notable genes no SKAT-O p-values could be calculated due to sparse data.

Chromosome	Analyzed regiona	Notable genesb	Pathway(s)c	# low frequency or common variants in survivors/non-survivorsd	# rare variants in survivors/non-survivorse	p-Valuef
1	68,884,506–68,925,642	RPE65	(Visual) signal transduction (by GPCRg), retinol metabolism	8/7	3/1	0.12
1	68,929,834–68,972,904	DEPDC1	No pathway known, GO-term: GO: 0007165; signal transduction	13/5	5/0	0.13
3	11,284,384–11,314,939	HRH1	7 super pathways, among them GPCRg activity and histamine binding	3/1	1/0	0.23
3	37,483,812–37,871,281	ITGA9	25 super pathways, among them GPCR g and signal transduction by L1	34/21	11/5	0.35
3	37,785,179–37,913,271	ITGA9-AS1	–	4/1	2/0	0.21
3	187,920,720–188,618,460	LPP	Stabilization and expansion of the E-cadherin adherens junction	14/13	16/9	0.48
3	194,008,988–194,040,593	LINC00887	–			0.84
4	833,064–936,174	GAK	Vesicle budding, membrane trafficking	37/30	28/14	0.14
9	79,782,360–80,007,921	VPS13A	–	55/24	22/3	0.04
12	22,768,075–22,807,349	ETNK1	3 super pathways, phospholipid metabolism	1/1	1/1	0.27
12	23,675,230–23,747,546	SOX5	4 super pathways, among them ERK signaling	6/5	3/0	0.33
13	27,319,338–27,344,922	GPR12	Peptide ligand binding receptors, GPCRsg, class A rhodopsin-like	3/3	0/0	0.32
13	27,630,286–27,756,033	USP12	Ubiquitin-proteasome dependent proteolysis	3/2	4/1	0.09
13	70,264,724–70,692,625	KLHL1	–	19/12	6/1	1.00
16	84,843,586–84,953,116	CRISPLD2	–	32/35	12/9	0.003
17	14,194,505–14,259,492	HS3ST3B1	Heparan sulfate biosynthesis/metabolism	3/1	1/0	0.45

According to GRCh37 (hg19) ± 10 kb.

See Table 1.

http://www.genecards.org/.

In the analyzed region - MAF > 0.005 as reported in at least one of three databases (ExAC non-Finnish European group or ESP Americans of European ancestry or dbSNP – see Taudien et al. 2016).

In the analyzed region - MAF < 0.005 or not reported in the three databases.

Exact two-sided p-values for the SKAT-O analyses.

G-protein coupled receptor.

Results

We identified 243 autosomal variants with suggestive evidence for an association with 28 day mortality in the discovery GWAS of 740 European patients with treated sepsis (p < 10− 5, Table 1, Fig. 1). These variants clustered in 14 GWAS loci (Table 2) and were further validated in two GWA studies and in an independent exome sequencing study (Taudien et al., 2016). In addition, we identified one locus in males on the X-chromosome (see Supplementary Fig. 15 and Supplementary Table 3).

Fig. 1

Discovery GWAS: a) Manhattan plot for the analysis of 28 day mortality in patients with treated sepsis (additive genetic model). The blue line represents the significance level for a suggestive association signal (i.e. p ≤ 1 × 10− 5) and the red line is the level for a genome-wide association signal (i.e. p ≤ 5 × 10− 8). b) Quantile–quantile (Q–Q) plot of all the GWAS p-values (λ = 1.0079) from the analysis of 28 day mortality in patients with treated sepsis (additive genetic model).

Table 2

Discovery GWAS top association signals (loci with SNP association signals that met a p-value < 10− 5) for 28 day mortality among patients with treated sepsis.

SNP	Chromo-some	Physical positiona	Variant type	Effect allele/other allele	Effect allele frequencyb	Model 1 (adjusted for sex, age, PCc)		Model 2 (adjusted for sex, age, PCc & APACHE II)		Notable genes
						Estimated odds ratio	p-Value	Estimated odds ratio	p-Value
rs382422	1	68,916,123	Intergenic	C/G	0.22	2.1	3.21 × 10− 6	2.6	8.98 × 10− 7	RPE65DEPDC1
rs58764888	3	11,217,691	Intronic	A/T	0.02	13.3	6.70 × 10− 7	15.0	3.60 × 10− 7	HRH1
rs72862231	3	37,853,059	Intronic; NCTd	A/T	0.05	4.4	1.73 × 10− 6	5.0	5.12 × 10− 7	ITGA9ITGA9-AS1
rs150062338	3	188,004,948	Intronic; regulatory region	T/C	0.01	38.6	2.32 × 10− 7	26.1	2.03 × 10− 6	LPP
rs10933728	3	194,027,568	Intronic; NCTd	G/A	0.03	7.0	5.62 × 10− 6	7.8	3.37 × 10− 6	LINC00887
rs115550031	4	856,102	Intronic; NCTd	A/G	0.02	13.8	2.45 × 10− 6	17.6	7.41 × 10− 7	GAK
rs62369989	5	117,409,248	Intronic; NCTd	G/T	0.26	2.1	7.98 × 10− 6	2.0	4.07 × 10− 5	LOC102467224
rs115036193	6	33,000,554	Intronic	T/C	0.01	16.2	2.21 × 10− 6	11.3	2.57 × 10− 5	HLA-DOAHLA-DPA1
rs117983287	9	80,020,874	Missense	A/C	0.01	18.2	8.16 × 10− 8	12.1	2.18 × 10− 6	VPS13A
rs150811371	12	23,661,042	Intergenic	A/G	0.08	3.4	2.93 × 10− 6	4.0	4.46 × 10− 7	ETNK1SOX5
rs945177	13	27,621,985	Intergenic	A/G	0.02	14.7	1.31 × 10− 6	12.4	5.46 × 10− 6	GPR12USP12
rs9529561	13	69,899,506	Intergenic	G/A	0.08	3.9	3.34 × 10− 7	3.6	1.68 × 10− 8	LINC00550KLHL1
rs2641697	16	84,885,777	Intronic; NCTd	G/C	0.36	2.0	5.99 × 10− 6	2.0	2.27 × 10− 5	CRISPLD2
rs7211184	17	14,257,083	Intergenic; regulatory region	C/G	0.72	2.0	9.43 × 10− 6	2.0	5.04 × 10− 6	HS3ST3B1CDRT7

According to GRCh37 (hg19).

Estimated effect allele frequency in all GWAS patients.

Principle components to address potential population stratification effects.

Non-coding transcript variant.

From the 14 lead SNPs that are markers of the 14 discovery GWAS loci (see Fig. 2 and Appendix for regional plots), a total of 9 are located within the coding region of genes and 8 have a minor allele frequency (MAF) ≤ 5% in our patient sample. Among them is the top associated missense variant (rs117983287; p = 8.16 × 10− 8) located at chromosome 9q21.2 in the VPS13A (vacuolar protein sorting 13 homolog A (S. cerevisiae)) gene (Fig. 2a, Supplementary Fig. 18). With a scaled C-score of 22.3 this variant is predicted to be among the 1% most deleterious substitutions that one can have in the human genome (Kircher et al., 2014). Our sensitivity analyses showed that the estimated odds ratios were robust with respect to treatment groups and were sensitive to the adjustment for the Acute Physiology and Chronic Health Evaluation (APACHE) II score (largest change of the estimated odds ratio (33% decrease) of all 14 discovery GWAS signals). Furthermore, we observed an impact of the risk variant on the 14 day average SOFA scores particularly for the cardiovascular and the renal sub-scores (Supplementary Table 1c). While no information for this variant was available in the meta-analysis data set by Rautanen et al. (2015) (see Table 3), the variant was supported by the independent PROGRESS GWAS data (uncorrected p = 0.027 for the association to the worst SOFA score within five days after enrolment; Table 3). Interestingly, the results for VPS13A in the independent exome sequencing study also indicated evidence for an association to 28 day mortality (pSKAT-O = 0.04; Table 4). GTEx and ENCODE annotations are provided in Supplementary Fig. 19.

Fig. 2

Regional association plot for a) the chromosome 9q21.2 locus (centered around the lead SNP rs117983287) and b) the chromosome 16 locus (centered around the lead SNP rs2641697) in the analysis of 28 day mortality in patients with treated severe sepsis/septic shock (additive genetic model). Colors indicate the correlation (r2 in 1000 Genomes data for Utah residents with northern or western European ancestry (CEU); phase 1, version 3) with the alleles of rs117983287/rs2641697.

Table 3

Validation of the autosomal SNP markers from the discovery GWAS in two independent GWA studies.

SNP	Chromo-some	Physical positiona	Effect allele/other allele	Meta-analysis of three discovery GWAS cohorts of Rautanen et al. (2015)				PROGRESS GWAS				Meta-analysis p-valueb
				Patients with sepsis caused by pneumonia or abdominal infections		Patients with sepsis caused by pneumonia		Death within 28 days or necessity of intensive care		Worst SOFA within five days after enrolment
				Odds ratio for the effect allele	p-Value	Odds ratio for the effect allele	p-Value	Odds ratio for the effect allele	p-Value	β for the effect allele	p-Value
rs382422	1	68,916,123	C/G	1.01	0.949	0.99	0.924	1.09	0.673	0.13	0.248	0.002
rs58764888	3	11,217,691	A/T	0.71	0.162	1.07	0.850	0.64	0.525	-0.44	0.224	0.090
rs72862231	3	37,853,059	A/T	1.08	0.646	1.25	0.252	0.57	0.230	0.18	0.445	0.003
rs150062338	3	188,004,948	T/C	–	–	–	–	–	–	-0.40	0.368	–
rs10933728	3	194,027,568	G/A	0.69	0.045	0.66	0.079	0.51	0.307	0.24	0.445	0.292
rs115550031	4	856,102	A/G	–	–	–	–	0.77	0.765	0.04	0.935	2.1 × 10− 5
rs62369989	5	117,409,248	G/T	0.91	0.252	0.87	0.186	1.30	0.219	0.09	0.464	0.020
rs115036193	6	33,000,554	T/C	–	–	–	–	1.15	0.821	0.22	0.554	3.1 × 10− 6
rs117983287	9	80,020,874	A/C	–	–	–	–	1.47	0.569	0.95	0.027	1.1 × 10− 7
rs150811371	12	23,661,042	A/G	1.03	0.803	1.03	0.871	0.83	0.547	-0.13	0.464	0.002
rs945177	13	27,621,985	A/G	–	–	–	–	1.09	0.873	0.11	0.761	2.9 × 10− 6
rs9529561	13	69,899,506	G/A	–	–	–	–	1.87	0.035	-0.12	0.508	2.4 × 10− 10
rs2641697	16	84,885,777	G/C	0.98	0.768	1.15	0.132	0.91	0.608	0.03	0.790	0.012
rs7211184	17	14,257,083	C/G	–	–	–	–	0.82	0.319	-0.04	0.689	3.8 × 10− 4

According to GRCh37 (hg19).

Meta-analysis based on one-sided p-values of the discovery GWAS (model 1) and both validation GWA studies (Rautanen et al. (2015): sepsis caused by pneumonia or abdominal infections; PROGRESS: death within 28 days or necessity of intensive care) using the weighted inverse normal (z-score) method (weighting by n; e.g. 740, 2534, 936).

Among the GWAS loci with more frequent lead SNPs (estimated MAF > 5%) and corresponding validation GWAS or exome sequencing association signals, we observed p = 3.34 × 10− 7 for an intergenic SNP (rs9529561) on chromosome 13q21.33 which was also supported by the PROGRESS data (uncorrected p = 0.035 for the association to death within 28 days or necessity of intensive care; p = 2.4 × 10− 10 in the meta-analysis). While again no data was available from the meta-analysis by Rautanen et al. (2015), the exome data for the KLHL1 (kelch like family member 1) gene results in 375 kb distance did not provide evidence for an association (pSKAT-O = 1.00; Table 4). For an intronic SNP (rs2641697) of the CRISPLD2 (cysteine rich secretory protein LCCL domain containing 2) gene on chromosome 16q24.1, we observed p = 5.99 × 10− 6 in the discovery GWAS and pSKAT-O = 0.003 for 28 day mortality in the exome data. However, this time none of the independent validation GWA studies supported this observation. GTEx and ENCODE annotations for the two loci are provided in Supplementary Figs. 20 and 21.

Finally, we compared the top association findings (p ≤ 10− 5) from Rautanen et al. (2015) with the results from our discovery GWAS and exome sequencing study (Supplementary Table 2). Within the reported 12 loci with variants associated with 28 day mortality among patients with sepsis caused by pneumonia or abdominal infections, we observed p = 0.01 for rs2096460 located in URB1 (URB1 ribosome biogenesis 1 homolog (S. cerevisiae)) on chromosome 21q22.11 and p = 0.04 for an intergenic variant (rs74438932) on chromosome 13 located 88 kb downstream of GPR12 (G protein-coupled receptor 12) gene on chromosome 13q12.13. However, the effect alleles in both GWAS were directionally inconsistent and none of the available exome sequencing association signals achieved a pSKAT-O ≤ 0.05. For details we refer to Supplementary Table 2.

Discussion

We report results of a GWAS in patients with treated sepsis which focused on common genetic variants associated with 28 day mortality. We validated our best findings using three independent data sets including two GWA and a whole-exome sequencing study. We applied the new Sepsis-3 definition (Singer et al., 2016, Seymour et al., 2016) in the discovery GWAS and the sequencing study which requires the presence of organ dysfunction for a diagnosis of sepsis. Furthermore, we followed-up the best discovery loci from the most recent and largest sepsis GWAS which applied a similar study design (Rautanen et al., 2015).

The SNP with the strongest GWAS signal was a missense and potentially deleterious variant located on chromosome 9q21.2 within VPS13A. This result was supported by the validation GWAS and the exome sequencing data. Recent experiments (Muñoz-Braceras et al., 2015) indicated an important regulatory role of VPS13A for autophagic degradation. Autophagy is a key component of our immune system and has also been associated to several human diseases (Cuervo and Macian, 2014, Schneider and Cuervo, 2014). However, this signal is located within a gene-rich region (9q21) that has been associated to many complex diseases like mental disorders, type 2 diabetes mellitus, some cancers and cardiovascular disease (An et al., 2012, Shimo et al., 2011). Thus, VPS13A might not be the only candidate. Notably, we also observed a similarly strong association GWAS signal for variants in GNA14 (guanine nucleotide binding protein (G protein), alpha 14), a gene which is located ~ 200 kb distal to VPS13A. GNA14 is member of the “G alpha Q signaling events”-pathway which is highlighted in the accompanying report by Taudien et al. 2016. In their report, rare variants in genes of the pathway, including GNA14, were found to have a putatively protective effect leading to a more favorable sepsis course. Identification of GNA14 in both studies is a strong argument that the G alpha Q signaling pathway might play an important role in sepsis. In contrast to our findings, Rautanen et al. (2015) did not list this 9q21 locus around VPS13A and GNA14 among their top findings. Potential explanations, apart from a false positive finding in our data sets, could be their exclusion of variants with minor allele frequencies < 2% and/or inclusion of less severely affected patients by applying the 2001 sepsis definition (Levy et al., 2003) or the general need for much larger samples size as also underlined by our power considerations.

Besides our strongest GWAS signal, we also found association signals which have previously been associated with complex diseases including sepsis phenotypes. Genomic variation in CRISPLD2 was reported to be associated with the presence of cleft lips (Mijiti et al., 2015), and recent work showed a decreased expression of CRISPLD2 in septic shock and an association with procalcitonin – one of the best validated biomarkers in sepsis research (Wang et al., 2013). Furthermore, we observed some converging support for an intergenic region on chromosome 13q21.33 which was previously reported to be associated with risk of chronic kidney disease (Köttgen et al., 2010) in the GWAS catalog. However, these authors described an intronic SNP within the DACH1 (Dachshund homolog 1) gene and the SNP alleles were in linkage equilibrium. Summarizing these considerations, the region on chromosome 9q21.2 including VPS13A, CRISPLD2 and to a lower extent the chromosome 13q21.33 locus are regions with a biologically plausible relationship to sepsis which are supported by both our data sets and external evidence.

When focusing on the top association findings from Rautanen et al. (2015), we already reported that we could not support the clinical implications suggested for the FER gene (Schöneweck et al., 2015). Moreover, our data did not strongly support other candidate genes either. However, Rautanen et al. (2015) focused on sepsis due to pneumonia which might in part explain these discrepancies.

Our study has several strengths and limitations. Firstly, our sample size was rather small so that the discovery GWAS was underpowered to detect smaller genetic effects while properly controlling the (family-wise) type I error rate. To protect against false positives, we assessed two independent GWA validation studies and a sample of sepsis patients by exome sequencing. Here it should also be noted that none of the markers formally replicated at a Bonferroni-corrected 0.05/14 ≈ 0.003 significance level. Secondly, the discovery GWAS patients included were selected patients from two RCTs (Brunkhorst et al., 2008, Brunkhorst et al., 2012) and treatment - though not affecting the outcome differently - might have had a differential genotype-dependent effect on 28 day mortality. To address this limitation we conducted sensitivity analyses adjusting for treatment group and obtained fairly similar i.e. robust findings. Thirdly, a GWAS has a focus on common variation scattered across the genome while exome sequencing detects both rare and common coding variants in the coding region. For the given reasons one might argue that mixing both strategies is a bad idea; others might argue that functionally relevant variation with stronger effects might still be most likely detectable in the exome. Fourthly, even after applying the new Sepsis-3 definition, sepsis remains a highly complex phenotype. Following Rautanen et al. (2015) more promising results might be generated by focusing on more homogenous subgroups such as sepsis patients with pneumonia. We agree that better defined subgroups (e.g. identical pathogens) might help overcome some of the challenges faced in sepsis research (Cohen et al., 2015). Here we applied the new Sepsis-3 definition to both the discovery GWAS and the sequencing study and we wanted to avoid multiplicity issues that arise if subgroups are defined post-hoc (Sun et al., 2014, Burke et al., 2015). However, the validation GWA studies addressed a slightly different phenotype spectrum (e.g. patients with CAP) and not exactly the same outcomes given that e.g. death is a rare event in CAP patients. Fifthly, 28 day mortality has been criticized as an relatively arbitrary and unspecific endpoint (Cohen et al., 2015) and clearly genetic studies of the disease course after sepsis would be strengthened by intermediated omic or endophenotype data (between genomic variation and clinical outcomes). Unfortunately, such additional omic data layers were not available to us in the same samples. Recent reviews on sepsis in the field of proteomics (Camprubí-Rimblas, 2015), metabolomics (Fanos et al., 2014) and transcriptomics (Almansa et al., 2015), however, did not imply any of the above mentioned genes. Given that omic data integration in sepsis research is still in its infancy, we think that this is an important facet of future sepsis-related genomic research as recently demonstrated by Davenport et al. (2016).

In summary, we performed a discovery GWAS in patients with treated sepsis which focused on common genetic variants associated with 28 day mortality and validated our best gene loci in independent GWA studies and an exome sequencing study. GWAS and exome data supported the VPS13A gene locus on chromosome 9q21.2. Furthermore, we identified one region on chromosome 13q21.33 and one candidate gene CRISPLD2 which should be a focus of future omic research activities in order to elucidate their biological influences. Future genome-wide studies in the field of sepsis will only be successful if more homogenous phenotype definitions in much larger samples are applied.

Funding Sources

The project was supported by the Paul-Martini-Sepsis Research Group, funded by the Thuringian Ministry of Education, Science and Culture (ProExcellence; grant PE 108-2); the public funded Thuringian Foundation for Technology, Innovation and Research (STIFT) and the German Sepsis Society (GSS); the Jena Center of Sepsis Control and Care (CSCC), funded by the German Ministry of Education and Research (BMBF; 01 EO 1002, 01 EO 1502). The VISEP and MAXSEP trials from the SepNet Study Group had been supported by a BMBF grant (01 KI 0106) and by unrestricted grants from B. Braun, HemoCue, Novo Nordisk, Astra Zeneca GmbH, Wedel, Germany and Bayer HealthCare, Leverkusen, Germany. The exome sequencing study was funded in part by the Hellenic Institute for the Study of Sepsis. The GenOSept study was supported by the European Union and benefits from the 6th framework programme of RTD funding. The PROGRESS study is supported by the German Federal Ministry of Education and Research, grant numbers 01KI07110 (Giessen), 01KI07111 (Jena), 01KI07113 (Leipzig), 01KI07114 (Berlin), 01KI1010I (Leipzig), and 01KI1010D (Greifswald).

None of the funders had a role in study design, data collection, data analysis, interpretation, writing of the report.

Conflicts of Interest

EGB declares that he received honoraria for providing scientific advice to AbbVie, Chicago IL, USA; Astella, Athens, Greece; Biotest AG, Dreieich, Germany; and ThermoFisher Scientific GmbH, Hennigsdorf, Germany. EGB also received unrestricted educational funding (paid to the University of Athens) by Biotest AG, Dreieich, Germany; Sanofi SA, Athens, Greece; ThermoFisher Scientific GmbH, Hennigsdorf, Germany and by the Seventh Framework European Program HemoSpec. All other author declare no conflicts of interest.

Author Contributions

Study concept and design: All authors.

Acquisition, analysis, or interpretation of data: All authors.

Drafting of the manuscript: Scherag, Schöneweck, Scholz, Brunkhorst.

Critical revision of the manuscript for important intellectual content: All authors.

Statistical Analysis: Scherag, Schöneweck, Kesselmeier, Trips, Scholz.

Obtained funding: Scherag, Taudien, Platzer, Kurzai, Slevogt, Brunkhorst.

Administrative, technical, or material support: Felder, Sponholz, Rautanen, Hill, Hinds, Hossain, Suttorp, Giamarellos-Bourboulis, Armaganidis.

Study supervision: Scherag.

Scherag and Schöneweck are joint first authors and Scholz and Brunkhorst are joint last authors.