네이처(2023)이 기사 인용
12,000회 액세스
293 알트메트릭
측정항목 세부정보
인간은 SARS-CoV-2 감염 후 상당한 개인간 임상적 다양성을 보이며1,2,3 이에 대한 유전적, 면역학적 근거가 해독되기 시작했습니다4. 그러나 SARS-CoV-2에 대한 면역 반응의 인구 차이의 정도와 동인은 아직 불분명합니다. 여기에서 우리는 SARS-CoV-2 또는 인플루엔자 A 바이러스로 자극된 다양한 조상의 건강한 기증자 222명의 말초 혈액 단핵 세포에 대한 단일 세포 RNA 서열 분석 데이터를 보고합니다. 우리는 SARS-CoV-2가 인플루엔자 A 바이러스에 비해 더 약하지만 더 이질적인 인터페론 자극 유전자 활동과 골수 세포에서 독특한 염증 유발 특징을 유도한다는 것을 보여줍니다. 바이러스에 대한 전사 반응은 주로 잠복 거대세포바이러스 감염과 관련된 림프구 분화 증가를 포함하여 세포 풍부도의 변화에 의해 주도되는 현저한 인구 차이를 나타냅니다. 발현 정량적 특성 유전자좌 및 중재 분석은 특정 유전자좌에 강력한 영향을 미치는 유전적 변이와 함께 면역 반응의 인구 불균형에 대한 세포 구성의 광범위한 효과를 보여줍니다. 또한, 우리는 자연 선택이 특히 동아시아인의 SARS-CoV-2 반응과 관련된 변종의 경우 면역 반응의 인구 차이를 증가시켰음을 보여주고, 네안데르탈인 유전자 이입이 반응과 같은 면역 기능을 변화시킨 세포 및 분자 메커니즘을 문서화합니다. 골수 세포에서 바이러스로. 마지막으로, 동일 위치화 및 전사체 전체 연관성 분석을 통해 SARS-CoV-2와 코로나19 중증도에 대한 면역 반응의 유전적 기초 사이의 중복이 밝혀져 현재 코로나19 위험 격차에 기여하는 요인에 대한 통찰력을 제공합니다.
코로나19 팬데믹의 주목할만한 특징은 무증상 감염부터 치명적인 질병에 이르기까지 SARS-CoV-2에 감염된 개인 사이에 상당한 임상적 차이가 있다는 점입니다1,2,3. 위험 요인에는 고령1, 남성 성별5, 동반질환6 및 숙주 유전학4,7,8이 포함됩니다. 또한, 선천적 오류 또는 I형 인터페론에 대한 자가 항체 중화를 포함한 선천적 면역의 변화9,10,1112,13,14는 임상 결과의 변화에 영향을 미치며, 역학적 및 유전적 데이터는 인구 집단 간의 차이를 시사합니다6,7,15,16 . 이는 면역 문제에 대한 전사 반응의 조상 관련 차이에 대한 보고와 함께17,18,19 전 세계 인구에 걸쳐 SARS-CoV-2에 대한 면역 반응의 규모와 변동 동인에 대한 조사가 필요합니다.
병원체가 부과하는 선택 압력은 인간 진화 과정에서 가장 중요했습니다. 선택적 청소 또는 고대 혼합물을 통한 RNA 바이러스에 대한 인간의 적응은 인구 유전적 분화의 원천으로 확인되었으며18,21,22 적응 신호는 동아시아인의 코로나바이러스 상호 작용 단백질에서 보고되었습니다23,24. 또한 고대 유전자 이입과 면역 사이의 연관성에 대한 증거도 있으며25, 네안데르탈인 일배체형은 코로나19 심각도와 관련이 있습니다26,27. 그러나 SARS-CoV-2에 대한 면역 반응에 대한 자연 선택과 고대 혼합물의 영향은 아직 조사되지 않았습니다.
우리는 중앙아프리카, 서유럽 및 동아시아 출신 개인의 말초 혈액 단핵 세포(PBMC)를 SARS-CoV-2에 노출시키고, 비교를 위해 인플루엔자 A 바이러스(IAV)에 노출시켜 이러한 문제를 해결했습니다. 단일 세포 RNA 염기서열 분석(scRNA-seq)과 정량적 및 집단 유전학 접근 방식을 결합하여 SARS-CoV-2에 대한 면역 반응의 인구 차이에 대한 환경적, 유전적 동인을 설명합니다.
우리는 세 가지 지리적 위치(중앙 아프리카, n = 80 남성, 서유럽, n = 80)에서 유래한 222명의 SARS-CoV-2 미경험 기증자로부터 PBMC에 대한 scRNA-seq 분석을 수행하여 SARS-CoV-2 및 IAV에 대한 전사 반응을 특성화했습니다. 남성, 동아시아, n = 여성 36명, 남성 26명) 및 유전적 조상이 다릅니다(보충 그림 1 및 보충 표 1). PBMC는 모의 대조군(비자극), SARS-CoV-2(조상 균주, BetaCoV/France/GE1973/2020) 또는 IAV로 6시간 동안 처리되었습니다(보충 참고 1, 보충 그림 2 및 보충 표 2). (H1N1/PR/8/1934). 우리는 1백만 개 이상의 고품질 단일 세포 전사체를 캡처했습니다(그림 1a, 보충 그림 3 및 보충 표 3a). 전사체 기반 클러스터를 시퀀싱(CITE-seq; 방법)에 의한 전사체 및 에피토프의 세포 색인화와 결합함으로써 골수성, B, CD4+ T, CD8+ T 및 자연살해(NK) 면역 계통에 걸쳐 22개의 세포 유형을 정의했습니다(그림 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
日本語
Русский
홈페이지