Evolution of immune genes is associated with the Black Death | Nature

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain
the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in
Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles
and JavaScript.
Subjects
Abstract
Infectious diseases are among the strongest selective pressures driving human evolution1,2. This includes the single greatest mortality event in recorded history, the first outbreak of the second pandemic of plague, commonly called the Black Death, which was caused by the bacterium Yersinia pestis3. This pandemic devastated Afro-Eurasia, killing up to 30–50% of the population4. To identify loci that may have been under selection during the Black Death, we characterized genetic variation around immune-related genes from 206 ancient DNA extracts, stemming from two different European populations before, during and after the Black Death. Immune loci are strongly enriched for highly differentiated sites relative to a set of non-immune loci, suggesting positive selection. We identify 245 variants that are highly differentiated within the London dataset, four of which were replicated in an independent cohort from Denmark, and represent the strongest candidates for positive selection. The selected allele for one of these variants, rs2549794, is associated with the production of a full-length (versus truncated) ERAP2 transcript, variation in cytokine response to Y. pestis and increased ability to control intracellular Y. pestis in macrophages. Finally, we show that protective variants overlap with alleles that are today associated with increased susceptibility to autoimmune diseases, providing empirical evidence for the role played by past pandemics in shaping present-day susceptibility to disease.
Main
Infectious diseases have presented one of the strongest selective pressures in the evolution of humans and other animals1,2. Not surprisingly, many candidates for population-specific positive selection in humans involve immune response genes, consistent with the hypothesis that exposure to new and/or re-emerging pathogens has driven adaptation5,6. However, it is challenging to connect signatures of natural selection with their causative pathogens unless the underlying loci are still associated with susceptibility to the same pathogen in modern populations7,8. Clarifying the dynamics that have shaped the human immune system is key to understanding how historical diseases contributed to disease susceptibility today.
We sought to identify signatures of natural selection in Europeans imposed by Yersinia pestis, the bacterium responsible for bubonic plague3. The first recorded plague pandemic began with the Plague of Justinian in ad 541 (refs. 9,10). Nearly 800 years later, the Black Death (1346–1350) marked the beginning of the second pandemic of plague, which spread throughout Europe, the Middle East and Northern Africa, reducing the population by up to 30–50%4,11. With no recent exposure to plague, Europeans living through the Black Death probably represented immunologically naive populations with little to no prior adaptation to Y. pestis. The high mortality rate suggests that genetic variants that conferred protection against Y. pestis infection might have been under strong selection during this time. Indeed, the nearly decadal plague outbreaks over the subsequent four hundred years of the second pandemic in Europe often (but not always) were associated with reduced mortality rates11,12, which could have been due to pathogen evolution or changing cultural practices, but potentially also linked to human genetic adaptation to Y. pestis.
Positive selection on immune genes
Genomic targets of selection imposed by Y. pestis during the Black Death, if present, have remained elusive13,14,15. To better identify such loci, we characterized genetic variation from ancient DNA extracts derived from individuals who died shortly before, during or soon after the Black Death in London and across Denmark. This unique sampling design differentiates, to the greatest extent possible, signatures due to Y. pestis from those associated with other infectious diseases or other selective processes (although we cannot exclude these entirely). From London, individuals were sampled from three cemeteries close to one another, tightly dated by radiocarbon, stratigraphy and historical records to before, during and after the Black Death (Fig. 1, Supplementary Table 1 andSupplementary Methods). From Denmark, individuals were sampled from five localities, geographically spread across the country, which were dated via archaeological means (such as burial arm positions), stratigraphy and historical records. We grouped all individuals into those that lived pre-Black Death (London: around ad 1000–1250, Denmark: around ad 850 to around ad 1350) and post-Black Death (London: ad 1350–1539, Denmark: around ad 1350 to around ad 1800). Within London, we also included individuals buried in the plague cemetery, East Smithfield, all of whom died during a two-year window of the Black Death between 1348 and 1349 (ref.16). Analysis of the mitogenomic diversity from these individuals identifies solely European mitogenomic haplotypes, avoiding a possible confound between natural selection and population replacement from non-European sources17.
Fig. 1: East Smithfield mass burial and sample locations, along with date ranges and final sample numbers used for the present study.
a, East Smithfield Black Death mass burial site from 1348–1349 (reproduced with permission of the Museum of London Archaeology, copyright MOLA). b, Site locations and archaeological site codes (in parentheses) for samples from London (inset, after ref. 48, Museum of London Archaeology) and from across Denmark. c, Top, population size estimates for London for around six centuries (data are from refs. 13,49,50) (Supplementary Table 1). Bottom, site locations with associated date ranges. Coloured boxes indicate date range for samples, numbers in boxes indicate samples meeting all criteria for inclusion in final analyses (see main text and Supplementary Information). Number in green star stems from East Smithfield and the dashed line refers to the time of the Black Death.
In total we screened 516 samples (n = 318 from London; n = 198 from across Denmark) for the presence of human DNA using a modified polymerase chain reaction (PCR) assay for the single copy nuclear c-myc gene17,18 and identified 360 with sufficient endogenous DNA content for downstream enrichment and sequencing of additional nuclear loci (Supplementary Methods). As many of our samples were poorly preserved and had low endogenous DNA content, we used hybridization capture to enrich for and sequence 356 immune-related genes, 496 genome-wide association study (GWAS) loci previously associated with immune disorders and 250 putatively neutral regions (1.5 kb each), as defined by Gronau and colleagues19 (Supplementary Table 2; Supplementary Methods). The targeted immune genes were manually curated on the basis of their role in immune-related processes, and include innate immune receptors, key immune transcription factors, cytokines and chemokines, and other effector molecules (Supplementary Table 3). To ensure that deamination and other forms of ancient DNA damage did not lead to spurious genotype calls, we trimmed 4 base pairs (bp) from the start and end of each sequencing read (Supplementary Fig. 1) and excluded all singleton variants (n = 106,757). Our final dataset contained 33,110 biallelic variants within the targeted regions (2,669 near GWAS loci, 19,972 in immune genes and 10,469 in putatively neutral regions), with a mean coverage of 4.6× reads per site per individual (see Supplementary Table 1 for per-individual coverage). We further filtered our results by excluding samples with missing genotype calls at more than 50% of those sites (retaining n = 206 individuals, Fig. 1c) and excluding variants with genotype calls for less than 10 individuals per time period and population. Using genotype likelihoods, we then calculated the minor allele frequency (MAF) per population at each time point. Finally, we retained only sites with a mean MAF (averaged across London and Denmark) greater than 5% (n = 22,868 sites), as our power to detect selection for variants below 5% is very low (Supplementary Fig. 2).
To detect alleles that may have conferred protection from, or susceptibility to, Y. pestis, we searched within candidate regions (immune genes and GWAS loci) for variants that exhibit unexpectedly large changes in allele frequency between pre- and post-Black Death samples. Specifically, we identified alleles for which the degree of differentiation (FST) was larger than expected by chance, when compared to variants in the putatively neutral genetic regions sampled from the same populations. We used the larger sample set from London as our discovery cohort. Burials in London were also more precisely dated and better geographically controlled than those from Denmark, improving our relative ability to detect selection in the cohort from London (Supplementary Methods). We found an enrichment of highly differentiated variants for all frequency bins with MAF greater than 10%, relative to a null expectation established using our neutral loci (Fig. 2a). Across these variants, differentiation at immune loci exceeded the 99th percentile of neutral variants at 2.4× the rate expected by chance (binomial test P = 7.89 × 10−12). For variants with an MAF greater than 30%, this enrichment was even more pronounced (3.9× the rate expected by chance; binomial test P = 1.16 × 10−14), probably due to increased power (Supplementary Fig. 2). Simulations show that differences in recombination rate and background selection between neutral and candidate loci are insufficient to explain the observed enrichments (Extended Data Fig. 1). To further validate the signatures of selection we observed among immune loci from our London sample, we performed the same analyses using the allele frequencies estimated from our Danish cohort. These samples were also enriched for highly differentiated sites relative to the expectation from neutral loci (1.6× the rate expected by chance, binomial test P = 9.21 × 10−4; Fig. 2b), further supporting evidence for plague-induced selection on immune genes.
Fig. 2: Positive selection at immune loci.
a,b, Enrichment of highly differentiated sites in functional regions relative to neutral regions when comparing the pre-Black Death (BD) population to the post-BD population in London (a) and Denmark (b). c, FST between London before and after the Black Death, limited to the 535 sites that show qualitative patterns consistent with natural selection (namely allele frequency changes in the same direction in both London and Denmark after the Black Death, and the opposite direction for individuals who died during the Black Death (Supplementary Table 4)). Manhattan plot showing loci with patterns indicative of positive selection. Point size and colour intensity (which alternates by chromosome) represents the –log10P value comparing populations in London before and after the plague, points coloured in orange represent the four positions and their associated genes, which are highly differentiated in Denmark as well. d–g, Patterns of allelic change over time for the four strongest candidates for positive selection. Error bars represent the standard deviation based on bootstrapping individuals from that population and each time point 10,000 times. Allele frequencies for London are shown in red and for Denmark are shown in blue. Modern allele frequencies are derived from 1000 Genomes data for Great Britain in London51.
To identify specific loci that represent the strongest candidates of selection, we applied a series of stringent criteria that leveraged the time periods and populations in this unique dataset. First, we identified 245 common variants (MAF greater than 10%) that were highly differentiated (FST greater than 95th percentile defined using neutral sites) when comparing pre- versus post-Black Death samples in London alone (Supplementary Table 4). Next, we reasoned that variants conferring increased susceptibility to, or protection from, Y. pestis should show opposing frequency patterns before, during and after the Black Death. Specifically, variants associated with susceptibility should increase in frequency among people who died during the Black Death and should decrease in frequency among individuals sampled post-Black Death (that is, the survivors and/or descendants of the survivors). Conversely, protective variants should show an inverse pattern. Using this reasoning, we narrowed down our list of putatively selected loci from 245 to 35 (Supplementary Table 4). Finally, we asked if these loci were also highly differentiated before and after the Black Death in our Danish replication cohort (that is, among the top 10% most highly differentiated sites, and in the same direction as seen in London). Four loci met these criteria, representing the strongest candidates for selection (Fig. 2c–g). We calculated the selection coefficient (s) for each of these variants using a Hidden Markov Model (HMM) framework (based on ref. 20, Supplementary Methods). Statistical support for non-neutral evolution (s ≠ 0) among our four candidate loci is strong when compared to that of neutral loci (P < 0.001 for each locus; Supplementary Table 5). Despite the large confidence intervals—an inherent limitation when trying to estimate s over a few generations—the absolute values for the point estimate of s range from 0.26 to 0.4, which are among the strongest selective coefficients yet reported in humans, to our knowledge (Extended Data Fig. 2 and Supplementary Table 5).
Functional dissection of candidate loci
None of our top candidate variants overlaps with (nor is in strong linkage disequilibrium with) coding variants, although one, near the ERAP2 gene, is strongly linked to a variant that affects splicing21,22. Their selective advantage may stem from an impact on gene expression levels, particularly in immune cell types that participate in the host response to Y. pestis infection. Macrophages in particular are recruited to sites of infection, where they interact with bacteria and contribute to plague resistance23,24. Macrophages phagocytize Y. pestis, but some bacteria survive and spread to the lymph node, where they replicate uncontrollably25,26. To test whether the four candidate loci we identified, or genes near them, are involved in the transcriptional response to Y. pestis, we incubated monocyte-derived macrophages from 33 individuals with heat-killed Y. pestis and compared their expression profiles to unstimulated control samples using RNA sequencing (Supplementary Methods). As expected, macrophages responded robustly to Y. pestis, such that principal component 1 of the gene expression data, which separates baseline versus Y. pestis stimulated conditions, explains 56% of the total variance (Extended Data Fig. 3).
Seven genes within 100 kb of our four candidate loci were expressed in this dataset: locus 1 (rs2549794): ERAP1, ERAP2, LNPEP; locus 3 (rs11571319): CTLA4, ICOS; and locus 4 (rs17473484): TICAM2, TMED7. NFATC1, the only gene within 100 kb of locus 2 (rs1052025), was not expressed in this dataset. With the sole exception of LNPEP, all of these genes were differentially expressed in response to Y. pestis stimulation (Fig. 3a), supporting their putative role in the host response. Macrophages from an additional panel of eight individuals infected with live and fully virulent Y. pestis showed similar directional changes in gene expression to those observed in response to heat-killed bacteria. This was true genome-wide (r = 0.88, P < 1 × 10−10, Extended Data Fig. 4a) and for genes near our four candidate loci (ERAP1 is an exception; it was upregulated in response to live bacteria but downregulated in response to heat-killed bacteria, Extended Data Fig. 4b). To investigate whether changes in gene expression were specific to Y. pestis or shared with other infectious agents, we analysed gene expression data from macrophages infected with live Listeria monocytogenes (a Gram-positive bacterium) and Salmonella typhimurium (a Gram-negative bacterium, like Y. pestis)27, as well as monocytes activated with bacterial and viral ligands targeting the Toll-like receptor (TLR) pathways (TLR1/2, TLR4 and TLR7/8) and live influenza virus28. These data show that all genes near our candidate loci (with the exception of CTLA4) respond to other pathogenic agents but that the direction of change in expression differs depending on the stimulus. For example, ERAP2 is downregulated in response to all live bacteria or bacterial stimuli, including Y. pestis, but is upregulated in response to viral ones (Extended Data Fig. 5).
a, Normalized log2 fold-change of genes within 100 kb of candidate variants in response to incubation of primary macrophages for four hours with heat-killed Y. pestis. Dark grey dots correspond to the fold-change observed for each of the 33 individuals tested; red dots and bars represent the mean ± standard deviation. With the exception of LNPEP, all associations are significant (10% false discovery rate). b,c, Effect of rs1747384 genotype upon TICAM2 expression (b) and rs2549794 genotype upon ERAP2 expression (c), split by macrophages stimulated for four hours with heat-killed Y. pestis and unstimulated macrophages. Red dots and bars represent the mean ± standard deviation. d, Comparison of ERAP2 expression among non-infected and infected cells with live Y. pestis for five hours, profiled using scRNA-sequencing data in individuals with homozygous rs2549794 genotypes. Colour intensity reflects the level of ERAP2 expression, standardized for unstimulated or infected cells. Major PBMC cell types are labelled and can be found clearly coloured in Extended Data Fig. 7; CD4+ and CD8+ T cells were analysed separately. NK, natural killer. e, Effects of rs2549794 genotype upon ERAP2 expression, split by infected and non-infected conditions, for each cell type. f, Illustration of the two haplotype-specific ERAP2 spliced forms for Haplotype A and Haplotype B with start (green) and stop (brown) codons. Below we show the effect of rs2248374 genotype upon total mRNA expression of ERAP2 (left) and the specific expression of the isoform (Haplotype B) encoding the truncated version of ERAP2 (right). For all panels: ***P < 0.001; **P < 0.01; *P < 0.05.
Having established that genes near our candidate loci show a transcriptional response to Y. pestis in macrophages, we next asked whether genetic variation at each of the four candidate loci is associated with gene expression levels at nearby genes (Fig. 3b,c). We identified an association between rs17473484 genotype and TICAM2 expression in which the protective allele is associated with higher expression of the gene in the unstimulated condition (Fig. 3b; P = 2.5 × 10−6) but not after Y. pestis stimulation (P = 0.24). This effect is intriguing as TICAM2 encodes an adaptor protein for TLR4. In vivo, TLR4 detects Y. pestis via the recognition of lipopolysaccharides (LPS) on the bacterial outer membrane29. Y. pestis attempts to circumvent this detection by deacylating surface LPS, thereby reducing the binding affinity for TLR4 (ref. 30). TICAM2 ushers LPS-bound TLR4 into endosomes and activates type I interferon responses31. It is therefore possible that increased TICAM2 expression confers protection against Y. pestis by increasing sensitivity to LPS and promoting an effective immune response.
The strongest association we identified was between rs2549794 and ERAP2 expression, in which the protective allele (C) is associated with a fivefold increase in expression of ERAP2 relative to the putatively deleterious T allele (Fig. 3c), in both unstimulated (P = 4.4 × 10−10) and Y. pestis-challenged (P = 8.7 × 10−7) cells. We observed similarly strong associations in macrophages and monocytes infected with other pathogens (Salmonella, Listeria, influenza) or stimulated with TLR-activating ligands (all P < 1 × 10−10; Extended Data Fig. 6). This pattern also generalizes to Y. pestis-infected peripheral blood mononuclear cells (PBMCs) more broadly, suggesting that rs2549794 is associated with ERAP2 expression levels regardless of the infectious/inflammatory stimulus or cell type. In detail, we generated single-cell RNA-sequencing data from PBMCs from ten individuals (five homozygous for the selectively favoured rs2549794 C allele and five homozygous for the alternate T allele), both at baseline and after infection with live, fully virulent Y. pestis. Across all immune cell types profiled—B cells, CD4+ T cells, CD8+ T cells, natural killer cells and monocytes (Extended Data Fig. 7)—we identified 5,570 genes for which infection with Y. pestis significantly alters gene expression levels (314 to 4,234 genes per cell type, 10% false discovery rate; Supplementary Table 6). Most genes near our candidate loci are differentially expressed in response to Y. pestis infection, but both the magnitude and direction of such effects is cell type-specific (Extended Data Fig. 8). For example, ERAP2 is upregulated upon stimulation in all PBMC cell types (Fig. 3d,e), but is downregulated in monocyte-derived macrophages. These differences could be due to differences in the transcription factors and enhancers active in each cell type, differences in our infection models (PBMCs versus monocyte-derived macrophages), or both. Notably, in all cell types and under all conditions, the selectively favoured rs2549794 C allele is associated with increased ERAP2 expression compared to the alternative T allele.
The ERAP2 locus is characterized by two haplotypes (A and B) that are common around the world21. Haplotype A encodes the canonical (full-length) ERAP2 protein consisting of 960 amino acids; Haplotype B is characterized by the presence of the G allele of rs2248374, a splice-site-altering variant that leads to the production of a splice isoform with an elongated exon 10 containing two premature stop codons. This ERAP2 isoform undergoes nonsense-mediated decay (NMD), resulting in undetectable ERAP2 protein levels21. Even when produced, the shorter protein encoded by Haplotype B has limited aminopeptidase activity32. The selectively favoured rs2549794 C allele is in strong linkage with allele A of rs2248374 (R2 > 0.8, D' = 1.0) of Haplotype A, which encodes the full-length ERAP2 protein, whereas the deleterious rs2549794 T allele is in linkage with the rs2248374 G allele of Haplotype B, which encodes the truncated version of ERAP2. Using real-time PCR specific to the mRNA encoded by Haplotype B, we found that, in macrophages, the decreased ERAP2 expression in individuals harbouring the deleterious rs2248374 G allele is coupled with a higher expression of the truncated isoform known to undergo NMD (Fig. 3f, P = 2.66 × 10−6). PCR amplification of exon 10 further confirmed that individuals homozygous for the protective rs2248374 A allele only express Haplotype A, whereas individuals homozygous for the rs2248374 G allele almost exclusively express the truncated Haplotype B (Extended Data Fig. 9).
ERAP1 and ERAP2 are aminopeptidases that work synergistically to trim peptides for presentation to CD8+ T cells by major histocompatibility complex (MHC) class I molecules33. Given their central role in antigen presentation, it is not surprising that polymorphisms near these genes, including rs2549794, have been associated with susceptibility to a variety of infectious agents34,35. ERAP2 deficiency leads to a significant remodelling of the repertoire of antigens that are presented by MHC to CD8+ T cells36,37, including MHC ligands that have high homology to peptide sequences derived from Yersinia species37. Having ERAP2-mediated aminopeptidase activity plausibly helps to promote the presentation of a more diverse array of Yersinia-derived antigens to CD8+ T cells, which in turn has an important role in protection against infection38,39. Indeed, mice depleted of CD8+ T cells all die within one-week post-infection with a milder Yersinia spp., Y. pseudotuberculosis, whereas all wild-type mice survive39. Unfortunately, rodent models only possess a single ERAP aminopeptidase that is homologous to ERAP1, limiting our ability to directly test the function of ERAP2 on antigen presentation and response to Y. pestis in vivo.
In addition to its canonical role in antigen presentation and CD8+ T cell activation, ERAP2 is also involved in viral clearance and cytokine responses40. We therefore sought to test whether ERAP2 genotype is associated with variation in the cytokine response to Y. pestis infection. To do so, we infected a further set of monocyte-derived macrophages from 25 individuals (9 homozygous for the selectively favoured ERAP2 haplotype, 9 heterozygous and 7 homozygous for the deleterious haplotype) with live, virulent Y. pestis and measured the protein levels of ten cytokines involved in various aspects of the immune response, at baseline and at 24 h postinfection. No differences in cytokine levels existed at baseline (all P > 0.05; Supplementary Table 7), but four cytokines showed a significant association with ERAP2 genotype upon stimulation. Specifically, the levels of granulocyte colony-stimulating factor (P = 0.0155), interleukin (IL)-1β (P = 0.00262) and IL-10 (P = 0.00248) significantly decreased with the number of protective C alleles, whereas we observed the opposite pattern for levels of the chemokine CCL3 (P = 0.0216), which is involved in the recruitment of neutrophils upon infection (Fig. 4a–d; Supplementary Table 7). Although the exact mechanism linking ERAP2 function to variation in cytokine responses in response to Y. pestis remains to be determined, we speculate that it might be due to the role played by ERAP2 in the differentiation and/or activation of myeloid cells and in the activation of the caspase1/NLRP3 inflammasome pathway40. Finally, we assayed the ability of macrophages to control internalized Y. pestis replication as a function of ERAP2 genotype. We observed an additive association between ERAP2 genotype and the ability to limit Y. pestis infection, such that individuals with more copies of the selectively favoured allele were better able to restrict intracellular replication (Spearman’s rho = 0.42, P = 0.0155; Fig. 4e). This experiment directly implicates ERAP2 in the response to Y. pestis infection, but through pathways other than ERAP2’s canonical role in antigen presentation. Together, these results support the idea that changes in ERAP2 allele frequencies during the Black Death were probably due to Y. pestis-induced natural selection.
a–d, Effect of genotype upon cytokine levels for granulocyte colony-stimulating factor (G-CSF) (a), interleukin-1β (IL-1β) (b), interleukin 10 (IL-10) (c) and C-C motif chemokine ligand 3 (CCL3) (d). Remaining cytokines showed no significant effects and are included in Supplementary Table 7. e, Boxplots showing the percentage of bacteria killed (y axis) by macrophages infected for 24 h as a function of ERAP2 genotype (x axis). The percentage of bacteria killed was calculated as the CFU2 h − CFU24 h/CFU2 h, where CFU is colony-forming unit. The P value results from a linear model examining the association between ERAP2 single nucleotide polymorphism genotypes (SNP) (coded as the number of protective rs2549794 alleles found in each individual: 0, 1 or 2) and the percentage of bacteria killed.
Discussion
Our results provide strong empirical evidence that the Black Death was an important selective force that shaped genetic diversity around some immune loci. Although our candidate selected loci are generally involved in the response to pathogens, and not just Y. pestis, our unique sampling design minimized the degree to which other historical events (such as tuberculosis8 or famine41,42,43) could have affected the inference of selection. To support our ancient genomic data, we confirmed that the strongest candidates for positive selection are directly involved in the immune response to Y. pestis using functional data from in vitro infection experiments.
We identified four loci that were strongly differentiated before and after the Black Death in London and replicated in our Danish cohort as the strongest candidates of selection. However, given our small sample sizes and the low sequencing depth from some samples, our replication power is limited. Thus, some of the other 245 highly differentiated loci in London were also probably impacted by natural selection, although they did not survive our conservative filtering criteria. Increased sample sizes coupled with additional functional data will help to further dissect the evolutionary role played by these variants in the immune response to Y. pestis.
ERAP2 showed the most compelling evidence for selection, both from a genetic and functional perspective, with an estimated selection coefficient of 0.4 (95% confidence interval 0.19,0.62, Extended Data Figs. 2 and 10). This estimate suggests that individuals homozygous for the protective allele were about 40% more likely to survive the Black Death than those homozygous for the deleterious variant. This allele is associated with both increased expression of ERAP2 and production of the canonical full-length ERAP2 protein21,22. We suggest that this protein increases the presentation of Yersinia-derived antigens to CD8+ T cells, stimulating a protective immune response against Y. pestis38,39. Furthermore, we show that macrophages from individuals possessing the selected ERAP2 allele engage in a unique cytokine response to Y. pestis infection and are better able to limit Y. pestis replication in vitro.
In general, individuals with more copies of the selectively advantageous haplotype displayed a weaker cytokine response to infection but a better ability to limit bacterial growth. For example, levels of IL-1β, a key proinflammatory cytokine often associated with pyroptotic cell death44, were threefold lower in individuals homozygous for the advantageous ERAP2 genotype when compared to individuals homozygous for the putatively deleterious one. Therefore, subjects with the advantageous haplotype are both more efficient at controlling internalized bacteria and at resisting Y. pestis-induced cell death than subjects with the deleterious haplotype—abilities that may help to reduce bystander tissue damage during infection. However, as our experiments were done in vitro we were unable to directly evaluate the impact of ERAP2 genotype on tissue damage, immune cell recruitment and survival.
More broadly, our results highlight the contribution of natural selection to present-day susceptibility towards chronic inflammatory and autoimmune disease. We show that ERAP2 is transcriptionally responsive to stimulation with a large array of pathogens, supporting its key role in the regulation of immune responses. Therefore, selection imposed by Y. pestis on ERAP2 probably affects the immune response to other pathogens or disease traits. Consistent with this hypothesis, the selectively advantageous ERAP2 variant is also a known risk factor for Crohn’s disease45, and ERAP2 variation has also been associated with other infectious diseases34,35. Thus, selection for pathogen defence in the presence of pathogens such as Y. pestis may be counterbalanced against the costs of immune disorders, resulting in a long-term signature of balancing selection21,22. Likewise, another of our top candidate loci (rs11571319 near CTLA4) is associated with an increased risk of rheumatoid arthritis46 and systemic lupus erythematosus47, such that retaining the putatively advantageous allele during the Black Death confers increased risk for autoimmune disease in present-day populations. To date, most of the evidence for an association between autoimmune risk alleles and adaptation to past infectious diseases remains indirect, primarily because the aetiological agents driving selection remain hidden. Our ancient genomic and functional analyses suggest that Y. pestis has been one such agent, representing empirical evidence connecting the selective force of past pandemics to present-day susceptibility to disease.
Reporting summary
Data availability
Hybridization capture data from the ancient individuals have been deposited in the NCBI Sequence Read Archive (SRA) under BioProject PRJNA798381. Expression data have been deposited into the NCBI Gene Expression Omnibus (GEO) under project GSE194118 (for macrophages) and the NCBI SRA under project accession PRJNA871128 (for PBMCs). Cytokine data is available in Supplementary Table 8 and CFU data in Supplementary Table 9.
Acknowledgements
We thank all members of the Barreiro laboratory and the Poinar laboratory for their constructive comments and feedback. We thank J. Tung for her comments and edits to the manuscript. Computational resources were provided by the University of Chicago Research Computing Center. Sequencing was performed at the Farncombe Sequencing Facility McMaster University. We thank the Cytometry and Biomarkers platform at the Institut Pasteur for support in conducting this study, with a special thanks to C. Petitdemange for help running the Luminex assay. We thank X. Zhang for assistance in simulating allele frequency changes under neutral evolution. This work was supported by grant R01-GM134376 to L.B.B., H.P. and J.P.-C., a grant from the Wenner-Gren Foundation to J.F.B. (8702), and the UChicago DDRCC, Center for Interdisciplinary Study of Inflammatory Intestinal Disorders (C-IID) (NIDDK P30 DK042086). The SSHRC Insight Development Grant supported the collection of the Danish samples (430-2017-01193). H.N.P. was supported by an Insight Grant no. 20008499 from the Social Sciences and Humanities Research Council of Canada (SSHRC) and The Canadian Institute for Advanced Research under the Humans and the Microbiome programme. T.P.V. was supported by NIH F32GM140568. X.C. and M. Steinrücken were supported by grant R01GM146051. We also thank the University of Chicago Genomics Facility (RRID:SCR_019196), especially P. Faber, for their assistance with RNA sequencing. H.P. thanks D. Poinar for continued support and manuscript suggestions and editing.
Contributions
L.B.B. and H.N.P. directed the study. J.K. designed the enrichment assays and generated ancient genomic data. T.P.V. led all data and computational analyses, with contributions from J.-C.G. X.C. performed all analyses to estimate selection coefficients under the supervision of M. Steinrücken. J.F.B., R.S. and V.Y. performed challenge experiments with macrophages and heat-killed Y. pestis. M. Shiratori and A. Dumaine performed the infection experiments on PBMCs and generated the single-cell RNA-sequencing data, with assistance from D.E. and D.M. C.E.D. and J.P.-C. performed and designed the infection experiments with live Y. pestis on macrophages and generated both cytokine and CFU data, with assistance from J.M. and R.B. C.J.Y. and M.B. designed the probes to quantify the isoform encoding the short version of ERAP2 and M.I.P. performed the experiments. R.R., S.N.D., J.A.G. and J.L.B. provided access to samples, archaeological information, including dating, and other relevant information. A.C. and N.V. provided insights into historical context. K.E. and G.B.G. provided additional sampling and bioinformatic processing and cluster maintenance. A. Devault and J.-M.R. provided insight on targeted enrichment and modified versions of baits used for immune enrichment. G.A.R. provided genomic input on loci and contributed financially to the sequencing of targets. T.P.V., J.K., H.N.P. and L.B.B. wrote the manuscript, with input from all authors.
Extended data figures and tables
Comparison of recombination rate (A) and background selection levels (B) between neutral loci and our candidate regions. Candidate regions were stratified into those which were tested and those which were candidates for positive selection based on high differentiation in London pre- vs post-BD. (C) Forward simulations matched for the rates of recombination and background selection of the regions targeted in our study show a slight enrichment of highly differentiated sites in candidate regions, but far from the level of enrichment observed in our collected data (D), replicated from Fig. 2a for comparison. For example, whereas our data differentiation at immune loci exceeded the 99th percentile of neutral variants at 2.4x the rate expected by chance (among variants with a MAF > 10%), the same enrichment is less than 1.2x in the simulated data.
(A) Distributions of \({\hat{s}}_{{\rm{MLE}}}\) for the four SNPs when replicates are simulated with the corresponding bootstrapped allele frequency distributions as initial conditions and bootstrap-corrected estimates \({\widetilde{s}}_{{\rm{MLE}}}\). Whiskers on the violin plots label the 2.5-, 50-, and 97.5-percentiles of their respective distributions. (B) ROC and (C) Precision-Recall curves for the estimation procedure to distinguish replicates under selection from those under neutrality.
Effect size of Y. pestis stimulation compared between heat-killed Y. pestis (x-axis, n = 33 individuals) and live Y. pestis (y-axis, n = 8 individuals). (A) shows all genes, with a blue line representing the best fit line (r = 0.88). (B) compares effect sizes at genes near candidates for positive selection profiled in both expression datasets (red: heat-killed; purple: live bacteria). Error bars represent the standard error in estimating the effect size. The direction of effect is consistent except for LNPEP (which is not significant in either analysis) and ERAP1. CTLA4 and TICAM2 are not shown because they were not expressed at sufficiently high levels in the macrophages from the 8 individuals infected with live Y. pestis. Asterisks placed near the point estimate of each value represent the significance: *** p < 0.001; ** p < 0.01; * p < 0.05.
Data are derived from Nedelec et al27. and Quach et al28 Nedelec et al. measured the gene expression response of monocyte-derived macrophages to infection with two live intracellular bacteria: Listeria monocytogenes (a Gram-positive bacterium) and Salmonella typhimurium (a Gram-negative bacterium). Quach et al. characterized the transcriptional response at 6 h of primary monocytes to bacterial and viral stimuli ligands activating Toll-like receptor pathways (TLR1/2, TLR4, and TLR7/8) and live influenza virus. The data for Y. pestis are the fold change responses observed in response to heat-killed bacteria. A negative estimate in plot (purple) indicates that the gene is downregulated and a positive value (red) indicates that the gene is up-regulated. The statistical support for the reported changes is given by the associated p values. Larger circle sizes represent smaller p values and empty circles refer to cases where the gene was not expressed in that dataset.
Effect of genotype at nearby loci on the expression of ERAP2 (top) and TICAM2 (bottom), across experimental conditions in this study and previously published26,27. For ERAP2, the protective “C” haplotype increases expression in all conditions (p < 0.001). For TICAM2, the protective reference haplotype decreases expression only in the unstimulated condition (p = 2.5x10−6; p > 0.05 in all other conditions).
UMAP projection of single-cell RNA sequencing data of non-infected cells and cells infected with live Y. pestis for five hours, after integrating samples. Major immune cell types cluster separately and cells are colored by the cell type to which they were assigned.
For each cell type profiled using single-cell RNA sequencing, we show the effect of Y. pestis infection upon gene expression. A negative estimate (purple) indicates that the gene is downregulated and a positive value (red) indicates that the gene is up-regulated. The statistical support for the reported changes is given by the associated p values. Larger circle sizes represent smaller p values and empty circles refer to cases where the gene was not expressed in that cell type.
Bioanalyzer traces showing the results of the PCR amplification of cDNA across the exon 10 splice junction from macrophages of 10 individuals with different genotypes for the slice variant rs2248374. The genotype of each individual is shown on top. A negative PCR control was also performed using water. The G allele at rs2248374 is predicted to produce an elongated exon 10 containing two premature stop codons (red rectangles), leading to nonsense mediated decay.
The time axis serves as an approximate reference of the relative sampling times for the empirical samples. Dashed vertical lines indicate the relative sampling time for each group of samples considered in the analysis, and the floating boxes with orange and blue dots represent pools of samples from a bi-allelic locus. Above and below the time axis are sketches that respectively correspond to the simulation scheme and the likelihood computations. The shaded red horizontal tree represents the population-continuity along approximate time (x-axis), with the Black Death pandemic occurring in the dark shaded period. The shortened branch with a skull at the end represents people who died of the disease. In each simulated replicate, ∆p and −∆p mark the respective changes of allele frequency during the pandemic in the mid-pandemic and post-pandemic sample pools. In the inference schematics, each horizontal straight line represents a sampling scheme from which a likelihood was computed. Lightning bolts labeled with s or −s represent the selection coefficients.
Rights and permissions
Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.