This study was approved by the Research Ethics Committee of the participating center (Hospital Universitario Nuestra Señora de Candelaria: CHUNSC_2018-16). We obtained written informed consent from all participants or an appropriate proxy. We used a case-control approach with DNA samples from 763 subjects. Controls included 416 individuals from the Cardiovascular, Diabetes and Cancer cohort study (24). All individuals declared at least two generations of ancestors born in the Canary Islands, as has been described elsewhere (21). A total of 347 patients with sepsis from the GEN-SEP study admitted into a network of Spanish ICUs, with at least two generations of ancestors born in the Canary Islands, were used as cases. Sepsis was defined according to the Third International Consensus Definitions for Sepsis (1) (see Supplementary Methods in Supplementary Material for details).

The Axiom® Genome-Wide Human CEU 1 Array (Thermo Fisher Scientific, Waltham, MA) was used for genotyping 587,352 variants in DNA from donors with the support of the National Genotyping Center (CeGen), Universidad de Santiago de Compostela Node (Spain). Quality control procedures of data were performed with R 3.2.2 and PLINK v1.07 (25). Samples with a genotype call rate <95% or family relationships (PIHAT>0.2) with others were removed, resulting in a total of 753 individuals for further analyses, (343 were patients with sepsis and 410 population controls). Moreover, SNP filtering based on genotyping rate <95%, minor allele frequency (MAF) <0.01, and large deviations from Hardy-Weinberg equilibrium expectations (p < 1 × 10⁻⁶) left a total of 494,390 variants.

To obtain the ancestry estimates, and to maximize the intersection of autosomal SNPs in subsequent analyses from the datasets of cases, controls, and reference populations, we followed the methods described elsewhere (21). We extracted EUR and SSA datasets from the 1000 Genomes Project (1KGP) Phase 3 data (26). The NAF representation was gathered from 125 samples with origins in North and South Morocco, Western Sahara, Algeria, Tunisia, Egypt, and Libya that were previously genotyped with the Genome-Wide Human SNP Array 6.0 (Affymetrix, Santa Clara, CA) (27). Genotyping quality controls were performed using PLINK v1.07. EUR, NAF, and SSA individuals with genotype call rates <95% were removed from the analysis, as well as SNPs with >5% missing rate or deviating from the Hardy-Weinberg equilibrium expectations (p < 1 × 10⁻⁶) in each population. The intersection and post-filtering resulted in 113,414 SNPs genotyped for downstream analyses in cases, controls, and reference populations (see Supplementary Methods in Supplementary Material for details).

ADMIXTURE v1.3 (28) was used to estimate proxies for global ancestry proportions and verify that cases and controls were similar in these terms. Briefly, cross-validation error was lowest for k = 4, allowing to differentiate the SSA, NAF and EUR components clearly. The latter was detected as two separate ancestries (29) that were considered in aggregate in the analyses for simplicity. ADMIXTURE results were represented with CLUMPAK (30) (see Supplementary Methods in Supplementary Material for details).

Admixture mapping analyses were based on local ancestry block estimates across autosomes, which were inferred using Efficient Local Ancestry Inference (ELAI v1.0) (31), and assuming three admixing populations (EUR, NAF, and SSA) as has been described elsewhere (32). Association testing was performed with EPACTS v3.2.6 software (33) between local ancestry scores at each genomic position and sepsis susceptibility for each ancestry separately. Resulting p-values were corrected using a genomic control strategy based on λ calculation (34). The significance threshold was declared as p = 1.82 × 10⁻⁴ after Bonferroni correction for an average of 276 ancestry blocks, as has been inferred for this population (21). Regional plots of association results were represented with LocusZoom (35). The Power Analysis in Multiancestry Admixture Mapping (PAMAM) (36) was used to calculate the statistical power of the study. Based on PAMAM, assuming a design of 300 cases and 400 controls, a significance threshold of 1.82 × 10⁻⁴, and the known EUR admixture, this study achieved >80% power to detect an Odds Ratio (OR) <0.55 (Supplementary Figure 1).

We performed a fine-mapping study of the significant admixture mapping region. For this, SNP imputation of the whole chromosome containing the region was conducted with the Michigan Imputation Server tool (37) using the Haplotype Reference Consortium (HRC) version r1.1 2016 as the reference panel (38) and estimating haplotypes with Shape-IT v2.r790 (39). EPACTS v3.2.6 (33) was used to perform SNP association testing with allele dosage data of those variants with MAF >0.01 and imputation quality (Rsq) >0.3 lying in the region. A Bonferroni correction was performed to prioritize the most significant variant within the 1.2 Mb significant admixture mapping region. Based on a total of 449 independent tests examined in the region, significance was established at p < 1.57 × 10⁻⁴ (suggestive significance at p < 3.13 × 10⁻³). Conditional logistic regressions were assessed using R 3.6.0 to reveal the independent SNP associations near the significant admixture mapping region.

We used iSAFE v1.0.4 (40) to provide evidence of a selective sweep embedded in the significant admixture mapping region and to pinpoint the most likely favored variant. iSAFE exploits the evolutionary contributions hidden in the flanking regions surrounding the region under selection to provide a ranking of variants (iSAFE-score) based on their contribution to the overall signal of selection.

For the analysis, we used phased data from 59 unrelated Iberians (IBS) from 1KGP and a random selection of 10 YRI subjects (out of the 108 available) to represent a non-target or outgroup population in GRCh37/hg19 coordinates. Data from chromosome 8p23.1 from 59 unrelated NAF samples were also used for comparisons. Since 1KGP does not include NAF datasets, we accessed data from whole-genome sequence from donors originating across NAF populations (Algeria [n = 3], Berber [n = 6], Egypt [n = 2], Libya [n = 2], Morocco [n = 8], Western Sahara [n = 4], Tunisia [n = 32], North Africa [undisclosed country, n = 2]). Seventeen of these samples were sequenced to a mean coverage >25X with Illumina HiSeq 2000 using 101 bases with paired end reads (41). Another 42 samples were sequenced using 150 paired-end reads either with Illumina HiSeq 4000 (n = 24, mean coverage >36X) or Illumina NovaSeq 6000 (n = 18, mean coverage >20X). Variant calling of this data was performed following best practices with an in-house pipeline based on the Burrows-Wheeler Alignment Tool (BWA-MEM v.0.7.12-r1039) and GATK v4 using the GRCh37/hg19 as a reference. Resulting calls were filtered to keep biallelic SNPs flagged as PASS. This callset was phased into haplotypes without providing any reference data using Eagle v2.4.1 (42).

iSAFE was executed in the region of interest enabling the IgnoreGaps flag and the default MaxFreq value (0.95). Ancestral fasta sequences for Homo sapiens (GRCh37) were downloaded from ENSEMBL release 75 (http://ftp.ensembl.org/pub/release-75/fasta/ancestral_alleles/). Given the exceptional performance of iSAFE in prioritizing the most likely favored variant in 94% of the times among the SNPs with the highest iSAFE scores (40), we performed the functional annotation of the 20 top ranked SNPs with the best scores using the Variant-to-Gene (V2G) pipeline aggregation from Open Target Genetics (43).

To assess the functional role of the most significant variant and its best proxies (r² >0.8), we accessed empirical data from different integrated online software tools, including LDLink v.4.1.0 webtool (44), The Open Targets Post-GWAS webtool (43), GTEx v.7 (45), HaploReg v4.1 (46), RegulomeDB v.2.0 (47), SNPDelScore (48), CHiCP (49), and the 3D Genome Browser (50). We also accessed three publicly available gene expression data, GSE32707 (51), GSE57065 (52, 53), and GSE28750 (54), from the Gene Expression Omnibus (GEO) data repository to assess differential expression of the genes residing in the significant admixture mapping region (plus 1 Mb on each side). The gene expression of all genes located within the region were compared using a Student's t-test. Finally, we performed a meta-analysis through a Fisher test of these GEO data sets using ImaGEO (55) (see Supplementary Methods in Supplementary Material for details).

Basic clinical and demographical data of cases and controls are provided in Table 1. Further details are included in the supplementary material (Supplementary Results; Supplementary Table 1; Supplementary Figure 2).

The admixture mapping study was performed using EUR, NAF, and SSA local ancestry estimates for 113,414 variant positions across the genomes of all cases and controls. No inflation effects were identified in the association results for any ancestry (λEUR = 0.96, λNAF = 0.96, and λSSA = 0.99). Association testing of the three ancestry scores with sepsis susceptibility (Figure 1) revealed 114 consecutive positions at 8p23.1 (chr8: 8,155,475–9,318,404; GRCh37/hg19 coordinates) associated with sepsis protection. The strongest significance was obtained for rs17149618 (p = 1.37 × 10⁻⁴; OR [95%CI] = 0.51 [0.40–0.66]) and it was related to a single local EUR ancestry peak of significance (Figure 1; Supplementary Figure 3; Supplementary Table 2). For rs17149618, the NAF ancestry associated with sepsis risk (p = 1.93 × 10⁻⁴; OR [95%CI] = 2.01 [1.57–2.58]). The SSA ancestry blocks were unrelated with this result and lacked significant associations with sepsis in the study population (Figure 1).

The region of 8p23.1 associated with sepsis in this study is broad (~1.2 Mb) and harbors several genes. In particular, most of the 114 significant positions were intergenic, although 27 of them corresponded to positions of the genes encoding the Malignant Fibrous Histiocytoma Amplified Sequence 1 (MFHAS1), Exoribonuclease 1 (ERI1), and PEAK1 Related, Kinase-Activating Pseudokinase 1 (PRAG1) (Supplementary Table 2). Strikingly, we noted a region telomeric to the significant admixture mapping region of 8p23.1 showing a gap of mapped genetic variation (Figure 2) that we interpreted as a consequence of repetitive elements associated with the cluster of defensin genes mapping in the region (DEFB103B, DEFB103A, DEFB109P1B, DEFB4A, DEFB4B, and DEFB130). Therefore, this admixture mapping study was unable to assess the genetic variation of these key effector molecules of innate immunity (55).

A joint analysis model of the local ancestry estimates and allele dosages at the most significant locus (rs17149618) had a very subtle impact on the EUR ancestry association (p = 4.11 × 10⁻⁵; OR [95%CI] = 0.46 [0.27–0.78]) (Table 2), suggesting that genetic variation at rs17149618 does not explain this admixture mapping peak. A fine-mapping assessing the genetic association of SNP variants of the 8p23.1 region was then performed to test other nearby variants as alternative explanations of the admixture mapping result (Figure 2). This analysis revealed two independent genetic variants with a suggestive association with sepsis protection, the most significant variant located within intron 1 of MFHAS1 (rs13249564; p = 9.94 × 10⁻⁴; OR = 0.65; 95%CI = 0.50–0.84) and an intergenic variant located between LINC00599 and MSRA (rs7820910; p = 1.42 × 10⁻³; OR = 0.49; 95%CI = 0.31–0.76) (Figure 2). A joint analysis of rs13249564 and local EUR ancestry did not alter the significance of the association of rs13249564 with sepsis (p = 2.20 × 10⁻⁴, OR [95% CI] = 0.59[0.45–0.78]), suggesting an independence between rs13249564 and EUR ancestry (Table 2). As a nested analysis, we then compared the two independent genetic variants between the population controls and the patients with sepsis caused by Gram-positive (N = 86, including those that were polymicrobial with mixed Gram-positive) or separately with those with Gram-negative bacteria (N = 115, including those that were polymicrobial with mixed Gram-negative). While both were significantly associated with Gram-negative bacterial sepsis (rs13249564, p = 1.39 × 10⁻³, OR [95% CI] = 0.50[0.33–0.77]; and rs7820910, p = 0.028, OR [95%CI] = 0.46[0.23–0.92]), none of the two was associated with Gram-positive bacterial sepsis (rs13249564, p = 0.113, OR [95% CI] = 0.70[0.45–1.09]; and rs7820910, p = 0.364, OR [95% CI] = 0.74[0.39–1.42]).

Furthermore, to evidence if a selective sweep was embedded in the admixture mapping region, an iSAFE scan was performed. The iSAFE scores in NAF ranged from 4.55 × 10⁻⁴ to 0.04, whereas in the IBS scan the results ranged from 2.05 × 10⁻⁴ to 0.104. Thus, the iSAFE scores in the region were highest for the IBS, with the top-ranking variants corresponding to positions 8.3–8.6 Mb (iSAFE score≥0.09), near PRAG1 and Claudin 23 (CLDN23) coding genes and telomeric to MFHAS1 (Figure 3; Table 3).

When we explored the potential functional implications of the prioritized variant in the fine-mapping study, rs13249564, we observed high evidence that it is an important regulatory variant, featuring DNase QTL and expression QTL (eQTL) in several cell lines (Supplementary Table 3). A RegulomeDB score of 0.609 and a RegulomeDB category of 4 were determined for it, indicating a weak effect in a DNase I hypersensitive site and transcription factor binding site, respectively. We found that this variant has high CellulAr dePendent dEactivating (CAPE) scores for DNase QTL and eQTL in several cell lines (Supplementary Table 3). According to GTEx results, rs13249564 is an eQTL in 13 different tissues (Supplementary Table 3), in addition to 4,836 significant Single-Tissue eQTLs and 126 significant Single-Tissue splicing QTL for MFHAS1 in all tissues. Likewise, there was also evidence of long-distance chromatin interactions among nearby candidate genes of the region, such as ERI1, CLDN23 and Tankyrase (TNKS) genes in lymphoblastoid cells. On the other hand, some of its 44 proxies (r² > 0.8) (Supplementary Figure 3) also have high CAPE scores for DNase QTL and eQTL in several cell lines, and many were significant eQTLs across different tissues in GTEx (Supplementary Table 4). Epigenome imputation, using peaks from H3K4me1 and H3K4me3 (as enhancers and promoters, respectively), and using peaks from H3K27ac and H3K9ac (as enhancers and promoters, respectively), revealed marks linked to most of the variants. Finally, we also annotated the functional implications of rs7820910, which also passed the suggestive threshold in the fine-mapping study (Figure 2). However, it lacked functional activity evidence (Supplementary Table 3).

The expression analysis of the genes mapping to the region 8p23.1 associated with sepsis revealed that besides MFHAS1 and CLDN23, many of the top-most significant differential expression (q-value <0.001) between sepsis cases and controls were concentrated in defensin genes (DEFA4 was the top ranked) (Supplementary Table 5). The results were similar in the three gene expression datasets comparing untreated and septic patients (GSE32707), healthy controls with those of septic shock patients (GSE57065), or with those from patients with sepsis (GSE28750).

In summary, we show the first admixture mapping study of sepsis. This allowed us to assess the association of local EUR ancestry with a protective effect against sepsis and suggested a sepsis risk association with local NAF ancestry at the 8p23.1 region, which harbors a novel promising genetic locus for sepsis susceptibility.