Further information and requests for reagents may be directed to, and will be fulfilled by the corresponding author Lluis Quintana-Murci ( quintana@pasteur.fr ).

We recruited 100 healthy, male donors of self-reported European descent (EUB) and 100 of self-reported African descent (AFB), all living in Belgium, at the Center for Vaccinology (CEVAC) of Ghent University Hospital (Ghent, Belgium). Samples were collected after written informed consent had been obtained, and the study was approved by the local ethics committee (Ethics Committee of the Ghent University), the Ethics Board of Institut Pasteur (EVOIMMUNOPOP-281297) and the relevant French authorities (CPP, CCITRS and CNIL). Inclusion was restricted to donors between 19 and 50 years of age, nominally healthy at the time of sample collection. A case report form was obtained for all donors, including information on vital sign measurements, medication, medical history and travel. No overrepresentation of any particular disease was observed relative to official report statistics published by the World Health Organization or in epidemiological studies. Serological testing was performed for all donors at the CEVAC, and those with serological signs of past or ongoing infection with human immunodeficiency virus (HIV), hepatitis B virus (HBV) or hepatitis C virus (HCV) were excluded.

The 200 subjects studied were genotyped for a total of 4,301,332 SNPs on the Illumina HumanOmni5-Quad BeadChips. Whole-exome sequencing was carried out for the same individuals with the Nextera Rapid Capture Expanded Exome kit, on the Illumina HiSeq 2000 platform, with 100-bp paired-end reads. This kit delivers 62 Mb of genomic content per individual, including exons, untranslated regions (UTR), and microRNAs.

Genomic DNA was extracted from the CD14-negative cell fraction (i.e. non-monocyte cells) by a standard phenol/chloroform protocol followed by ethanol precipitation. The DNA was quantified by Nanodrop spectrometry and with the Quant-iT PicoGreen dsDNA Assay Kit.

RNA was obtained from 978 of the 1000 samples, and was sequenced on an Illumina HiSeq2000. The quality and quantity of all samples was reassessed before sequencing. Samples were then randomized before library preparation in order to obtain a balanced number of samples across ethnicity and cellular conditions per sequencing batch/lane/machine/index. Standard reagents were used for transcriptome sequencing: TruSeq RNA Sample Prep Kit v2 for mRNA library construction, TruSeq SR Cluster Kit v3-HS for cluster generation and TruSeq SBS kit v3-HS for sequencing. We pooled six samples per lane to generate outputs of around 30 million 101-bp single-end reads per sample (ranging from 27.7 to 94.8 million reads, mean 34.4) ( Figure S2 A).

Total RNA was extracted with the Nucleospin miRNA kit from Macherey Nagel, including the enzymatic digestion of genomic DNA. Extractions were performed in batches of 30 samples (i.e. 5 conditions for 3 Africans and 3 Europeans), and RNA quality and quantity were assessed with a Nanodrop spectrometer and the Agilent Bioanalyzer RNA 6000 nano kit. We generated a final set of 978 samples from the 200 donors fulfilling the quality and quantity criteria (RIN > 7, quantity > 2.5 μg) for high-throughput RNA-sequencing, including 200, 188, 197, 193 and 200 samples for the non-simulated, LPS, Pam 3 CSK 4 , R848 and IAV conditions, respectively.

Monocytes were exposed to five different conditions for 6 hr, in order to capture transcriptional signatures from both an early response and the beginning of a late response, i.e., an “intermediate response” (). The choice of this time point was also based on a pilot study on the kinetics of gene expression of several key inflammatory and antiviral response genes (IL1A, IL23A, IL6, IL8, TNF, IRF1 and STAT2) upon immune activation. Our results showed that 6 hr of stimulation was the best time point to capture simultaneously expression signals from early, intermediate and late response genes, with respect to other time points at 2, 4, 8 and 24 hr (data not shown). One monocyte flask was left untreated as a baseline control, while the others were each exposed to one of four different immune stimuli. These stimuli included synthetic ligands specifically activating three Toll-like receptor (TLR) signaling pathways: 1 ng/ml ultrapure LPS from E. coli, 0.2 μg/ml synthetic triacylated lipoprotein PamCSK, and 0.3 μg/ml imidazoquinoline compound R848. Monocytes were also infected with strain A/USSR/90/1977(H1N1) of the human seasonal influenza A virus (IAV) at a MOI = 1, and IAV particles were produced as previously described (). After stimulation, cells were collected by centrifugation, lysed in a guanidinium thiocyanate solution provided in the Nucleospin miRNA kit, according to the manufacturer’s instructions, and stored at −80°C until RNA extraction. Cellular assays were performed per batch of 30 samples from 6 individuals, including 3 Africans and 3 Europeans, across all 5 conditions.

Purity and cell death of the isolated monocytes were assessed for all donors on a fraction of 10 5 CD14 + monocytes stained, according to the manufacturer’s instructions, with fluorescent APC-conjugated anti-CD14 antibodies and propidium iodide, respectively. Samples were then analyzed on a MACSQuant Analyzer 10 benchtop flow cytometer (Miltenyi Biotec) and using FlowJo vX.0.6 software. The mean values obtained for all samples were 96.8% for monocyte purity and 2.1% for initial cell death rates.

For each donor, 300 × 10 6 PBMCs were thawed, washed twice and resuspended in pre-warmed RPMI-1640 Glutamax medium, supplemented with 10% FCS and penicillin/streptomycin (complete medium). Monocytes were then positively selected with magnetic CD14 microbeads, according to the manufacturer’s instructions. The number of monocytes was determined with a Kova Glasstic Slide 10 with a grid in the presence of trypan blue. For each donor, 30 × 10 6 monocytes were split between five 25 cm 2 non-treated flasks (i.e. one flask per condition and five conditions per donor), each containing 6 × 10 6 monocytes in 9 ml of complete medium. Monocytes were allowed to rest for one hour at 37°C under 5% CO 2 before stimulation.

For each participant, we collected 300 ml of whole blood into anticoagulant EDTA-blood collection tubes and peripheral blood mononuclear cells (PBMCs) were isolated on Ficoll-Paque density gradients. PBMCs were frozen in 90% fetal calf serum (FCS) and 10% dimethyl sulfoxide, at a density of 50 × 10 6 PBMCs/ml and transported in dry shipper from CEVAC to Institut Pasteur. Vials were then cryopreserved in liquid nitrogen until use.

Quantification and Statistical Analysis

RNA-Sequencing Analysis Kim et al., 2013 Kim D.

Pertea G.

Trapnell C.

Pimentel H.

Kelley R.

Salzberg S.L. TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. 3 CSK 4 (EUB: 100, AFB: 96), 191 R848 (EUB: 98, AFB: 93), and 199 IAV samples (EUB: 99, AFB: 100). Reads were assessed for multiple quality metrics, including number of reads, nucleotide distribution and sequencing quality, and the last base of all reads was trimmed due to a fall in sequencing quality. RNA reads were then mapped onto the human GRCh37 genome with TopHat (), resulting in the successful mapping of 89.9% of reads per sample on average (minimum 67.3%; maximum 93.7%). We used the RSeQC package to assess the alignment of reads with various genomic features, GC content, and gene body coverage ( Figures S2 B–S2E). Samples with uneven gene body coverage were found to be more likely to be outliers. We used gene body coverage regularity as an indicator of library quality, removing eight samples due to irregular gene body coverage. The remaining 970 samples were used for subsequent analyses and consisted of 200 non-stimulated (EUB: 100, AFB: 100), 184 LPS (EUB: 96, AFB: 88), 196 PamCSK(EUB: 100, AFB: 96), 191 R848 (EUB: 98, AFB: 93), and 199 IAV samples (EUB: 99, AFB: 100). Trapnell et al., 2012 Trapnell C.

Roberts A.

Goff L.

Pertea G.

Kim D.

Kelley D.R.

Pimentel H.

Salzberg S.L.

Rinn J.L.

Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. L o g 2 ( 1 + F P K M i j ) = α + β c o n d + β p o p + ∑ c β c C o v a r i a t e c + ∑ k β k , j b a t c h + ε i j

where α is the intercept, β c o n d and β p o p are the fixed effect of the condition and population on sample j, β c are the fixed effect of continuous covariates on sample j, the β k , j b a t c h are the random effects of batch covariate k, on sample j, and ε i j are the residuals. Cufflinks/CuffDiff (v2.0.2) () was used to quantify expression levels in FPKM (fragments per kilobase of transcript per million mapped reads) for each annotated transcript of the genome in Ensembl (v.70), and FPKM values for which Cufflinks returned FAIL status (< 0.5% of quantified transcripts) were set to missing values. Gene expression data were filtered to remove genes with low levels of expression (mean FPKM < 1 in all conditions) and their quality was checked by principal component analysis (PCA). PCA captured differences between conditions and populations on the first two axes, but we tested for additional causes of technical variability, by fitting, for each gene, a mixed model of gene expression as a function of condition, population, and technical covariates, including total RNA concentration, RIN, percentage of high-quality bases (Q30), mean GC content, library concentration, 5′/3′ coverage bias (measured as the mean difference in coverage between the 5′ and 3′ ends of the gene) as continuous covariates, and date of experiment, library preparation batch, sequencing batch, sequencer used, sequencing index, and sequencing lane as putative batch effects. Putative batch effects were modeled as random effects to prevent the loss of degrees of freedom, whereas all other covariates (condition, population and continuous covariates) were included as fixed effects, giving the following model for gene i and sample j:whereis the intercept,andare the fixed effect of the condition and population on sample j,are the fixed effect of continuous covariates on sample j, theare the random effects of batch covariate k, on sample j, andare the residuals. Johnson et al., 2007 Johnson W.E.

Li C.

Rabinovic A. Adjusting batch effects in microarray expression data using empirical Bayes methods. The proportions of genes affected by each factor are reported in Figure S2 F for various levels of explained variance. We observed that GC content, 5′/3′ bias, date of the experiment and library batch were among the strongest confounding factors, and accordingly corrected the data for these factors before analysis, following the pipeline detailed in Figure S3 . First, we adjusted the data for GC content and 5′/3′ bias using linear models. Then, we imputed missing values by K-nearest neighbor imputation, and adjusted for experiment date and library batch by sequentially running ComBat () for each batch effect, with condition and population as covariates. Batch-corrected gene expression levels, in FPKM, were then recalculated from the adjusted transcript level estimates. Refitting our linear mixed model confirmed that correction was satisfactory for most of the technical covariates ( Figure S2 F).

Assessment of Technical and Biological Variability −16; The reproducibility of our RNA-Seq experiments was assessed by performing technical and biological replicates on seven independent donors (4 AFB and 3 EUB) across the five experimental conditions. We showed that (i) the coefficients of variation of genes within technical replicates were consistently, and significantly, smaller in magnitude and less variable than those within biological replicates (Wilcoxon Rank-Sum Test, p < 10 Figure S2 G), and (ii) technical replicates exhibit higher correlation coefficients (r) between samples with respect to the distribution of r values calculated from pairwise comparisons between biological replicates ( Figure S2 H).

Modules of Correlated Genes Langfelder and Horvath, 2008 Langfelder P.

Horvath S. WGCNA: an R package for weighted correlation network analysis. Modules of genes presenting correlated expression patterns, extracted from log-transformed FPKM data, were defined by weighted correlation network analysis (WGCNA) (). In our setting of immune response activation, this analysis detects modules of correlated genes that can reflect either shared regulation by common transcription factors, or regulation by independent transcription factors with similar patterns of activation upon stimulation. Tukey’s biweight correlation was used as a measure of gene relatedness to reduce the influence of outliers, and correlations were measured across all 970 samples. The scale-free topology of the networks was assessed for various values of the β shrinkage parameter, according to WGCNA user manual, and the default value of β = 6 appeared to give a satisfactory fit to scale-free topology. Signed clustering of genes (grouping only positively correlated genes) was used to simplify the interpretation of the extracted modules. We also found that varying the level of shrinkage (β = 5 or 6) or the depth of the clustering (deepsplit parameter set to 3 or 4) had only a mild impact on the number of clusters or the enrichments obtained, confirming the robustness of these analyses. Roider et al., 2009 Roider H.G.

Manke T.

O’Keeffe S.

Vingron M.

Haas S.A. PASTAA: identifying transcription factors associated with sets of co-regulated genes. Thomas-Chollier et al., 2011 Thomas-Chollier M.

Hufton A.

Heinig M.

O’Keeffe S.

Masri N.E.

Roider H.G.

Manke T.

Vingron M. Transcription factor binding predictions using TRAP for the analysis of ChIP-seq data and regulatory SNPs. Mathelier et al., 2014 Mathelier A.

Zhao X.

Zhang A.W.

Parcy F.

Worsley-Hunt R.

Arenillas D.J.

Buchman S.

Chen C.Y.

Chou A.

Ienasescu H.

et al. JASPAR 2014: an extensively expanded and updated open-access database of transcription factor binding profiles. For each module, we used PASTAA () to identify transcription factor binding site motifs overrepresented within the annotated proximal promoters of the genes within each module. We first defined the proximal promoter region for each gene as the region extending 200 bp on either side of the transcription start site (TSS) of the most abundant transcript on the basis of Cufflinks FPKM estimates. We then used the transcription factor affinity prediction (TRAP) method () to measure the binding affinities of each transcription factor present in the Jaspar core vertebrate database () with the proximal promoters of the 12,578 expressed genes, and these affinities were then used as the input for PASTAA enrichment analysis. We reported only enrichments significant at a false discovery rate (FDR) of 0.05 with a fold-change (i.e. observed/expected) greater than 1.2. For each module, we represent the transcription factor binding sites with the highest value for the lower limit of the odds ratio confidence interval.

Differential Expression Analysis 2 fold change between populations – | log 2 F C p o p | – greater than 0.2 and at FDR < 0.05. We then calculated the fold-change in expression after stimulation relative to the basal state, and used t tests to determine whether there was a differential response. Population differential response genes (popDRGs) were then defined as popDEGs for which there was a differential response between populations under stimulated conditions, at FDR < 0.05 (i.e. the transcriptional response to treatment, relative to the basal state, differed between populations), resulting in a larger difference in expression after stimulation. | log 2 F C p o p s t i m u l a t e d | > | log 2 F C p o p b a s a l |

Differential expression was assessed directly from log-transformed FPKM, using t tests for each condition. FDR was then calculated jointly for all conditions, with the R package fdrtool, and genes differentially expressed between populations (popDEGs) were defined as genes presenting an absolute logfold change between populations –– greater than 0.2 and at FDR < 0.05. We then calculated the fold-change in expression after stimulation relative to the basal state, and used t tests to determine whether there was a differential response. Population differential response genes (popDRGs) were then defined as popDEGs for which there was a differential response between populations under stimulated conditions, at FDR < 0.05 (i.e. the transcriptional response to treatment, relative to the basal state, differed between populations), resulting in a larger difference in expression after stimulation.

Gene Ontology Enrichment Analysis Ashburner et al., 2000 Ashburner M.

Ball C.A.

Blake J.A.

Botstein D.

Butler H.

Cherry J.M.

Davis A.P.

Dolinski K.

Dwight S.S.

Eppig J.T.

et al. The Gene Ontology Consortium

Gene ontology: tool for the unification of biology. All Gene Ontology (GO) enrichment analyses were performed with GOSeq package (), using the default settings, with the 12,578 expressed genes as the background set. Only enrichments significant at FDR of 0.05 and with a fold-change (i.e. observed/expected) greater than 1.2 are reported.

SNP Genotyping Data Analysis Chang et al., 2015 Chang C.C.

Chow C.C.

Tellier L.C.

Vattikuti S.

Purcell S.M.

Lee J.J. Second-generation PLINK: rising to the challenge of larger and richer datasets. −3 in AFB or EUB populations (N = 4,007). After quality-control filtering, we retained a total of 3,489,108 SNPs. The SNP call rate for the 200 individuals was 99.8% on average, ranging from 98.89% to 99.99%. No evidence was found for 2nd-degree cryptic relatedness (kinship coefficient > 0.07) in KING ( Manichaikul et al., 2010 Manichaikul A.

Mychaleckyj J.C.

Rich S.S.

Daly K.

Sale M.

Chen W.M. Robust relationship inference in genome-wide association studies. Using PLINK v1.9 (), we removed SNPs that: (i) were typed with probes mapping to several genomic locations (N = 12,440), (ii) presented a poor genotype clustering (GenTrain score < 0.35; N = 809), (iii) had the same chromosomal position as another SNP in dbSNP b138 (N = 6,968), (iv) were not reported in dbSNP b138 (N = 5,311), (v) presented a call rate < 95% (N = 79,310), (vi) were monomorphic in our sample (N = 652,385), (vii) were located on the sex chromosomes (N = 50,994), and (viii) presented a Hardy-Weinberg p < 10in AFB or EUB populations (N = 4,007). After quality-control filtering, we retained a total of 3,489,108 SNPs. The SNP call rate for the 200 individuals was 99.8% on average, ranging from 98.89% to 99.99%. No evidence was found for 2-degree cryptic relatedness (kinship coefficient > 0.07) in KING (), or for sex mismatch, for any of the individuals. Two AFB individuals presented an excess of heterozygosity (< ± 3SD from the population average), as a result of their moderate levels of non-African ancestry, as estimated using ADMIXTURE.

Whole-Exome Data Analysis Li and Durbin, 2009 Li H.

Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. DePristo et al., 2011 DePristo M.A.

Banks E.

Poplin R.

Garimella K.V.

Maguire J.R.

Hartl C.

Philippakis A.A.

del Angel G.

Rivas M.A.

Hanna M.

et al. A framework for variation discovery and genotyping using next-generation DNA sequencing data. −3 in AFB or EUB populations (N = 4,510). The application of these quality-control filters resulted in the retention of 471,740 SNPs. Read-pairs were processed according to the GATK Best Practice recommendations. Read-pairs were first mapped onto the human GRCh37 genome with BWA v.0.7.7 (), and reads duplicating the start position of another read were marked as duplicates with Picard Tools v.1.94 ( http://broadinstitute.github.io/picard ). We used GATK v.3.2.2 () for base quality score recalibration (“BaseRecalibrator”), insertion/deletion realignment (“IndelRealigner”), and SNP and insertion/deletion discovery for each sample (“Haplotype Caller”). Individual variant files were combined with “GenotypeGVCFs” and filtered with “VariantQualityScoreRecalibration.” Individual coverage was 52.32 × on average, ranging from 33.84 to 100.59 ×, and individual breadth of coverage at 5 × was 92.42%, ranging from 83.5% to 95.0%. We removed those of the 540,990 SNPs obtained that: (i) were triallelic (N = 11,925), (ii) presented a call rate < 95% (N = 44,716), (iii) were located on the sex chromosomes (N = 8,369), and (iv) presented a Hardy-Weinberg p < 10in AFB or EUB populations (N = 4,510). The application of these quality-control filters resulted in the retention of 471,740 SNPs.

Imputation of Genome-wide SNP and Exome Data Before merging the Omni5 and exome datasets, we first checked genotype concordance for 169,406 SNPs common to the two platforms. We flipped alleles for 8,025 SNPs with incompatible allelic states, and removed 119 SNPs with alleles that remained incompatible after allele flipping. The total concordance rate was 99.66%. The concordance rates for each of the 200 individuals exceeded 99%, confirming an absence of errors during DNA sample processing. Of the 8,155 SNPs with discordance rates > 1%, 296 were due to C/G or A/T SNPs, and high genotype concordance between the two DNA typing technologies was restored by allele flipping. The remaining 7,881 SNPs were removed. The entire Omni5 and exome datasets (3,489,108 and 471,740 SNPs, respectively) were then merged, yielding a final concordance rate of 99.93%, for a total of 3,782,260 SNPs. Delaneau et al., 2013 Delaneau O.

Zagury J.F.

Marchini J. Improved whole-chromosome phasing for disease and population genetic studies. Howie et al., 2009 Howie B.N.

Donnelly P.

Marchini J. A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. Before imputation, we phased the data with SHAPEIT2 (), using 500 conditioning haplotypes, 50 MCMC iterations, 10 burn-in and 10 pruning iterations. SNPs and allelic states were then aligned with the 1,000 Genomes Project imputation reference panel (Phase 1 v3.2010/11/23). We removed 8,705 SNPs with identical positions in our data and in the reference panel but incompatible alleles, even after allele flipping, and 4,137 SNPs with C/G or A/T alleles. Genotype imputation was performed with IMPUTE v.2 (), considering 1-Mb windows and a buffer region of 1 Mb. Of the 38,098,530 SNPs obtained after imputation, we removed SNPs that: (i) presented an information metric below 0.8 (N = 18,085,215), (ii) had a duplicate (N = 59,914), (iii) presented a call rate < 90% (N = 329,910), and (iv) were monomorphic (N = 4,053). The final imputed dataset included 19,619,457 SNPs. 2 between true genotypes (i.e., obtained by Illumina array genotyping or exome sequencing) and imputed genotypes for the same SNPs (i.e., obtained by artificially removing genotyped SNPs from the data before imputation and then imputing them). In very good agreement with recent studies ( Auton et al., 2015 Auton A.

Brooks L.D.

Durbin R.M.

Garrison E.P.

Kang H.M.

Korbel J.O.

Marchini J.L.

McCarthy S.

McVean G.A.

Abecasis G.R. 1000 Genomes Project Consortium

A global reference for human genetic variation. To evaluate imputation accuracy, we estimated correlation coefficients Rbetween true genotypes (i.e., obtained by Illumina array genotyping or exome sequencing) and imputed genotypes for the same SNPs (i.e., obtained by artificially removing genotyped SNPs from the data before imputation and then imputing them). In very good agreement with recent studies (), the average correlation coefficient was 95.6% across all genotyped SNPs with information metric > 0.8 (93.6% for SNPs with MAF < 0.10 and 97.7% for SNPs with MAF > 0.10). This shows that our stringent quality filters ensure that only accurately imputed SNPs are analyzed.

eQTL Mapping Shabalin, 2012 Shabalin A.A. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. 10 p-values of eQTL > 0.95). For each gene, SNPs were considered “local,” likely cis-acting, if they were located less than 1Mb away from the gene transcription start or end site. They were otherwise considered to be trans-acting. eQTL mapping was performed separately for each population and condition, and false-positives due to outliers were prevented by discarding, from the analysis, eQTL associations that did not pass a p-value threshold of 10−3 for local eQTLs, and 10−5 for trans-eQTLs in Kruskal-Wallis rank tests. For expression quantitative trait loci (eQTL) mapping, only variants with a minor allele frequency (MAF) ≥ 0.05 in the population studied were retained in the analysis, resulting in a set of 10,278,745 SNPs (i.e., corresponding to the merged genotyping and imputed SNP dataset; 8,913,090 SNPs in Africans and 6,178,808 SNPs in Europeans,). We mapped eQTLs with the MatrixEQTL R package (). PC1 and PC2 of the genotype matrix were included in the model to account for possible population stratification. The inclusion of additional PC in the model (up to PC6) was tested and showed highly consistent results (i.e., correlation of -logp-values of eQTL > 0.95). For each gene, SNPs were considered “local,” likely cis-acting, if they were located less than 1Mb away from the gene transcription start or end site. They were otherwise considered to be trans-acting. eQTL mapping was performed separately for each population and condition, and false-positives due to outliers were prevented by discarding, from the analysis, eQTL associations that did not pass a p-value threshold of 10for local eQTLs, and 10for trans-eQTLs in Kruskal-Wallis rank tests. For both cis- and trans-eQTLs, FDR was computed by mapping eQTLs on 100 datasets with genotypes permuted within each population. We then kept, after each permutation, the most significant p-value per gene, across all conditions and populations. Finally, we computed the false discovery rate associated with each p-value threshold in cis or in trans, and subsequently selected the p-value threshold that provided a 5% FDR, leading to p = 7.67 × 10−7 and p = 2.7 × 10−12 for cis- and trans-eQTLs, respectively. For local eQTLs, we report only the SNP at which the strongest association was observed (i.e., eQTL peak-SNP). When multiple SNPs in perfect LD fell within the peak, only one SNP is reported. eQTLs for which the eQTL peak-SNP had an allelic effect size (|β eQTL |) below 0.2 were discarded from further analysis. We next mapped fold-changes between the basal and stimulated states using MatrixEQTL, and defined response eQTLs (reQTLs) as stimulated eQTLs associated to a significant difference in response to stimulation (p < 10−3, | β e Q T L s t i m | > | β e Q T L b a s a l | ). For trans-eQTLs, we reported, within each 1Mb window, the SNP for which we observed both the largest number of trans-associated genes and strongest p-value of association. Furthermore, for each trans-eQTL that passed genome-wide significance at p = 2.7 × 10−12 (FDR of 5%), we performed a SNP-based analysis to identify genes regulated in trans by the eQTL at a Bonferroni p < 0.05, correcting for the 12,578 genes tested within the condition where the eQTL was found.

Population Differences Attributable to Genetics 2 > 0.5) with the eQTL peak-SNP, in the population where the eQTL was discovered and fine map the eQTL signal by fitting across populations the following linear model: expression j = α + β . SNP j + γ . Pop j + ε j

where SNP j is the genotype of the individual j for the variant under study, Pop j is a binary variable indicating the population origin (0 for Europeans and 1 for Africans), and ε j is a random, normally distributed residual. In this model, α is the intercept, β reflects the effect of the derived allele of the SNP on gene expression, and γ estimates the fold change in expression between populations observed for individuals with identical genotype (i.e. gene expression differences that are not explained by genetics). We next focused on the SNP showing the strongest association p-value with gene expression across populations, and estimated the difference in population expression that is attributable to the SNP as: FC SNP = FC pop − γ ’

with γ’ representing γ set to ensure that the ratio of FC SNP /FC pop is between 0 and 1, i.e. γ’ = 0, if the sign of γ differs from that of FC pop ; γ’ = FC pop , if |γ| > |FC pop |; and γ’ = γ otherwise. The percentage of population differences in expression that is attributable to genetics is then given by the ratio FC SNP /FC pop . To estimate the fraction of population differences in gene expression that can be attributed to genetic variants, we used a two-step strategy. First, we consider the set of all SNPs in LD (r> 0.5) with the eQTL peak-SNP, in the population where the eQTL was discovered and fine map the eQTL signal by fitting across populations the following linear model:where SNPis the genotype of the individual j for the variant under study, Popis a binary variable indicating the population origin (0 for Europeans and 1 for Africans), and εis a random, normally distributed residual. In this model, α is the intercept, β reflects the effect of the derived allele of the SNP on gene expression, and γ estimates the fold change in expression between populations observed for individuals with identical genotype (i.e. gene expression differences that are not explained by genetics). We next focused on the SNP showing the strongest association p-value with gene expression across populations, and estimated the difference in population expression that is attributable to the SNP as:with γ’ representing γ set to ensure that the ratio of FC/FCis between 0 and 1, i.e. γ’ = 0, if the sign of γ differs from that of FC; γ’ = FC, if |γ| > |FC|; and γ’ = γ otherwise. The percentage of population differences in expression that is attributable to genetics is then given by the ratio FC/FC

Defining Population-Specific eQTLs We aimed at distinguishing population specific eQTLs (i.e., SNPs present at similar frequencies in both populations but having an effect on gene expression in one population only) from eQTLs detected in one population only due to population differences in allelic frequencies. To do so, we first focused on the 1,109 genes associated with an eQTL (including 363 genes associated with a reQTL) where all SNPs in LD (r2 > 0.5) with the eQTL peak-SNP were present at frequency > 5% in both populations. We then tested these eQTLs for replication at a relaxed threshold of p < 0.05 across all SNPs at the locus, to decrease the false negative rate, and focused on the 127 genes for which the eQTL was not replicated (including 28 genes with a reQTL). 2 > 0.5 with the eQTL peak-SNP), the following linear model: expression j = α + β . SNP j + γ . Pop j + δ . SNP j ∗ Pop j + ε j

where SNP j is the genotype of the individual j for the variant under study, Pop j is a binary variable indicating the population origin (0 for Europeans and 1 for Africans), and ε j is a random, normally distributed residual. In this model, β reflects the effect of the derived allele of the SNP on gene expression, γ estimates the fold change in expression between populations observed for individuals with identical genotype, and δ captures the differences in eQTL effect size between populations. Such a model allows to test for a difference in eQTL effect size between populations by testing the null hypothesis, δ = 0 (interaction test). Finally, we considered as population-specific, eQTLs whose effect size was significantly different between populations. To do so, we fit, for each SNP at the locus (r> 0.5 with the eQTL peak-SNP), the following linear model:where SNPis the genotype of the individual j for the variant under study, Popis a binary variable indicating the population origin (0 for Europeans and 1 for Africans), and εis a random, normally distributed residual. In this model, β reflects the effect of the derived allele of the SNP on gene expression, γ estimates the fold change in expression between populations observed for individuals with identical genotype, and δ captures the differences in eQTL effect size between populations. Such a model allows to test for a difference in eQTL effect size between populations by testing the null hypothesis, δ = 0 (interaction test). 2 > 0.5) with the eQTL peak-SNP presented a significant interaction p-value. P i n t e r a c t i o n ( l o c u s ) = Max s n p ∈ l o c u s P i n t e r a c t i o n ( s n p ) .

To be conservative and to account for the uncertainty in detecting the causal variant at the eQTL, we considered an eQTL as population specific if all SNPs in LD (r> 0.5) with the eQTL peak-SNP presented a significant interaction p-value. We then considered eQTLs (or reQTLs) as being population specific when the interaction p-value at the locus was lower than 10−3 (corresponding to FDR < 0.01), leading to a final set of 16 population-specific eQTLs (including 5 reQTLs).

Regulatory Elements and Transcription Factor Binding Sites Zerbino et al., 2015 Zerbino D.R.

Wilder S.P.

Johnson N.

Juettemann T.

Flicek P.R. The ensembl regulatory build. Ernst and Kellis, 2015 Ernst J.

Kellis M. Large-scale imputation of epigenomic datasets for systematic annotation of diverse human tissues. Kundaje et al., 2015 Kundaje A.

Meuleman W.

Ernst J.

Bilenky M.

Yen A.

Heravi-Moussavi A.

Kheradpour P.

Zhang Z.

Wang J.

Ziller M.J.

et al. Roadmap Epigenomics Consortium

Integrative analysis of 111 reference human epigenomes. Ernst and Kellis, 2015 Ernst J.

Kellis M. Large-scale imputation of epigenomic datasets for systematic annotation of diverse human tissues. Regulatory features were extracted from Ensembl Regulatory Build v80 (), which contains regulatory element predictions based on open chromatin regions and histone marks from ENCODE and the Roadmap Epigenomics datasets (). SNPs overlapping a regulatory element were then classified into four categories: promoter, promoter flanking, enhancer, and CTCF binding sites. Similarly, ENCODE uniformly processed transcription factor binding site (TFBS) clusters (V3) () were downloaded from UCSC, and their overlap with the physical position of all SNPs was determined. We then used Fisher’s exact test to assess the eQTL enrichment of specific TFBS or regulatory elements, considering the peak-SNP at each locus, or a randomly selected SNP if multiple SNPs in perfect LD were found. All SNPs with a MAF ≥ 0.05 located less than 1Mb away from an expressed gene were used to constitute the background set. In each condition (or combination of conditions), only the TFBS with the highest values for the lower limit of the odds ratio confidence intervals are reported.

Quantification of Allelic Imbalance Li and Durbin, 2009 Li H.

Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. van de Geijn et al., 2015 van de Geijn B.

McVicker G.

Gilad Y.

Pritchard J.K. WASP: allele-specific software for robust molecular quantitative trait locus discovery. Li, 2011 Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. For the quantification of allele-specific imbalance, we focused on exonic SNPs genotyped as heterozygous in our exome data, excluding SNPs with discordant genotypes in the Omni5 data. We used BWA mem (v.0.7.7) () to remap RNA-seq reads onto the hg19 genome for all 970 samples, and extracted all reads aligned with a genetic variant. We reduced mapping bias, by using WASP () to exclude reads overlapping with known variants (based on dbSNP138) likely to alter the read mapping location. Briefly, for each read overlapping one or more dbSNP variants, WASP creates alternative reads consisting of all possible combinations of reads given these SNPs. It then remaps the alternative reads to the genome, and keeps the original read only if all alternative versions of the read map to the same position. Finally, SAMtools mpileup (), with option -d 10000, was used to count the number of reads mapping to each allele at heterozygous loci. The allelic ratio (AR) was defined for each site as the proportion of minor alleles among all reads, and the allelic imbalance (AI) was defined as the absolute deviation from a balanced ratio of 0.5 (i.e. AI = |AR-0.5|).

aseQTL and asrQTL Mapping We mapped allele-specific expression QTLs (aseQTLs), by estimating the allelic ratio on the subset of eQTL-genes with sufficient expression coverage at heterozygous exonic SNPs (N ≥ 10 reads) in at least five individuals of each eQTL genotype (heterozygous/homozygous). We extracted the phase information between the strongest local eQTL and the exonic SNP, and tested the correlation between the AR and the phased eQTL genotypes (coded 0 for homozygous, and ± 1 for heterozygotes with variants in phase or in the opposite phase), in a gene-, condition- and population-specific manner. Each exonic SNP was considered as an independent observation. Similarly, allele-specific response QTLs (asrQTLs) were mapped by assessing the correlation between the phased reQTL genotypes and the change in AR at the exonic site after stimulation. The power to detect aseQTLs was computed for various eQTL effect sizes |β|, number of observations n and number of reads per exonic SNP N. We assumed the same number of observations for heterozygous and homozygous genotypes at the eQTL, and equal coverage across all exonic SNPs. Power was then computed for a standard t-test assuming a mean allelic ratio N alternative /N reference of 0.5 in homozygous individuals and 2β/(1+ 2β) in heterozygous individuals. Residual variance was set to 0.25/N to match that of a binomial distribution with parameters (0.5, N).

ASE Analysis at the Individual Level To ensure sufficient power when exploring ASE within single individuals, we considered a higher coverage of heterozygous exonic SNPs (N ≥ 30 reads), and used a binomial test to evaluate allelic imbalance. We also excluded sites at which one allele accounted for less than 2% of the reads or less than 3 reads in total, as such sites might be subject to genotyping errors or systematic mapping biases. The FDR was first calculated across all SNPs, individuals and conditions, using fdrtool, and ASE was defined as the combination of significance at FDR = 0.05 and an absolute log 2 fold change of expression between alleles of more than 0.2 (|log 2 (N alt /N ref )| > 0.2). For each significant ASE event in stimulated conditions, we checked for differences in allelic imbalance relative to the non-stimulated condition, and defined allele-specific response as the subset of ASE displaying significantly higher allelic imbalance (p < 10−3, Fisher’s exact test) after stimulation with respect to the basal state. Finally, we used simulations to evaluate FDR among the set of genes with at least one ASE/ASR event. We generated 1,000 null datasets, by randomly reassigning reads to the alternative and reference allele with equal probability, and estimated the number of genes with at least one significant ASE or ASR event at each p-value threshold. We then computed FDR as the ratio of the average number of genes with ASE in our resampling to the observed number of genes with ASE at the same p-value threshold.

ASE Enrichment in Rare Coding Variants McLaren et al., 2010 McLaren W.

Pritchard B.

Rios D.

Chen Y.

Flicek P.

Cunningham F. Deriving the consequences of genomic variants with the Ensembl API and SNP Effect Predictor. For each exonic SNP for which we quantified ASE, we used Variant Effect Predictor (VEP) () with Ensembl v.70 Transcript Annotation to identify the set of transcripts overlapping the variant, and Cufflinks FPKM to identify the most strongly expressed overlapping transcript in the individual/condition concerned. VEP annotation was then used to classify variants, according to the most abundant transcript, as synonymous (synonymous_variant/ non_coding exon_variant), missense (missense_variant) or nonsense (stop_gained, stop_lost, initiator_codon_variant). Enrichment in rare coding variants was then assessed using Fisher’s exact test comparing each category with synonymous variants.

Natural Selection Analysis: Neutrality Statistics ST and iHS, which detect signals of population-specific positive selection, i.e., mutations that provided a selective advantage to a specific human population. F ST measures population differentiation by comparing the variance of allele frequencies within and between populations ( Holsinger and Weir, 2009 Holsinger K.E.

Weir B.S. Genetics in geographically structured populations: defining, estimating and interpreting F(ST). ST is a population pairwise comparison, we derived a directional F ST , equal in absolute value to the pairwise F ST but with a positive sign if the derived allele was more frequent in the population studied, and a negative sign otherwise. This enables to distinguish selection events that likely occurred in Africans from those that likely occurred in Europeans. The integrated haplotype score (iHS) measures the degree of extended haplotype homozygosity of the putatively selected allele over that of the putatively neutral allele ( Voight et al., 2006 Voight B.F.

Kudaravalli S.

Wen X.

Pritchard J.K. A map of recent positive selection in the human genome. We used two metrics, Fand iHS, which detect signals of population-specific positive selection, i.e., mutations that provided a selective advantage to a specific human population. Fmeasures population differentiation by comparing the variance of allele frequencies within and between populations (), as local positive selection tends to increase allele frequency differences between populations. As Fis a population pairwise comparison, we derived a directional F, equal in absolute value to the pairwise Fbut with a positive sign if the derived allele was more frequent in the population studied, and a negative sign otherwise. This enables to distinguish selection events that likely occurred in Africans from those that likely occurred in Europeans. The integrated haplotype score (iHS) measures the degree of extended haplotype homozygosity of the putatively selected allele over that of the putatively neutral allele (), as the long-range associations of the selected mutation with neighboring SNPs are not disrupted by recombination. ST and iHS. The CSS was designed to identify variants with both a higher derived allele frequency in one population (positive value of directional F ST ), and a longer haplotype length around the derived allele of the variant in that population (characterized by a negative iHS value). It was computed for all variants with derived allele frequency 0.2 ≤ DAF ≤ 0.95 from genome-wide ranks of both directional F ST ( R F s t ) and iHS ( R i H S ) , attributing the highest rank to positive values of F ST and negative values of iHS, respectively. We defined the CSS as following: C S S = r a n k ( R F s t . R i H S ) N o b s

with N o b s being the total number of variants with 0.2 ≤ DAF ≤ 0.95 in the population studied. CSS ranges from 0 to 1, and increases with the strength of positive selection targeting the derived allele. Furthermore, we used a composite selection score (CSS) allowing to capture signals of recent, strong selective events, by combining Fand iHS. The CSS was designed to identify variants with both a higher derived allele frequency in one population (positive value of directional F), and a longer haplotype length around the derived allele of the variant in that population (characterized by a negative iHS value). It was computed for all variants with derived allele frequency 0.2 ≤ DAF ≤ 0.95 from genome-wide ranks of both directional Fand iHSattributing the highest rank to positive values of Fand negative values of iHS, respectively. We defined the CSS as following:withbeing the total number of variants with 0.2 ≤ DAF ≤ 0.95 in the population studied. CSS ranges from 0 to 1, and increases with the strength of positive selection targeting the derived allele. Chen et al., 2010 Chen H.

Patterson N.

Reich D. Population differentiation as a test for selective sweeps. Finally, we used the cross-population composite likelihood ratio score, XP-CLR, a region-based metric detecting extended regions where the allele frequencies of multiple contiguous markers are distorted from the prediction under neutrality (). XP-CLR detects classical selective sweeps as well as selection events on pre-existing alleles (standing variation). XP-CLR was scored every 2,000 bp, using windows of 0.2 cM and downsampling to 200 SNPs per window.

Enrichment Tests for Natural Selection Signals To map selection signals at haplotypes containing eQTLs, we determined, for each statistic (iHS, F ST , or CSS), the strongest signal of selection on derived alleles of all SNPs in high LD (r2 > 0.8) with the eQTL peak-SNP. To assess significance, we then compared, for each population and condition, the mean of these values across all eQTLs/reQTLs, with the expected distribution obtained from resampling 10,000 sets of random SNPs matched for MAF (using bins of MAF of 0.05) and the number of SNPs in LD (r2 > 0.8, using bins of 0-2, 3-5, 6-10, 11-20, 21-50, and > 50 SNPs in LD). Similarly, for XP-CLR, we compared the mean of XP-CLR scores at eQTLs/reQTLs (considering the region that contains the eQTL peak-SNP), to the expected distribution obtained from resampling 10,000 sets of random SNPs matched for MAF and LD patterns.

Detection of Candidate eQTLs under Selection ST and iHS at the genome-wide level, focusing on signals consistent with selection on derived alleles, within each population separately. To support the adaptive nature of candidate eQTLs, we computed neutral p values for each statistic using simulations based on validated demographic models of Africans and Europeans ( Grossman et al., 2013 Grossman S.R.

Andersen K.G.

Shlyakhter I.

Tabrizi S.

Winnicki S.

Yen A.

Park D.J.

Griesemer D.

Karlsson E.K.

Wong S.H.

et al. 1000 Genomes Project

Identifying recent adaptations in large-scale genomic data. Grossman et al., 2013 Grossman S.R.

Andersen K.G.

Shlyakhter I.

Tabrizi S.

Winnicki S.

Yen A.

Park D.J.

Griesemer D.

Karlsson E.K.

Wong S.H.

et al. 1000 Genomes Project

Identifying recent adaptations in large-scale genomic data. Voight et al., 2006 Voight B.F.

Kudaravalli S.

Wen X.

Pritchard J.K. A map of recent positive selection in the human genome. ST or iHS (1% threshold) was computed from SNPs with DAF ≥ 0.2 in a 100kb window around each putatively selected locus. Significance was assessed from a beta binomial distribution fitted, in each population separately, to the observed genome-wide distribution of the proportion of outliers, to account for variations in the number of SNPs at each locus. To identify candidate eQTLs under selection, we used an outlier approach where we computed the top 1% values of Fand iHS at the genome-wide level, focusing on signals consistent with selection on derived alleles, within each population separately. To support the adaptive nature of candidate eQTLs, we computed neutral p values for each statistic using simulations based on validated demographic models of Africans and Europeans (). Furthermore, we tested for local enrichment of outliers (top 1% signals) within a 100kb-window around each eQTL (50kb on each side), similarly to previous work (). The proportion of outliers of For iHS (1% threshold) was computed from SNPs with DAF ≥ 0.2 in a 100kb window around each putatively selected locus. Significance was assessed from a beta binomial distribution fitted, in each population separately, to the observed genome-wide distribution of the proportion of outliers, to account for variations in the number of SNPs at each locus.

Archaic eQTLs and Enrichment Analyses Prüfer et al., 2014 Prüfer K.

Racimo F.

Patterson N.

Jay F.

Sankararaman S.

Sawyer S.

Heinze A.

Renaud G.

Sudmant P.H.

de Filippo C.

et al. The complete genome sequence of a Neanderthal from the Altai Mountains. Auton et al., 2015 Auton A.

Brooks L.D.

Durbin R.M.

Garrison E.P.

Kang H.M.

Korbel J.O.

Marchini J.L.

McCarthy S.

McVean G.A.

Abecasis G.R. 1000 Genomes Project Consortium

A global reference for human genetic variation. Sankararaman et al., 2014 Sankararaman S.

Mallick S.

Dannemann M.

Prüfer K.

Kelso J.

Pääbo S.

Patterson N.

Reich D. The genomic landscape of Neanderthal ancestry in present-day humans. Chang et al., 2015 Chang C.C.

Chow C.C.

Tellier L.C.

Vattikuti S.

Purcell S.M.

Lee J.J. Second-generation PLINK: rising to the challenge of larger and richer datasets. 2 > 0.8. Among these, 9,677 tagged all aSNPs in Europeans and are referred to here as “archaic tagging SNPs.” They were not necessarily aSNPs themselves, reflecting the fact that haplotypes inherited from Neandertals can harbor a mixture of different variants (i.e. variants that appeared in the Neandertal lineage, and ancient variants pre-existing in both lineages before admixture, but for which one allele is carried almost exclusively by Neandertal haplotypes in modern Europeans). We determined the level of Neandertal ancestry of the detected eQTLs, by first defining an “archaic eQTL” as an eQTL for which regulatory variants were introduced into European genomes by introgression from archaic hominins. We identified such eQTLs using the complete genome sequence of Neandertal from Altai (). Briefly, the 1000 Genomes phase 3 variants () were considered as of putative archaic origin (archaic SNPs, or aSNPs) if the Neandertal allele was present in at least one non-African individual and absent from the Yoruba population. According to this definition, 230,779 aSNPs were detected in the 100 European individuals analyzed here. We rendered the analysis more conservative, by further restricting the definition of aSNPs to those in regions of the modern human genome for which Neandertal ancestry has been predicted with a high degree of confidence (marginal probability of Neandertal ancestry ≥ 0.9 and a genetic length ≥ 0.02cM) (). This resulted in a final set of 197,959 aSNPs, of which 77,823 presented a MAF > 0.05. More than 96% of these aSNPs had an archaic allele frequency below 1% in our African samples (who are slightly differentiated from the Yoruba of 1000 Genomes phase 3), consistent with a strong enrichment in true archaic variants. To account for LD between aSNPs and characterize haplotypes that were inherited from Neandertal, we used PLINK () to extract a set of 924,362 genome-wide SNPs tagging all European variants at an r> 0.8. Among these, 9,677 tagged all aSNPs in Europeans and are referred to here as “archaic tagging SNPs.” They were not necessarily aSNPs themselves, reflecting the fact that haplotypes inherited from Neandertals can harbor a mixture of different variants (i.e. variants that appeared in the Neandertal lineage, and ancient variants pre-existing in both lineages before admixture, but for which one allele is carried almost exclusively by Neandertal haplotypes in modern Europeans). We explored the effect of introgression from Neandertals on the immune repertoire of Europeans, by counting, for each condition, the number of eQTLs overlapped by at least one archaic tagging SNP (or for which the archaic tagging SNP overlapped the reQTL in stimulated conditions), referred to here as archaic eQTLs. We then compared the number of archaic eQTLs detected with the number of SNPs expected to overlap the eQTL, when resampling SNPs tagging non-archaic haplotypes, at random from genic regions (< 1 Mb from a gene). We resampled 1,000 sets of 9,677 tagSNPs with the same allele frequency spectrum as the 9,677 archaic tagging SNPs found in Europeans, and determined their overlap with (r)eQTLs, to assess the significance of our observations. We finally report only the archaic (r)eQTLs for which at least 2 aSNPs were found to be in high LD with the eQTL peak-SNP, and (ii) the haplotype containing the largest number of archaic alleles within the eQTL was sufficiently long for the formal exclusion of incomplete lineage sorting. Huerta-Sánchez et al., 2014 Huerta-Sánchez E.

Jin X.

Asan

Bianba Z.

Peter B.M.

Vinckenbosch N.

Liang Y.

Yi X.

He M.

Somel M.

et al. Altitude adaptation in Tibetans caused by introgression of Denisovan-like DNA. Dannemann et al., 2016 Dannemann M.

Andrés A.M.

Kelso J. Introgression of Neandertal- and Denisovan-like haplotypes contributes to adaptive variation in human Toll-like receptors. The presence of aSNPs in present-day humans can be explained either by introgression or by incomplete lineage sorting (ILS). ILS occurs when an ancestral variant predating the split between humans and Neandertals is retained in both lineages, but lost from a specific human population (i.e. the African population; Figure S7 A). Given the time since introgression (47,000-65,000 years ago), the haplotypes containing alleles resulting from ILS would be expected to be shorter than those containing an aSNP introgressed from Neandertals. We distinguished between these two scenarios by first defining the core archaic haplotype for each eQTL as the haplotype within the eQTL carrying the longest stretch of archaic alleles, and then determining whether its size exceeded the expected length of haplotypes assuming an ILS model. We used the approach described by () and the most conservative parameters for the age of Altai Neandertal and Denisovan bones reported by (). We used the mean recombination rate calculated for a region composed of the core archaic haplotype in a region of 1Mb surrounding the haplotypes (500 kb on either side of the eQTL) in the 1000 Genomes CEU individuals (phase 1). p values were adjusted for multiple testing with the Benjamini-Hochberg procedure.