Природа (2023)Цитировать эту статью
12 тысяч доступов
293 Альтметрика
Подробности о метриках
Люди демонстрируют значительную межиндивидуальную клиническую изменчивость после заражения SARS-CoV-21,2,3, генетическая и иммунологическая основа которой начала расшифровываться4. Однако степень и причины популяционных различий в иммунных реакциях на SARS-CoV-2 остаются неясными. Здесь мы сообщаем данные одноклеточного РНК-секвенирования мононуклеарных клеток периферической крови — от 222 здоровых доноров различного происхождения, — которые были стимулированы SARS-CoV-2 или вирусом гриппа А. Мы показываем, что SARS-CoV-2 индуцирует более слабую, но более гетерогенную, стимулируемую интерфероном генную активность по сравнению с вирусом гриппа А, а также уникальную провоспалительную характеристику в миелоидных клетках. Транскрипционные ответы на вирусы демонстрируют заметные популяционные различия, в первую очередь обусловленные изменениями численности клеток, включая усиление лимфоидной дифференцировки, связанное со латентной цитомегаловирусной инфекцией. Локусы количественных признаков экспрессии и анализ медиации показывают широкое влияние клеточного состава на популяционные различия в иммунных ответах, при этом генетические варианты оказывают сильное влияние на определенные локусы. Кроме того, мы показываем, что естественный отбор увеличил популяционные различия в иммунных реакциях, особенно для вариантов, связанных с реакцией на SARS-CoV-2 у жителей Восточной Азии, и документируем клеточные и молекулярные механизмы, с помощью которых интрогрессия неандертальцев изменила иммунные функции, такие как реакция миелоидных клеток к вирусам. Наконец, анализ колокализации и общетранскриптомных ассоциаций выявляет перекрытие между генетической основой иммунных ответов на SARS-CoV-2 и тяжесть COVID-19, что дает представление о факторах, способствующих текущим различиям в риске заражения COVID-19.
Примечательной особенностью пандемии COVID-19 является значительная клиническая вариабельность среди лиц, инфицированных SARS-CoV-2, варьирующаяся от бессимптомной инфекции до смертельного заболевания1,2,3. Факторы риска включают пожилой возраст1, а также мужской пол5, сопутствующие заболевания6 и генетику хозяина4,7,8. Кроме того, вариации врожденного иммунитета9,10,11, включая врожденные ошибки или нейтрализующие аутоантитела против интерферонов I типа12,13,14, способствуют вариациям клинических результатов, а эпидемиологические и генетические данные свидетельствуют о различиях между популяциями6,7,15,16 . Это, наряду с сообщениями о связанных с происхождением различиях в транскрипционных реакциях на иммунные вызовы17,18,19, требует исследования масштабов и факторов вариаций иммунных ответов на SARS-CoV-2 среди населения во всем мире.
Давление отбора, вызванное патогенами, имело первостепенное значение в ходе эволюции человека20. Адаптация человека к РНК-вирусам посредством избирательного охвата или архаичных примесей была идентифицирована как источник генетической дифференциации популяций18,21,22, а сигналы адаптации были зарегистрированы во взаимодействующих с коронавирусом белках у жителей Восточной Азии23,24. Существуют также доказательства связи между архаической интрогрессией и иммунитетом25, при этом неандертальские гаплотипы связаны с тяжестью заболевания COVID-1926,27. Однако влияние естественного отбора и архаичной примеси на иммунные реакции на SARS-CoV-2 еще предстоит изучить.
Мы ответили на эти вопросы, подвергнув мононуклеарные клетки периферической крови (РВМС) лиц центральноафриканского, западноевропейского и восточноазиатского происхождения воздействию SARS-CoV-2 и, для сравнения, вирусу гриппа А (IAV). Объединив секвенирование одноклеточной РНК (scRNA-seq) с подходами количественной и популяционной генетики, мы определяем экологические и генетические факторы популяционных различий в иммунных ответах на SARS-CoV-2.
Мы охарактеризовали транскрипционные ответы на SARS-CoV-2 и IAV, выполнив scRNA-seq анализ РВМС от 222 доноров, ранее не болевших SARS-CoV-2, происходящих из трех географических регионов (Центральная Африка, n = 80 мужчин; Западная Европа, n = 80). мужчины; Восточная Азия, n = 36 женщин и 26 мужчин) и с разным генетическим происхождением (дополнительный рисунок 1 и дополнительная таблица 1). РВМС обрабатывали в течение 6 часов (дополнительное примечание 1, дополнительный рисунок 2 и дополнительная таблица 2) с помощью ложного контроля (нестимулированного), SARS-CoV-2 (предковый штамм, BetaCoV/Франция/GE1973/2020) или IAV. (H1N1/PR/8/1934). Мы захватили более 1 миллиона высококачественных одноклеточных транскриптомов (рис. 1a, дополнительный рисунок 3 и дополнительная таблица 3a). Объединив кластеры на основе транскриптома с клеточной индексацией транскриптомов и эпитопов путем секвенирования (CITE-seq; Методы), мы определили 22 типа клеток в миелоидных, B, CD4+ T, CD8+ T и естественных киллерах (NK) иммунных линиях (рис. 1b). , дополнительный рисунок 4 и дополнительная таблица 3b – d). После воздействия вируса большинство типов клеток показали умеренные изменения в численности, при этом самые сильные изменения наблюдались в миелоидной линии после лечения IAV (дополнительное примечание 2 и дополнительная таблица 3e).
0.5; out of 12,655 with detectable expression; Supplementary Table 3f). These responses were highly correlated across lineages and featured a strong induction of interferon-stimulated genes (ISGs) (Extended Data Fig. 1a). However, myeloid responses were markedly heterogeneous, with SARS-CoV-2 inducing a transcriptional network enriched in inflammatory-response genes (Gene Ontology (GO): 0006954; fold-enrichment (FE) = 3.4, FDR < 4.9 × 10−8; Supplementary Table 3g). For example, IL1A, IL1B and CXCL8 were highly and specifically upregulated in response to SARS-CoV-2 (log2[FC] > 2.8, FDR < 2.3 × 10−36), consistent with in vitro and in vivo studies28,29./p> 0.2) and 898 and 652 displaying differential responses between populations (popDRGs; FDR < 0.01, |log2[FC]| > 0.2) after stimulation with SARS-CoV-2 and IAV, respectively (Fig. 2b and Supplementary Table 4b,c). popDRGs included key immunity regulators, such as the IFN-responsive GBP7 and the gene coding for the macrophage inflammatory protein MIP-3, CCL23, both of which were more strongly upregulated in Europeans (Fig. 2c). The GBP7 response was common to both viruses and all lineages (log2[FC] > 0.88, Student’s t-test, adjusted P (Padj) < 1.4 × 10−3), but that of CCL23 was specific to SARS-CoV-2-stimulated myeloid cells (log2[FC] = 0.72, Student’s t-test, Padj = 5.3 × 10−4). We estimated that population differences in cellular composition accounted for 15–47% of popDEGs and for 7–46% of popDRGs, with the strongest impact on NK cells (Fig. 2b,d and Extended Data Fig. 3c). Variation in cellular composition mediated pathway-level differences in response to viral stimulation between populations (Supplementary Table 4d). For example, in virus-stimulated NK cells, genes involved in the promotion of cell migration, such as CSF1 or CXCL10, were more strongly induced in Europeans (normalized enrichment score > 1.5, gene set enrichment analysis, Padj < 0.009). However, the loss of this signal after adjustment for cellular composition (Fig. 2e) indicates that fine-scale cellular heterogeneity drives population differences in immune responses to SARS-CoV-2./p> 3.2, Fisher’s exact test, P < 9.4 × 10−4), with more than 98% of replicated eQTLs affecting expression in the same direction (Supplementary Note 8, Supplementary Fig. 8 and Supplementary Table 5e). The correlation of reQTL effect sizes across ontogenetically related cell types was weaker than for eQTLs (r = 0.36 and 0.50, respectively, Wilcoxon’s rank-sum test, P < 5.6 × 10−13; Extended Data Fig. 4d). Furthermore, the proportion of virus-dependent reQTLs differed across cell types. In lymphoid cells, only 7.7% of reQTLs differed in effect size between viruses (interaction P < 0.01; Fig. 3b,c), whereas 49% of myeloid reQTLs were virus dependent (interaction P < 0.01), with 46 and 185 reQTLs displaying specific, stronger effects after SARS-CoV-2 and IAV stimulation, respectively. The strongest SARS-CoV-2 reQTL (rs534191, Student’s t-test, P = 1.96 × 10−16 (SARS-CoV-2) and P = 0.05 (IAV); Fig. 3d) was identified in myeloid cells at MMP1, encoding a biomarker of COVID-19 severity38. These analyses reveal that the effects of virus-induced reQTLs are cell-type dependent and highlight the virus specificity of the genetic basis of the myeloid response./p> 1.7 across lineages after stimulation). Thus, population variation in immune responses is driven largely by cellular heterogeneity, but genetic variants with marked allele frequency variation contribute to population differences at specific loci./p> 3) targeting the same (IFITM2, IFIT5) or different (ISG20, IFITM3, TRIM14) eQTLs at highly differentiated genes, suggesting repeated adaptations targeting IFN-mediated antiviral immunity (Supplementary Note 10, Supplementary Table 7c and Supplementary Fig. 10). We determined whether selection had altered gene expression in specific cell types or in response to SARS-CoV-2 or IAV by testing for increased population differentiation (PBS) at (r)eQTLs within each cell type, relative to random single-nucleotide polymorphisms (SNPs) matched for allele frequency, linkage disequilibrium (LD) and distance to the nearest gene. In the basal state, eQTLs were more strongly differentiated in Europeans, the strongest signal observed for γδ T cells (Extended Data Fig. 6a). Among popDEGs for which genetics mediates more than 50% of the differences between Africans and Europeans, 34% presented signals of rapid adaptation in Europeans (versus 21% in Africans, Fisher’s exact test, P = 7.7 × 10−6). For example, population differences at GBP7 have been driven by a frequency increase, over the last 782–1,272 generations, of the rs1142888-G allele in Europeans (maximum |Z| > 4.3, Extended Data Fig. 6b)./p> 4.2, Fisher’s exact test, P < 2.3 × 10−6; Supplementary Note 6 and Supplementary Table 7d). Furthermore, among SARS-CoV-2-specific reQTLs, 28 reQTLs (5.3%) displayed signals of adaptation in East Asians starting 770–970 generations ago (around 25,000 years)—a timeframe associated with genetic adaptation at SARS-CoV-2-interacting proteins23 (OR relative to other populations = 2.6, Fisher’s exact test, P = 7.3 × 10−4; Fig. 4a and Extended Data Fig. 7a–c). An example is the immune mediator LILRB1, which has a SARS-CoV-2-specific reQTL (rs4806787) in pDCs (Extended Data Fig. 7d). However, the selection events making the largest contribution to the differentiation of SARS-CoV-2 responses in East Asia (top 5% PBS) began before this period (more than 970 generations ago, OR = 1.94, Fisher’s exact test, P = 0.019; Fig. 4b). For example, the rs1028396-T allele (80% frequency in East Asia versus 16–25% elsewhere), associated with a weaker response of SIRPA to SARS-CoV-2 in CD14+ monocytes, presents a selection signal beginning more than 45,000 years ago (Fig. 4b and Extended Data Fig. 7e). SIRPα inhibits infection by endocytic viruses, including SARS-CoV-241. These results suggest recurrent genetic adaptation targeting antiviral immunity over the last 50,000 years, contributing to present-day population differences in immune responses to SARS-CoV-2./p> 3). Each horizontal line represents a variant, sorted in descending order of time to onset of selection. The area shaded in purple highlights the period (770–970 generations ago) associated with genetic adaptation at host coronavirus-interacting proteins in East Asians23. Several immunity-related genes are highlighted. b, Allele frequency trajectories of two SARS-CoV-2 reQTLs (rs1028396 at SIRPA and rs11645448 at NOD2) in Africans (YRI, green), Europeans (CEU, yellow) and East Asians (CHS, purple). The full lines indicate the maximum a posteriori estimate of allele frequency at each epoch and shaded areas indicate the 95% CIs. The dendrograms show the estimated unrooted population phylogeny for each eQTL based on PBS (that is, the branch length between each pair of populations is proportional to −log10[1 − FST])./p> 1.2, one-sided permutation test, P < 2 × 10−2; Extended Data Fig. 8a and Supplementary Table 8a–c). Enrichment was strongest in SARS-CoV-2-stimulated CD16+ monocytes from Europeans, suggesting that archaic haplotypes altering myeloid responses have been preferentially retained in their genomes. Archaic haplotypes with eQTLs are generally present at higher frequencies compared with archaic haplotypes without eQTLs (Δf(introgressed allele) >3.2%, Student’s t-test, Padj < 8 × 10−3; Extended Data Fig. 8b and Supplementary Table 8d,e), even after adjustment for minor allele frequency (MAF) to ensure similar power for eQTL detection, supporting the adaptive nature of Neanderthal regulatory alleles./p> 3.8, respectively, one-sided resampling, P < 10−4), and a specific enrichment in reQTLs at severity loci (FE > 3.7, one-sided resampling, P < 3 × 10−3; Fig. 5a). This trend was observed across most cell lineages (Extended Data Fig. 10a). Colocalization analyses identified 40 genes at which there was a high probability of (r)eQTL colocalization with COVID-19 hits (posterior probability that both traits are linked to the same SNP (PPH4 ) > 0.8) and transcriptome-wide association studies (TWASs) linked predicted gene expression with COVID-19 risk for 30 of these genes (FDRTWAS < 0.01; Supplementary Table 9a). These included direct regulators of innate immunity, such as IFNAR2 in non-stimulated CD4+ T cells, IRF1 in non-stimulated NK and CD8+ T cells, OAS1 in lymphoid cells stimulated with SARS-CoV-2 and IAV, and OAS3 in SARS-CoV-2-exposed CD16+ monocytes (Fig. 5b and Extended Data Fig. 10b,c). These results support a contribution of immunity-related (r)eQTLs to COVID-19 risk./p> 0.8) and (2) presented positive selection signals (top 1% PBS, n = 13 eQTLs) or evidence of archaic introgression (n = 7 eQTLs), 14 of which regulate genes of which the expression is correlated with COVID-19 susceptibility and/or severity (FDRTWAS < 0.01) (Fig. 6). For example, two variants in high LD at DR1 (rs569414 and rs1559828, r2 > 0.73) displayed extremely high levels of population differentiation, probably due to selection outside Africa (DAF = 0.13 in Africa versus higher than 0.62 in Eurasia; Extended Data Fig. 10d). DR1 suppresses type I IFN responses49 and the selected alleles, which decrease COVID-19 severity, reduce DR1 expression in most immune cells (Fig. 6). Likewise, an approximately 39 kb Neanderthal haplotype, spanning the MUC20 locus in Eurasians, contains the rs2177336-T allele that increases MUC20 expression in SARS-CoV-2-stimulated cells, particularly for CD4+ T cells, and decreases COVID-19 susceptibility (Fig. 6). Together, these results reveal how past selection or Neanderthal introgression have impacted immune responses that contribute to present-day disparities in COVID-19 risk./p> 0.8) and presenting either strong population differentiation (top 1% PBS genome-wide) or evidence of Neanderthal introgression. a, Effects of the target allele on gene expression across immune lineages and stimulation conditions. b, Clinical and functional annotations of associated genes. c, Present-day population frequencies of the target allele. d, The effects of the target allele on COVID-19 risk (infection, hospitalization and critical state), colocalization probability and the lineage and condition in which gene expression most likely affects COVID-19 risk as detected by transcriptome-wide association (TWA) analyses. For expression or COVID-19 associations, the arrows indicate increases/decreases in expression or disease risk with each copy of the target allele, and the opacity reflects the strength of association (two-sided Student’s t-test −log10[P]). For the TWA analysis, the arrows indicate the effect of an increase in gene expression on the risk of COVID-19. In a and d, the arrow colours indicate stimulation conditions (non-stimulated (grey), SARS-CoV-2-stimulated (red), IAV-stimulated (blue)) and the background colour indicates the lineage (myeloid (pink), B (purple), CD4+ T (blue), CD8+ T (green), NK (light green)). For each eQTL, the target allele is defined as (1) the derived allele for highly differentiated eQTLs or (2) the allele that segregates with the archaic haplotype for introgressed eQTLs. When the ancestral state is unknown, the minor allele is used as a proxy for the derived allele. Note that, in some cases (for example, OAS1), the introgressed allele can be present in Africa, which is attributed to the reintroduction in Eurasia of an ancient allele by Neanderthals46. C, critical; H, hospitalized; R, reported./p>90% being born in either Cameroon or the Democratic Republic of Congo. For this study, 71 additional individuals of East Asian descent (ASH) were included, of whom 62 were retained after quality control (see the ‘scRNA-seq library preparation and data processing’ section). ASH individuals were recruited at the School of Public Health, University of Hong Kong, and were included in a community-based sero-epidemiological COVID-19 study (research protocol number JTW 2020.02). Inclusion for the study described here was restricted to nominally healthy ASH individuals (30 men and 41 women) aged between 19 and 65 years of age and seronegative for SARS-CoV-2. Samples were collected at the Red Cross Blood Transfusion Service (Hong Kong) where the PBMC fraction was isolated and frozen. Target sample sizes were determined to ensure >80% power for the detection of eQTLs with R2 higher than 0.2, at a P < 5 × 10−9 threshold./p> 5%) having an r2 > 0.8 for the correlation between observed and imputed genotypes (>95% concordance for 96% of common variants). After imputation, variants with a MAF < 1% or with a low predicted quality of imputation (that is, DR2 < 0.9) were excluded, yielding a final dataset of 13,691,029 SNPs for downstream analyses./p> 0.82, P < 7.6 × 10−13) (Supplementary Fig. 5b,c)./p> 0.2 were considered to be differentially expressed between populations (that is, ‘raw’ popDEGs). We adjusted for cellular composition within each lineage L by introducing into model (10) a set of variables \({({F}_{j\cdot })}_{j\in L}\) encoding the frequency in the PBMC fraction of each cell type j comprising the lineage (for example, naive, effector and regulatory subsets of CD4+ T cells)./p> 1% or |βa| < 0.2) and (2) displayed significant differences between the raw and adjusted effect sizes (|TΔβ| > 1.96) imputable to the effect of cellular composition./p> 5% in at least one of the three populations considered, resulting in a set of 10,711,657 SNPs, of which 4,164,060 were located <100 kb from a gene. We used MatrixEQTL (v.2.3)84 to map eQTLs in a 100 kb region around each gene and obtain estimates of eQTL effect sizes and their standard error. eQTL mapping was performed separately for each immune lineage/cell type and condition, based on rank-transformed gene expression values. eQTL analyses were performed adjusting for population, age, chromosomal sex, cell composition (within each lineage), as well as cell mortality and total number of cells in the sample, and a data-driven number of surrogate variables included to capture unknown confounders and remove unwanted variability. Specifically, for each immune lineage/cell type and condition, surrogate variables were obtained using the sva function from the sva R package (v.3.40.0)85 with option method=‛two-steps’, providing all other covariates as known confounders (mod argument). The number of surrogate variables to use in each lineage/cell type and condition was determined automatically based on the results from num.sv function with method=‛be’85./p>106. These runs were therefore discarded, and the associated eQTLs were assigned a null Z-score during FDR computation (see below). For each eQTL, the index SNP was defined as the SNP with the highest posterior inclusion probability (that is, the α parameter in the output of SuSiE) for that eQTL, and the 95% credible interval was obtained as the minimal set of SNPs S such that αs > 0.01 for all \(s\in S\) and \({\sum }_{s\in S}{\alpha }_{s} > 0.95\). Only eQTLs with a log-Bayes factor (lbf) > 3 were considered for further analyses./p> 10 in the lung and eQTLs with MAF > 5% in the GTEx dataset. We considered any eQTL with (1) P < 0.01 in the lung and (2) the same effect direction between lung and the lineage/cell type/condition in which it is the most significant in our study as replicated. As a comparison, we evaluated the amount of eQTLs that would be replicated when selecting SNPs at random, matching for MAF in GTEx (bins of 5%) and the distance between the eQTL index SNPs and the nearest gene (that is, bins of 0–1, 1–5, 5–10, 10–20, 20–50 and 50–100 kb), and computed the fold-enrichment in replicated eQTLs as the ratio between the observed and expected number of replicated eQTLs./p>100 kb). For each population and set of eQTLs, we defined the fold-enrichment (FE) in positive selection as the ratio of observed/expected values for mean PBS and extracted the mean and 95% confidence interval of this ratio across all resamplings. One-sided resampling P values were calculated as the number of resamplings with a FE > 1 divided by the total number of resamplings. Resampling P values were then adjusted for multiple testing by the Benjamini–Hochberg method./p> 3 to constitute evidence of selection and we inferred the onset of selection of a variant as the first generation in which |Z| > 3./p>150,000 were considered to be introgressed. This cut-off score has been shown to provide a good trade-off between power and accuracy based on simulations of introgression under realistic demographic scenarios98. For both calling methods (that is, CRF and S′), we used the recombination map from the 1KG Project Phase 3 data release56./p> 0.8) with at least two other aSNPs and, to exclude incomplete lineage sorting, comprising an LD block of >10 kb. This yielded a set of 100,345 high-confidence aSNPs (Supplementary Table 8a). We further categorized aSNPs as of Neanderthal origin, Denisovan origin or shared origin according to their presence/absence in the Vindija Neanderthal and Denisovan Altai genomes. Finally, we considered any site that was in high LD with at least one aSNP in the same population in which introgression was detected to be introgressed, and classified introgressed haplotypes as of Neanderthal origin, Denisovan origin or shared origin according to the most frequent origin of aSNPs in the haplotype. For introgressed SNPs, we defined the introgressed allele as (1) the allele rare or absent from individuals of African ancestries if the SNP was an aSNP; and (2) for non-aSNPs, the allele most frequently segregating with the introgressed allele of linked aSNPs. In each population, introgressed alleles with a frequency in the top 1% for introgressed alleles genome-wide were considered to present evidence of adaptive introgression./p>95% sharing at the lineage level). Within each cell type/stimulation condition, we considered the set of all (r)eQTLs for which the index SNP displayed at least a marginal association (Student’s t-test, P < 0.01) with gene expression. For each population and (r)eQTL set, we then grouped (r)eQTLs in high LD (r2 > 0.8), retaining a single representative per group, and counted the total number of (r)eQTLs for which the index SNP was in LD (r2 > 0.8) with an aSNP (that is, introgressed eQTLs). We then used PLINK (v.1.9) --indep-pairwise (with a 500 kb window, 1 kb step, an r2 threshold of 0.8, and a MAF > 5%)57 to define tag-SNPs for each population, and we determined the expected number of introgressed SNPs by resampling tag-SNPs at random with the same distribution for MAF, LD scores and distance to the nearest gene. We performed 10,000 resamplings for each (r)eQTL set and population. One-sided resampling-based P values were calculated as the frequency at which the number of introgressed SNPs among resampled SNPs exceeded the number of introgressed SNPs among (r)eQTLs. Resampling-based P values were then adjusted for multiple testing using the Benjamini–Hochberg method./p> 5% in each population (retaining a single representative per LD group) and compared the frequency of the introgressed allele with that of introgressed tag-SNPs genome-wide. We modelled r(Freq), the (rank-transformed) frequency of introgressed tag-SNPs according to the presence/absence of a linked eQTL (IeQTL), and the mean MAF of the SNP across the three populations (giving a higher power for eQTL detection)./p> 0.8) and 0 otherwise; \(\overline{{MAF}}\) is the mean MAF calculated separately for each population; α is the intercept of the model; β measures the difference in rank r(Freq) between eQTLs and non eQTLs; and γ is a nuisance parameter capturing the relationship between \(\overline{{MAF}}\) and r(Freq). Under this model, the difference in frequency between eQTLs and non-eQTLs can be tested directly in a Student’s t-test with \({{\mathscr{H}}}_{0}:\beta =0\)./p> 200, Fisher’s test, P = 4.2 × 10−40) between loci associated with susceptibility (C2) and severity (A2 or B2)8, 81 out of 105 COVID-19 associated eQTLs (at nominal P < 10−4) are associated specifically with either susceptibility (n = 19) or severity (n = 62), supporting the relevance of considering these traits separately in our analysis./p> 0.8 was considered to display significant colocalization./p>
Русский
English
Français
Deutsch
Español
Italiano
Português
日本語
한국어
Дом