Natura (2023) Cita questo articolo
12k accessi
293 Altmetrico
Dettagli sulle metriche
Gli esseri umani mostrano una sostanziale variabilità clinica interindividuale dopo l’infezione da SARS-CoV-21,2,3, la cui base genetica e immunologica ha iniziato a essere decifrata4. Tuttavia, l’entità e i fattori che determinano le differenze della popolazione nelle risposte immunitarie alla SARS-CoV-2 rimangono poco chiari. Qui riportiamo i dati di sequenziamento dell’RNA di singole cellule per le cellule mononucleate del sangue periferico, provenienti da 222 donatori sani di diverse origini, che sono state stimolate con SARS-CoV-2 o virus dell’influenza A. Mostriamo che SARS-CoV-2 induce un’attività genetica stimolata dall’interferone più debole, ma più eterogenea, rispetto al virus dell’influenza A, e una firma proinfiammatoria unica nelle cellule mieloidi. Le risposte trascrizionali ai virus mostrano marcate differenze nella popolazione, principalmente guidate da cambiamenti nell'abbondanza cellulare, inclusa una maggiore differenziazione linfoide associata all'infezione latente da citomegalovirus. Le analisi dei loci dei tratti quantitativi di espressione e delle analisi di mediazione rivelano un ampio effetto della composizione cellulare sulle disparità della popolazione nelle risposte immunitarie, con varianti genetiche che esercitano un forte effetto su loci specifici. Inoltre, mostriamo che la selezione naturale ha aumentato le differenze di popolazione nelle risposte immunitarie, in particolare per le varianti associate alla risposta SARS-CoV-2 negli asiatici orientali, e documentiamo i meccanismi cellulari e molecolari attraverso i quali l’introgressione dei Neanderthal ha alterato le funzioni immunitarie, come la risposta delle cellule mieloidi ai virus. Infine, le analisi di colocalizzazione e di associazione a livello del trascrittoma rivelano una sovrapposizione tra la base genetica delle risposte immunitarie alla gravità della SARS-CoV-2 e della COVID-19, fornendo approfondimenti sui fattori che contribuiscono alle attuali disparità nel rischio di COVID-19.
Una caratteristica notevole della pandemia di COVID-19 è la sostanziale variazione clinica tra gli individui infetti da SARS-CoV-2, che va dall’infezione asintomatica alla malattia fatale1,2,3. I fattori di rischio comprendono l'età avanzata1, il sesso maschile5, le comorbilità6 e la genetica dell'ospite4,7,8. Inoltre, la variazione dell'immunità innata9,10,11, inclusi errori congeniti o autoanticorpi neutralizzanti contro gli interferoni di tipo I12,13,14, contribuisce alla variazione dell'esito clinico e i dati epidemiologici e genetici suggeriscono differenze tra le popolazioni6,7,15,16 . Ciò, insieme alle segnalazioni di differenze legate agli antenati nelle risposte trascrizionali alle sfide immunitarie17,18,19, richiede indagini sull’entità e sui fattori di variazione nelle risposte immunitarie alla SARS-CoV-2 tra le popolazioni di tutto il mondo.
Le pressioni selettive imposte dagli agenti patogeni sono state fondamentali durante l’evoluzione umana20. L'adattamento umano ai virus a RNA, attraverso spazzamenti selettivi o mescolanze arcaiche, è stato identificato come fonte di differenziazione genetica della popolazione18,21,22 e segnali di adattamento sono stati segnalati nelle proteine che interagiscono con il coronavirus negli asiatici orientali23,24. Esistono anche prove di collegamenti tra introgressione arcaica e immunità25, con aplotipi di Neanderthal associati alla gravità del COVID-1926,27. Tuttavia, gli effetti della selezione naturale e della mescolanza arcaica sulle risposte immunitarie alla SARS-CoV-2 restano da studiare.
Abbiamo affrontato queste domande esponendo le cellule mononucleari del sangue periferico (PBMC) di individui di origine dell’Africa centrale, dell’Europa occidentale e dell’Asia orientale alla SARS-CoV-2 e, per confronto, al virus dell’influenza A (IAV). Combinando il sequenziamento dell’RNA a cellula singola (scRNA-seq) con approcci quantitativi e di genetica della popolazione, delineiamo i fattori ambientali e genetici delle differenze della popolazione nelle risposte immunitarie alla SARS-CoV-2.
Abbiamo caratterizzato le risposte trascrizionali a SARS-CoV-2 e IAV eseguendo l'analisi scRNA-seq dei PBMC di 222 donatori naive alla SARS-CoV-2 provenienti da tre località geografiche (Africa centrale, n = 80 maschi; Europa occidentale, n = 80 maschio; Asia orientale, n = 36 femmine e 26 maschi) e con diversi antenati genetici (Figura 1 supplementare e Tabella supplementare 1). I PBMC sono stati trattati per 6 ore (Nota supplementare 1, Figura 2 supplementare e Tabella supplementare 2) con un finto controllo (non stimolato), SARS-CoV-2 (ceppo ancestrale, BetaCoV/Francia/GE1973/2020) o IAV (H1N1/PR/8/1934). Abbiamo catturato oltre 1 milione di trascrittomi unicellulari di alta qualità (Fig. 1a, Figura 3 supplementare e Tabella 3a supplementare). Combinando cluster basati su trascrittoma con l'indicizzazione cellulare di trascrittomi ed epitopi mediante sequenziamento (CITE-seq; Metodi), abbiamo definito 22 tipi di cellule attraverso lignaggi immunitari mieloidi, B, CD4 + T, CD8 + T e natural killer (NK) (Fig. 1b , Figura supplementare 4 e Tabella supplementare 3b-d). Dopo l'esposizione al virus, la maggior parte dei tipi di cellule ha mostrato cambiamenti moderati in abbondanza, con i cambiamenti più forti osservati nella linea mieloide dopo il trattamento con IAV (Nota supplementare 2 e Tabella supplementare 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>
Italiano
English
Français
Deutsch
Español
Português
日本語
한국어
Русский
Casa