Genome-wide association analysis and multi-omic Mendelian randomization study exploring the immune response in vitiligo

doi:10.21203/rs.3.rs-5010438/v1

Download PDF

Article

Genome-wide association analysis and multi-omic Mendelian randomization study exploring the immune response in vitiligo

https://doi.org/10.21203/rs.3.rs-5010438/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

The etiology of vitiligo involves immune dysregulation, though its precise genetic underpinnings remain incompletely understood. Here, we employ a multi-omics approach, integrating genome-wide association studies (GWAS) and Mendelian randomization (MR), to investigate the association between immune response-related genes and vitiligo. We conducted a comprehensive meta-analysis of three GWAS encompassing generalized vitiligo cohorts from Jin et al., the FinnGen cohort, and the UK Biobank to identify novel genetic determinants of vitiligo susceptibility. Using the results from this meta-analysis, we employed Mendelian randomization (MR) and summary data-based MR (SMR) to discern immune response genes having a putative causal relationship with vitiligo on the level of plasma proteome. Additionally, we integrated summary data on immune response methylation and expression abundance levels for multi-omics validation. Further exploration involved assessing the differential abundance of immune response genes at the single-cell transcriptomic level and tracking their expression dynamics during cellular differentiation. Our meta-analysis unveiled 25 genome-wide significant vitiligo risk variants, six of which were previously unreported. Notably, the predicted protein levels of eight genes displayed associations with vitiligo, encompassing the methylation levels of CD160 and TYRO3, as well as the gene expression level of CD160. These genes were predominantly expressed in T cells and mononuclear phagocytes within vitiligo skin lesions, exhibiting distinct expression patterns and temporal changes across various disease states. Through the integration of GWAS and multi-omics MR approaches, this study identifies several immune response genes implicated in vitiligo pathogenesis, offering promising targets for future therapeutic and preventive strategies.

Biological sciences/Immunology/Autoimmunity

Health sciences/Diseases/Immunological disorders

Health sciences/Diseases/Skin diseases/Vitiligo

Biological sciences/Biological techniques/Genomic analysis/Genome wide association studies

Biological sciences/Biological techniques/Proteomic analysis

Vitiligo

Immune response

Meta-analysis

Mendelian randomization

Multi-omics

Vitiligo is a common acquired depigmenting disorder characterized by the selective loss of melanocytes. The pathogenesis of vitiligo is multifactorial, involving autoimmunity, oxidative stress, mitochondrial dysfunction, and genetic predisposition among other factors^{[1], [2], [3]}. Overactivation of the immune response plays a crucial role in vitiligo. Under the stimulus of localized oxidative stress, autoreactive cytotoxic CD8 + T cells, keratinocytes, fibroblasts, and other cells participate in the destruction of melanocytes by secreting cytokines such as interferon-gamma (IFN-γ) and CXCL9^{[1], [4], [5]}. Despite numerous reports on immune response-related genes in vitiligo, no studies have comprehensively and systematically established their potential causal relationship with the disease.

Mendelian randomization (MR) uses genetic variants as instrumental variables for exposures^[6]. This approach helps reduce residual confounding, strengthen causal inference, and minimize the potential for reverse causation^[7]. Large-scale GWAS and molecular quantitative trait loci (QTL) data have been utilized to identify causal relationships between diseases and methylation, gene expression, and protein abundance^{[8], [9]}.

This study conducted a meta-analysis of three GWAS datasets from different sources to identify novel vitiligo-associated risk loci and annotate genes. By integrating multi-omics QTL, utilizing MR and summary data-based MR (SMR), we explored potential associations between immune response-related methylation, gene expression, as well as protein abundance and vitiligo. Single-cell data validated the variations of relevant genes at vitiligo lesions and their cell-specific enrichment.

Study design

Figure 1 illustrates the overall design of the study. We performed a meta-analysis from three publicly available vitiligo GWAS for gene annotation and functional enrichment of re-identified risk loci. The genetic correlation between the results obtained from the meta-analysis and common risk factors for vitiligo was calculated. Subsequently, we utilized pQTL data from seven large-scale proteome studies to extract instrumental variables for immune response genes and performed Mendelian randomization (MR) and summary data-based MR (SMR) on vitiligo. We also extracted eQTL and mQTL data for immune response genes for validation at different multi-omics levels and repeated MR in two subtypes of vitiligo (early-onset and late-onset). Finally, we performed single-cell transcriptomic analysis on identified causal candidate proteins to investigate the biological credibility of our genetic findings.

Data sources of vitiligo data

Genetic associations with vitiligo were derived from GWAS conducted by Jin et al. on generalized vitiligo (4,680 cases and 39,586 controls) ^[10], the FinnGen cohort (131 cases and 207,482 controls)^[11], and the UK Biobank (95 cases and 337,064 controls)^[12]. All participants in the three GWAS were of European ancestry. Summary statistics for the two subtypes of vitiligo GWAS were obtained from Jin et al., who stratified vitiligo by age, with the early-onset subgroup consisting of 704 cases and 9,031 controls, and the late-onset subgroup consisting of 1,467 cases and 19,156 controls^[13].

Datasets of protein, expression, and methylation quantitative trait loci

The largest set of immune response genes, comprising 2,775 human immune response-related genes, was extracted from the MSigDB database. Protein quantitative trait loci (pQTL) data were obtained from seven published plasma proteome GWAS studies ^{[14], [15], [16], [17], [18], [19],[20]}. For duplicated proteins in the study, those with the maximum sum of R² (R² = 2×EAF×(1-EAF)×beta²) ^[21] were selected. Expression quantitative trait loci (eQTL) data were obtained from the eQTLGen Consortium, including 31,684 individuals^[22]. Methylation quantitative trait loci (mQTL) data were derived from a meta-analysis by McRae et al. of two European cohorts, the Brisbane Systems Genetics Study (n = 614) and the Lothian Birth Cohort (n = 1,366). Individual methylation probes were normalized using a generalized linear model^[23]. Tissue-specific expression eQTL data were retrieved from the Genotype-Tissue Expression (GTEx) portal. The GTEx v8 dataset includes 838 donors and 17,382 samples from 52 tissues and two cell lines, evaluating tissue-specific expression of target genes^[24]. In this analysis, we utilized eQTL data from skin tissue.

Meta-analysis of vitiligo GWAS

We conducted a sample size-weighted meta-analysis of vitiligo GWAS data from three sources using METAL^[25]. This encompassed a total of 589,038 individuals. Variants with a minor allele frequency (MAF) below 0.5% were removed, resulting in 14,608,392 associations. Genome-wide significant single nucleotide polymorphisms (SNPs) (P < 5.0 × 10^− 8) were identified and linkage disequilibrium (LD, r² < 0.001) pruned, isolating significant genome-wide signals for vitiligo.

Annotation of risk variants

We employed FUMA (https://fuma.ctglab.nl/) to annotate our findings. Three gene-mapping strategies were utilized in FUMA to identify genes associated with SNPs included in the MR analysis: positional mapping, gene mapping based on eQTLs, and gene mapping based on chromatin interactions. Top genes at each locus identified by Polygenic Priority Score (PoPS) were prioritized for our GWAS significant loci (p < 5 × 10^− 8).

The PoPS method is a novel gene prioritization approach that integrates GWAS summary statistics with gene expression, biological pathways, and predicted protein interactions to identify causal genes. We applied PoP scores because it has demonstrated higher predictive confidence in nominating causal genes at non-coding GWAS loci compared to other similarity-based or locus-based methods. PoPS analysis returned scores for a total of 18,383 genes across our GWAS datasets. We annotated our GWAS loci using Ensembl genes within a 500 Kb window and selected genes with the highest PoP scores at each locus as the prioritized genes.

Pathway enrichment analysis

We conducted enrichment analyses to identify biological pathways associated with vitiligo risk loci that met the GWAS meta-analysis p-value threshold. Two strategies were used to determine the genes for enrichment analysis. First, genes annotated by FUMA were used. Additionally, for each variant, cis-eQTLs within a 1 Mb region related to the variant were identified from peripheral blood and skin tissue in GTEx v8, and all genes with p < 5 × 10^− 4 were extracted. All retrieved genes were merged into a single gene set, which was then used to identify enriched pathways using Gene Ontology and KEGG pathway analyses.

Associations of vitiligo risk variants with risk factors for vitiligo

For genetic variants that surpassed the vitiligo GWAS threshold (p < 5 × 10^− 8), we assessed their genetic association with seven vitiligo risk factors identified from available GWAS data: depression, anxiety, insomnia, rheumatoid arthritis, autoimmune thyroiditis, smoking dependency, and alcohol dependency. GWAS data for these vitiligo risk factors were sourced from the Psychiatric Genomics Consortium (PGC) and the FinnGen cohort.

We used p < 0.005 (0.05 divided by the number of tested vitiligo secondary traits) to account for multiple testing. For associations that met the p-value threshold, we further evaluated whether the direction of the association with vitiligo risk factors was consistent with the association with vitiligo.

LD score regression

We utilized LD Score regression (LDSC) to estimate the genetic correlation (rg) between vitiligo and its risk factors. GWAS summary statistics were filtered using the HapMap3 reference panel. LDSC assesses the relationship between test statistics and linkage disequilibrium to quantify the contribution of inflation from a true polygenic signal or bias^[26]. This method evaluates genetic correlation from GWAS summary statistics without being affected by sample overlap^[27]. The genetic covariance normalized by SNP heritability represents genetic correlation. A p-value of < 0.005 was considered statistically significant, while 0.005 < p < 0.05 indicated suggestive evidence of genetic correlation^[28].

MR analysis

We screened for immune response proteins from the pQTL dataset. Selection criteria for pQTL instruments were as follows: (1) pQTLs should be significantly associated with any exposure (p < 1×10^− 8). (2) To avoid linkage disequilibrium (LD) effects, pQTLs should be independent (R² < 0.001 and genetic distance > 10,000 kb). (3) pQTLs should not be associated with any confounders or outcomes. The strength of genetic instruments was estimated using R² and the F statistic (F = R² × (N − 2) / (1 - R²))^[21], where R² is the proportion of protein level variability explained by each genetic instrument. Instruments were further restricted to cis-QTLs, defined as being within a 1 Mb region of the transcription start site of the protein-coding gene. Regions outside this area were defined as trans-QTLs^[16].

The TwoSampleMR package (version 4.2.2)^[29] was used to conduct MR analyses. For methylation genes, expression genes, and proteins with only one instrument, the Wald ratio method was used to estimate the log odds change in vitiligo risk associated with an SD increase in exposure levels represented by the instrumental variable. For proteins with multiple instruments, the inverse variance weighted (IVW) method was employed to obtain MR effect estimates. Heterogeneity tests were conducted to assess the heterogeneity of genetic instruments based on Q statistics. Additional analyses, including simple mode, weighted mode, weighted median, and MR-Egger, were performed to consider horizontal pleiotropy^[30]. MR-Egger results were used only when the intercept indicated the presence of horizontal pleiotropy. A false discovery rate (FDR) < 0.05 was considered significant. In the analysis of the two subtypes of vitiligo, replication MR analyses were conducted to validate the identified proteins. FDR < 0.05 was defined as a significant level. Additionally, we evaluated the correlation of identified proteins with previously reported vitiligo risk factors to explore whether the identified proteins have a potential mediating causal relationship with vitiligo. p < 0.05 was defined as a significant level.

Since no heterogeneity evidence based on the IVW model's Cochran Q statistic or MR-Egger intercept test was observed for individual cis-pQTLs, we conducted sensitivity analyses using a larger number of mixed pQTLs (cis + trans) for the proteins. p < 0.05 was defined as a significant level. We used the MR Steiger directionality test to check if the exposure has a directional causal relationship with the outcome.

SMR analysis

We conducted summary data-based Mendelian randomization (SMR) analysis to validate the causal relationship between the protein abundance of immune response genes and vitiligo^[31]. We screened for immune response genes across various QTL datasets, obtaining methylated genes, expressed genes, and proteins with available instruments from mQTL, eQTL, and pQTL datasets, respectively. Top associated cis-QTLs were selected by considering a window around the respective genes (± 1 Mb) and using a P-value threshold of 5.0 × 10^− 8. SNPs with allele frequency differences exceeding a specified threshold were excluded between any paired datasets, including LD reference samples, QTL summary data, and outcome summary data. SMR was implemented using the SMR software tool (SMR v1.3.1). P-values were adjusted using the Benjamini-Hochberg method to control the false discovery rate (FDR) at 0.05.

Single-cell transcriptomic expression analysis

Single-cell RNA-seq data from vitiligo lesion tissues and healthy control tissues in the Genome Sequence Archive (GSA) were used to further evaluate the cell type-specific expression of immune response genes with potential causal effects.

Data quality control, normalization, cell clustering, and dimension reduction were conducted using Seurat (version 4.4.0)^[32]. The top 3000 highly variable genes (HVGs) per sample were analyzed post-normalization using the variance stabilizing transformation (VST). Technical or batch effects were mitigated using the R package harmony^[33]. A principal component analysis (PCA) matrix incorporating the top 50 components was constructed using the RunPCA function. A nearest neighbor graph was constructed using the FindNeighbors function, followed by cell clustering using the FindClusters function with a resolution parameter set at 0.8. Clustered cells were projected into two-dimensional space for visualization using the non-linear dimension reduction methods RunUMAP and RunTSNE from the Seurat package. Cell annotation was performed using the expression of canonical marker genes. Trajectory analysis of T cell subpopulations and exploration of pseudo-temporal changes in the expression levels of immune response genes were conducted using the R package monocle^[34].

Genome-wide meta-analysis identifies significant vitiligo loci

A meta-analysis of vitiligo GWAS data from studies by Jin et al., the FinnGen cohort, and the UK Biobank (Supplementary Table 1) identified variants at a genome-wide significance level (p < 5 × 10^− 8) (Fig. 1). The quantile-quantile (Q-Q) plot of the meta-analysis is shown in Supplementary Fig. 1. We annotated potential causal genes for newly discovered vitiligo variants and conducted enrichment analysis of these variants. Additionally, we explored the associations and genetic correlations between these variants and seven vitiligo risk factors.

After quality control, the meta-analysis of vitiligo GWAS data from three studies yielded results for 14,608,392 genetic variants associated with vitiligo. Twenty-five variants showed genome-wide significant signals, six of which were located more than 100 Kb from previously reported index variants (Fig. 2 and Table 1).

Gene mapping of significant vitiligo genetic variants

We performed fine-mapping using GWAS summary statistics (Supplementary Fig. 2). Using FUMA, we mapped genes within a 500 Kb region of the index SNPs and identified genes with the highest Polygenic Priority Score (PoPS) in this region (Table 1). Notably, rs201529506 did not index any annotated genes within the 500 Kb region.

Table 1

Loci identified in the GWAS meta-analysis of vitiligo
rsID	Effect allele	Other allele	EAF	N	Beta	SE	P value	PoPs gene	PoPs score
Previously reported loci
rs11453006	D	I	0.9	44266	0.446	0.060	1.17E-15	TUBB3	0.289
rs35059390	D	I	0.53	44266	0.191	0.030	3.29E-11	TG	0.928
rs34976449	D	I	0.56	44266	0.199	0.030	5.38E-12	RERE	0.551
rs11434358	D	I	0.61	44266	0.174	0.030	1.75E-09	TICAM1	0.629
rs9981980	T	G	0.67	44266	0.285	0.030	2.51E-22	UBASH3A	0.222
rs148136154	T	C	0.88	44266	-0.315	0.040	6.89E-15	CD80	0.455
rs58473363	D	I	0.6	44266	-0.207	0.030	9.84E-13	CD44	0.941
rs113759362	A	G	0.98	44266	0.577	0.090	1.44E-11	HLA-DRB1	0.805
rs9279765	D	I	0.76	44266	-0.357	0.040	1.81E-23	HLA-DRB1	0.805
rs60949565	D	I	0.94	44266	-0.358	0.050	5.27E-12	HERC2	0.100
rs61710173	D	I	0.52	44266	-0.207	0.030	6.04E-13	CCR6	0.371
rs78567876	A	T	0.55	44266	0.531	0.030	1.34E-75	HLA-DRB1	0.805
rs35663695	D	I	0.69	44266	0.315	0.030	2.88E-21	NOX4	-0.123
rs28520266	C	G	0.5	44266	-0.211	0.030	2.43E-13	IFIH1	0.465
rs141105463	D	I	0.58	44266	0.223	0.030	6.75E-14	FOXP1	0.624
rs113284964	D	I	0.67	44266	0.199	0.030	2.08E-11	RHOH	1.424
rs10676540	D	I	0.79	44266	0.223	0.040	2.90E-09	TOB2	0.365
rs145282249	A	T	0.55	44266	-0.278	0.030	7.74E-22	LPP	0.588
rs67822292	T	G	0.9988	337159	0.004	0.001	3.01E-09	NA	NA
Novel loci
rs114259595	A	G	0.62	44266	0.166	0.030	8.27E-09	BACH2	0.595
rs79838180	T	G	0.62	44266	0.236	0.030	5.61E-15	ANKRD11	0.353
rs5845323	D	I	0.57	44266	-0.285	0.030	3.77E-23	NCF4	0.462
rs142543923	D	I	0.95	44266	-0.392	0.060	5.77E-13	CCDC88A	0.208
rs201529506	A	C	0.6348	544772	-0.001	0.000	4.09E-08	NA	NA
rs2017445	A	G	0.6185	251879	0.265	0.029	2.76E-09	BLOC1S1	0.321

EAF, effect allele frequency; SE, standard error; PoPs, Polygenic Priority Score.

Firstly, we conducted enrichment analysis using genes within the 500 Kb region of the index SNPs. GO analysis showed biological pathways enriched in immune-related pathways such as regulation of leukocyte cell-cell adhesion and antigen processing and presentation. KEGG analysis indicated that genes are associated with autoimmune diseases like autoimmune thyroid disease and Type I diabetes mellitus (Fig. 3, A and B). Subsequently, using the 25 vitiligo GWAS variants as cis-eQTLs in GTEx v8, we identified 38 significantly associated unique genes (p < 5 × 10^− 4) and performed further GO and KEGG analyses. The results revealed biological pathways including neutrophil aggregation, leukocyte aggregation and chemotaxis, and KEGG pathways like complement and coagulation cascades and the IL-17 signaling pathway (Fig. 3, C and D).

PoPS-suggested genes were used to characterize different variants. rs9279765 (HLA-DRB1) showed the most significant association values with autoimmune thyroiditis (p = 4.9 × 10^− 16) and rheumatoid arthritis (p = 2.5 × 10^− 9). Additionally, rs2017445 (BLOC1S1) was associated with rheumatoid arthritis (p = 4.3 × 10^− 4), rs10676540 (TOB2) with insomnia (p = 2.6 × 10^− 3), and rs5845323 (NCF4) with autoimmune thyroiditis (p = 3.2 × 10^− 3). Vitiligo variants showed the most associations with autoimmune thyroiditis, and the direction of association was consistent with the risk of vitiligo. No associations were found between variants and depression, anxiety, alcohol dependency, or smoking dependency (Fig. 3E).

Using LDSC regression analysis, we assessed the genetic correlation between seven vitiligo risk factors and vitiligo. As shown in Supplementary Fig. 2, LDSC indicated no suggestive correlation between the seven vitiligo risk factors and vitiligo (p > 0.05).

MR and SMR identify immune response genes for vitiligo

From the cis-pQTLs of immune response genes extracted from seven plasma proteome GWAS cohorts, 381 protein-SNP pairs were utilized as instrumental variables (Table S3) in a two-sample MR analysis with GWAS meta-analysis. Using the Wald ratio or IVW methods, 30 immune response proteins were significantly associated with vitiligo risk (FDR < 0.05) (Table S4, Fig. 4, A and B). To further validate the observed results, SMR analysis was also performed on the included immune response proteins, revealing 100 proteins significantly associated with vitiligo risk (FDR < 0.05) (Table S8). Eight proteins passed both MR and SMR tests. Elevated predicted levels of CXCL6 were associated with increased vitiligo risk, whereas other proteins (CDH17, LEAP2, CRP, IL9, CFHR5, CD160, and TYRO3) were negatively correlated with vitiligo risk (Table 2). SMR plots and effect diagrams for these eight proteins are displayed in Supplementary Figs. 3 and 4.

Table 2

Summary results from Mendelian randomization (MR) and SMR for 8 proteins passing all tests
Protein	Method	Number of SNPs	MR			SMR
			Beta	SE	FDR	Beta	SE	FDR
CD160	Wald ratio	1	-0.0092	0.0008	0.0000	0.0152	0.0061	0.0431
CRP	Wald ratio	1	-0.0125	0.0034	0.0397	0.0176	0.0027	0.0000
CFHR5	Wald ratio	1	-0.0027	0.0003	0.0000	0.0133	0.0003	0.0000
CXCL6	Wald ratio	1	0.0039	0.0004	0.0000	0.0053	0.0005	0.0000
LEAP2	Wald ratio	1	-0.0133	0.0005	0.0000	-0.0107	0.0018	0.0000
IL9	Wald ratio	1	-0.0192	0.0033	0.0000	-0.0308	0.0064	0.0000
CDH17	Wald ratio	1	-0.0113	0.0014	0.0000	-0.0114	0.0015	0.0000
TYRO3	Wald ratio	1	-0.0093	0.0003	0.0000	0.0200	0.0082	0.0485

SNP, single nucleotide polymorphisms; MR, Mendelian randomization; SMR, summary-based Mendelian randomization; SNP, single nucleotide polymorphisms; SE, standard error; FDR, false discovery rate.

In sensitivity analysis (limited to MR), mixed pQTLs (cis + trans) for these eight proteins were reselected for MR. Excluding CRP and IL9, the remaining six maintained significant associations with vitiligo risk (FDR < 0.05) using the Wald ratio or IVW methods, consistent with the direction in cis MR (Table S5). In the cis + trans MR, heterogeneity evidence based on the Cochran Q statistic in the IVW model was observed for CRP, IL9, and CDH17 (P_{heterogeneity} > 0.05). No horizontal pleiotropy was detected through MR-Egger intercept tests. All proteins' relationship with vitiligo was confirmed by Steiger directionality tests (Table S5).

In subtype analyses of vitiligo, all eight proteins passed replication MR analysis, showing associations with both early-onset and late-onset types (FDR < 0.05) (Tables S6 and S7). The results for CRP, CHFR5, CD160, and TYRO3 were directionally consistent across vitiligo and its subtypes.

Integrating multi-omics and tissue-specific evidence

Our aim was to integrate multi-omics evidence to identify immune response genes among the eight candidate proteins significantly associated with vitiligo risk at gene expression and methylation levels.

CD160 gene expression correlated with vitiligo risk at an FDR < 0.05 level. Higher predicted levels of CD160 (β = 0.010, se = 0.002) were positively associated with vitiligo risk (Table S9).

Results for the causal effects of methylation in immune response genes on vitiligo are stored in Table S10. After correcting for multiple tests, we identified nine CpG sites near two immune response genes, including CD160 (cg06095777, cg11743829, cg25221984, cg08614201, cg20975414, and cg12832565), and TYRO3 (cg24374636, cg15330654, and cg25906537) (Table S). The direction of effect estimates at different CpG sites within the same gene was not always consistent. For example, an increase in predicted methylation of CD160 at cg25221984 by one standard deviation was associated with decreased vitiligo risk (β = -0.003, se = 0.001), whereas increases at other CpG sites were associated with increased risk (Table S10).

We further explored whether the identified proteins in blood could replicate similar causal relationships in skin tissue. No significant associations were found between identified immune response genes and vitiligo risk in skin tissue. However, we identified several skin tissue-specific immune response genes associated with vitiligo risk (FDR < 0.05), which also met the FDR < 0.05 level in plasma protein MR or SMR analyses, including GNLY, IL1R2, TDGF1, MASP1, C6, SFTPD, PLCG2, and MIF (Table S11).

Phenome-wide association analysis

We evaluated the causal relationships between eight MR proteomic genes and seven vitiligo risk factors to explore potential mediation in their association with vitiligo risk.

We observed that CFHR5 (β = -0.0008, se = 0.0004) and IL9 (β = 0.004, se = 0.002) were associated with anxiety or panic attacks, and the direction of association was consistent. No causal relationships were observed between the remaining proteins and the seven vitiligo risk factors (Fig. 4C).

Additionally, for the eight MR proteomic genes, we searched the EpiGraphDB database and identified 62 reported associations with medical terms related to vitiligo, except for CDH17, IL9, and CD160, which were not listed in the database. The most frequently occurring terms were severe depression, bipolar disorder, and schizophrenia (Table S12).

Cell type-specific expression in skin tissue

We downloaded scRNA-seq data from five healthy controls and ten vitiligo patients from the GSA, totaling 48,887 cells after merging the data. Clusters were annotated as eight cell types based on gene markers, including keratinocytes, melanocytes, fibroblasts, endothelial cells, smooth muscle, T & NK cells, langerhans cells, and mononuclear phagocytes (M. phagocytes) (Fig. 5, A and B). The relative abundance of the eight main cell groups across 15 samples and two groups is shown in Fig. 5C. We found that the proportion of immune cells increased in vitiligo patients compared to healthy controls, while the numbers of keratinocytes and melanocytes decreased and the rest increased.

In the top five genes with the highest PoP scores mapped from newly discovered vitiligo GWAS variants, as well as the eight proteins related to vitiligo, CRP, CFHR5, and IL9 were not detected in the scRNA-seq data. BACH2, CCDC88A, BLOC1S1, NCF4, ANKRD11, CD160, LEAP2, CDH17, and TYRO3 were expressed in the skin tissue, while CXCL6 was not detected. Compared to healthy controls, CCDC88A, BLOC1S1, NCF4, LEAP2, and CDH17 showed increased expression in vitiligo lesions, whereas BACH2, ANKRD11, CD160 and TYRO3 showed reduced levels (Fig. 5D). Additionally, some genes showed cell type-specific enrichment in vitiligo lesions. BACH2, CD160, and CDH17 were specifically enriched in T & NK cells, LEAP2 showed increased expression in langerhans cells, and TYRO3 was specifically enriched in melanocytes (Fig. 5E). CXCL6 and CDH17 were not detected in normal skin tissue (Supplementary Fig. 5A).

We further investigated the heterogeneity of gene expression in T & NK subgroups. We extracted T & NK from the scRNA-seq dataset and identified six subgroups with distinct gene expression profiles through clustering and annotation, including five T cell clusters and one NK cell cluster. T cell clusters were divided into CD4 + Teff, CD4 + Treg, CD8 + Tem, CD8 + Trm, and CD8 + Tex (Fig. 6A). Compared to healthy controls, the percentages of CD8 + Tem, CD8 + Trm, and CD8 + Tex increased in vitiligo, while CD4 + Teff decreased (Fig. 6B). Each T cell subgroup exhibited enrichment in different biological processes (Supplementary Fig. 5C). We found that BACH2 and CDH17 were enriched in CD8 + Tex and upregulated in vitiligo. CD160 was enriched in NK cells and downregulated in vitiligo (Fig. 6C and Supplementary Fig. 5, D and E). Pseudotime analysis revealed CD8 + Tex and NK in the terminal stages of T cell differentiation (Supplementary Fig. 5F), indicating that BACH2, CD160, and CDH17 expression levels increased during T cell differentiation, with changes in CDH17 appearing in later stages (Fig. 6D).

Similarly, we performed further clustering and annotation of myeloid cells from the samples, identifying six subgroups: three mononuclear phagocyte subgroups (Mφ_CXCL8, Mφ_APOE, and Mono_FCN1) and three langerhans subgroups (LC1, LC2, and aLC) (Fig. 6E). The number of myeloid cell subtypes increased to varying degrees in vitiligo, particularly Mφ_APOE (Fig. 6F). Each myeloid cell subgroup showed enrichment in different biological processes (Supplementary Fig. 5H). We found that CCDC88A was highly expressed in Mφ_CXCL8 and Mono_FCN1, BLOC1S1 was highly expressed in Mono_FCN1, LC1, and LC2, and NCF4 was highly expressed in Mφ_APOE (Fig. 6G).

In this study, we conducted a meta-analysis of GWAS data from 589,083 vitiligo patients, identifying 25 distinct variants associated with vitiligo. By integrating methylation, gene expression, and protein abundance pQTL data, we identified eight putative genes related to immune response that may contribute to vitiligo susceptibility. Finally, using GWAS and MR analysis, we integrated candidate genes and validated the differential expression of nine genes in single-cell data from vitiligo lesions. Our study fills a gap in the genomic understanding of genes associated with vitiligo susceptibility and extends our knowledge of the biological pathways linked to all known vitiligo risk loci to date.

We employed various strategies to ensure the biological credibility of the 25 GWAS variants identified. Four of these variants showed genetic associations with risk factors for vitiligo. Additionally, we discovered that variants rs10676540, rs5845323, and rs9279765 were associated with autoimmune thyroiditis, and KEGG pathway analysis revealed enrichment in autoimmune thyroiditis among other autoimmune diseases. Vitiligo is commonly linked with various autoimmune diseases, with autoimmune thyroiditis being the most frequently occurring comorbidity^{[35], [36]}. Increased incidence of autoimmune thyroiditis has been reported in children with vitiligo^[37]. The gene HLA-DRB1, encoded by rs9279765, is involved in directing antigen-specific T-helper effector functions, including antibody-mediated immune responses and macrophage activation^{[38], [39], [40]}. Previous research has also indicated many common genetic susceptibility loci between generalized vitiligo and autoimmune thyroid diseases, with HLA-DR3 showing the strongest association^[41]. Our findings further validate the shared roles of both MHC and non-MHC candidate genes in the pathogenesis of these diseases. Moreover, numerous studies have indicated that vitiligo patients are susceptible to psychiatric disorders^{[42], [43]}. However, neither our association analysis nor LDSC analysis found genetic correlations between them. Likewise, Wang et al.'s MR analysis did not demonstrate a genetic association between generalized vitiligo and psychiatric disorders^[44]. These conclusions suggest that the correlation might be influenced by a complex interplay of multifactorial pathogenic mechanisms.

In this study, among the six newly identified significant vitiligo loci, five were successfully annotated. BACH2, CCDC88A, BLOC1S1, NCF4, and ANKRD11 were annotated as the priority genes. Several of these genes have already been extensively reported to be closely related to autoimmune conditions. The genetic polymorphisms in BACH2 are associated with many autoimmune and allergic diseases, including vitiligo^[45]. BACH2 inhibits the differentiation of CD4 + T cells into Th2 cells, suppresses the production of Th2 cytokines, promotes the differentiation of Tregs (regulatory T lymphocytes), and enhances Treg-mediated immunity^{[46], [47], [48]}. NCF4 mediates intracellular ROS levels and modulates autoreactivity and arthritogenic T cell activation in the regulation of autoimmunity and chronic inflammation^{[49], [50]}. Additionally, NCF4 mutations regulate the intrinsic oxidative burst in B cells, driving plasma cell formation and altering their CXCR3/CXCR4 expression^[51]. In this study, we also found that NCF4 is highly expressed in Mφ_APOE, which further supplements the evidence of NCF4 involving in the innate immune pathways in vitiligo. Based on these findings, we expanded the evidence and confirmed the association between BACH2 and NCF4 with the risk of vitiligo.

We conducted MR and SMR analyses using the plasma proteome pQTLs, and performed SMR validation at other QTL levels and in skin tissue. We found that eight immune response-related genes are associated with susceptibility to vitiligo, among which CD160 has a causal relationship with vitiligo at the levels of protein abundance, gene expression, and methylation. CD160 regulates the immune system and participates in autoimmunity, which is a critical gene abnormally expressed in autoimmune diseases^[52]. Our MR study has refined the potential role of CD160 in vitiligo and expanded the regulation of CD160 expression levels by several CpGs such as cg06095777.

We also found that CD160 is specifically enriched in T and NK cells in vitiligo patients, and its expression is significantly reduced in vitiligo lesions. CD160 plays contradictory roles in regulating T cells. On one hand, the surface protein CD160 is essential for the formation of memory in CD8 + T cells^[53]. Interaction of CD160 with MHC I enhances the function of CD8 + cytotoxic T lymphocytes (CTLs)^{[54], [55]}. In some studies, on the other hand, CD160 also negatively regulates CD4 + T cells and NKT cells^{[56], [57], [58]}. Additionally, CD160 activates NK cell toxicity through the PI3K/Akt/mTORC1 signaling pathway^{[59], [60]} and enhances the production of cytokines including IFN-γ, TNF-α, and IL-6^{[60], [61], [62]}. We speculate that CD160 predominantly exerts a T cell-suppressing function in vitiligo. However, more basic researches are needed to elucidate the effects of CD160 on the development and progression of vitiligo.

Our study boasts several strengths. Firstly, we acquired the largest sample size of vitiligo whole-genome data through GWAS meta-analysis, enabling us to validate previously reported vitiligo variants and discover new ones under these conditions. We utilized both Mendelian Randomization (MR) and Summary Mendelian Randomization (SMR) methods to estimate the causal effects of immune response-related genes using genetic variants. Additionally, we further validated our candidate genes by conducting MR analysis on vitiligo subtypes. Furthermore, by integrating evidence from multiple omics levels, we fortified the causal relationship between immune response-related genes and vitiligo risk factors. We also provided insights into the potential pathogenic roles of candidate genes in vitiligo through single-cell analysis.

However, there are still some limitations in our study. The current analysis is limited to European populations, and whether these findings are applicable to other ancestries needs further validation. Due to the limited number of immune response-related proteins in the pQTL dataset, and our subsequent analysis was based on candidate proteins identified through proteomic analysis, there might be some omissions in the candidate immune response-related genes obtained in this study. Additionally, although single-cell analysis identified the cellular type enrichment of most candidate genes in the disease, we did not replicate the results in tissue-specific eQTL analysis, which may also fail to reflect the actual role of these markers in the development of vitiligo.

In conclusion, we have discovered novel vitiligo-associated variants and several vitiligo-pathogenic genes through GWAS and multi-omics MR analysis, some of which have biologically plausible evidence. These findings of new potential causal relationships warrant further experimental and clinical research to provide new insights into the etiology of vitiligo and to offer new targets for the development of vitiligo biomarkers and therapeutic drugs.

Acknowledgments

We are grateful to the physicians and nurses of the Department of Dermatology at Jiangsu Provincial People's Hospital and The First Affiliated Hospital of Nanjing Medical University for their support of this study.

Author contributions

All authors contributed to the study conception and design. All authors read and approved the final manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (No. 82273549) and the Natural Science Foundation of Jiangsu Province (BK20221414).

Data Availability

The data we used are publicly available summary statistics and can be obtained upon reasonable request from the corresponding author.

Ethics approval and consent to participate

This study was conducted using published studies that provided publicly available summary statistics. All original studies were approved by the appropriate ethical review boards and participants provided informed consent. Therefore, approval from a new ethics review board was not required.

Consent for publication

All authors approved the submitted version.

Additional Information

Competing interests

The authors declare that they have no competing interests.

Bergqvist, C., Ezzedine, K. & Vitiligo A focus on pathogenesis and its therapeutic implications. J. Dermatol. 48, 252–270. 10.1111/1346-8138.15743 (2021).
Frisoli, M. L., Essien, K., Harris, J. E. & Vitiligo Mechanisms of Pathogenesis and Treatment. Annu. Rev. Immunol. 38, 621–648. 10.1146/annurev-immunol-100919-023531 (2020).
Luo, L., Zhu, J., Guo, Y. & Li, C. Mitophagy and immune infiltration in vitiligo: evidence from bioinformatics analysis. Front. Immunol. 14 10.3389/fimmu.2023.1164124 (2023).
Xu, Z. et al. Anatomically distinct fibroblast subsets determine skin autoimmune patterns. Nature. 601, 118–124. 10.1038/s41586-021-04221-8 (2022).
Gellatly, K. J. et al. scRNA-seq of human vitiligo reveals complex networks of subclinical immune activation and a role for CCR5 in T(reg) function. Sci. Transl. Med. 13, eabd8995. 10.1126/scitranslmed.abd8995 (2021).
Boef, A. G., Dekkers, O. M. & le Cessie, S. Mendelian randomization studies: a review of the approaches used and the quality of reporting. Int. J. Epidemiol. 44, 496–511. 10.1093/ije/dyv071 (2015).
Pagoni, P., Dimou, N. L., Murphy, N. & Stergiakouli, E. Using Mendelian randomisation to assess causality in observational studies. Evid. Based Ment. Health. 22, 67–71. 10.1136/ebmental-2019-300085 (2019).
Xu, S. et al. Oxidative stress gene expression, DNA methylation, and gut microbiota interaction trigger Crohn's disease: a multi-omics Mendelian randomization study. BMC Med. 21, 179. 10.1186/s12916-023-02878-8 (2023).
Chen, J. et al. Multi-omic insight into the molecular networks of mitochondrial dysfunction in the pathogenesis of inflammatory bowel disease. EBioMedicine 99, 104934, doi: (2024). 10.1016/j.ebiom.2023.104934
Jin, Y. et al. Genome-wide association studies of autoimmune vitiligo identify 23 new risk loci and highlight key pathways and regulatory variants. Nat. Genet. 48, 1418–1424. 10.1038/ng.3680 (2016).
Dang, Y. & Zhang, S. Causal relationship between endometriosis and inflammatory bowel disease: A Mendelian randomization analyses. Clin. translational Med. 14, e1496. 10.1002/ctm2.1496 (2024).
Bycroft, C. et al. The UK Biobank resource with deep phenotyping and genomic data. Nature. 562, 203–209. 10.1038/s41586-018-0579-z (2018).
Jin, Y. et al. Early-onset autoimmune vitiligo associated with an enhancer variant haplotype that upregulates class II HLA expression. Nat. Commun. 10, 391. 10.1038/s41467-019-08337-4 (2019).
Pietzner, M. et al. Mapping the proteo-genomic convergence of human diseases. Sci. (New York N Y). 374, eabj1541. 10.1126/science.abj1541 (2021).
Ferkingstad, E. et al. Large-scale integration of the plasma proteome with genetics and disease. Nat. Genet. 53, 1712–1721. 10.1038/s41588-021-00978-w (2021).
Sun, B. B. et al. Genomic atlas of the human plasma proteome. Nature. 558, 73–79. 10.1038/s41586-018-0175-2 (2018).
Suhre, K. et al. Connecting genetic risk to disease end points through the human blood plasma proteome. Nat. Commun. 8, 14357. 10.1038/ncomms14357 (2017).
Folkersen, L. et al. Genomic and drug target evaluation of 90 cardiovascular proteins in 30,931 individuals. Nat. metabolism. 2, 1135–1148. 10.1038/s42255-020-00287-2 (2020).
Yao, C. et al. Genome-wide mapping of plasma protein QTLs identifies putatively causal genes and pathways for cardiovascular disease. Nat. Commun. 9, 3268. 10.1038/s41467-018-05512-x (2018).
Suhre, K. Genetic associations with ratios between protein levels detect new pQTLs and reveal protein-protein interactions. Cell. genomics. 4, 100506. 10.1016/j.xgen.2024.100506 (2024).
Papadimitriou, N. et al. Physical activity and risks of breast and colorectal cancer: a Mendelian randomisation analysis. Nat. Commun. 11, 597. 10.1038/s41467-020-14389-8 (2020).
Võsa, U. et al. Large-scale cis- and trans-eQTL analyses identify thousands of genetic loci and polygenic scores that regulate blood gene expression. Nat. Genet. 53, 1300–1310. 10.1038/s41588-021-00913-z (2021).
McRae, A. F. et al. Identification of 55,000 Replicated DNA Methylation QTL. Sci. Rep. 8, 17605. 10.1038/s41598-018-35871-w (2018).
Consortium, G. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Sci. (New York N Y). 369, 1318–1330. 10.1126/science.aaz1776 (2020).
Willer, C. J., Li, Y. & Abecasis, G. R. METAL: fast and efficient meta-analysis of genomewide association scans. Bioinf. (Oxford England). 26, 2190–2191. 10.1093/bioinformatics/btq340 (2010).
Bulik-Sullivan, B. K. et al. LD Score regression distinguishes confounding from polygenicity in genome-wide association studies. Nat. Genet. 47, 291–295. 10.1038/ng.3211 (2015).
Bulik-Sullivan, B. et al. An atlas of genetic correlations across human diseases and traits. Nat. Genet. 47, 1236–1241. 10.1038/ng.3406 (2015).
Wielscher, M. et al. Genetic correlation and causal relationships between cardio-metabolic traits and lung function impairment. Genome Med. 13, 104. 10.1186/s13073-021-00914-x (2021).
Hemani, G. et al. The MR-Base platform supports systematic causal inference across the human phenome. eLife 7, doi: (2018). 10.7554/eLife.34408
Bowden, J., Davey Smith, G., Haycock, P. C. & Burgess, S. Consistent Estimation in Mendelian Randomization with Some Invalid Instruments Using a Weighted Median Estimator. Genet. Epidemiol. 40, 304–314. 10.1002/gepi.21965 (2016).
Wu, Y. et al. Integrative analysis of omics summary data reveals putative mechanisms underlying complex traits. Nat. Commun. 9, 918. 10.1038/s41467-018-03371-0 (2018).
Stuart, T. et al. Comprehensive Integration of Single-Cell Data. Cell. 177, 1888–1902e1821. 10.1016/j.cell.2019.05.031 (2019).
Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods. 16, 1289–1296. 10.1038/s41592-019-0619-0 (2019).
Trapnell, C. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat. Biotechnol. 32, 381–386. 10.1038/nbt.2859 (2014).
Rios-Duarte, J. A., Sanchez-Zapata, M. J. & Silverberg, J. I. Association of vitiligo with multiple cutaneous and extra-cutaneous autoimmune diseases: a nationwide cross-sectional study. Arch. Dermatol. Res. 315, 2597–2603. 10.1007/s00403-023-02661-y (2023).
Gill, L. et al. Comorbid autoimmune diseases in patients with vitiligo: A cross-sectional study. J. Am. Acad. Dermatol. 74, 295–302. 10.1016/j.jaad.2015.08.063 (2016).
Uncu, S., Yaylı, S., Bahadır, S., Okten, A. & Alpay, K. Relevance of autoimmune thyroiditis in children and adolescents with vitiligo. Int. J. Dermatol. 50, 175–179. 10.1111/j.1365-4632.2010.04665.x (2011).
Galperin, M. et al. CD4(+) T cell-mediated HLA class II cross-restriction in HIV controllers. Sci. Immunol. 3 10.1126/sciimmunol.aat0687 (2018).
Kwok, W. W. et al. Frequency of epitope-specific naive CD4(+) T cells correlates with immunodominance in the human memory repertoire. Journal of immunology (Baltimore, Md.: 188, 2537–2544, doi: (1950). 10.4049/jimmunol.1102190 (2012).
Muehling, L. M. et al. Circulating Memory CD4 + T Cells Target Conserved Epitopes of Rhinovirus Capsid Proteins and Respond Rapidly to Experimental Infection in Humans. Journal of immunology (Baltimore, Md.: 197, 3214–3224, doi: (1950). 10.4049/jimmunol.1600663 (2016).
Jacobson, E. M., Huber, A. & Tomer, Y. The HLA gene complex in thyroid autoimmunity: from epidemiology to etiology. J. Autoimmun. 30, 58–62. 10.1016/j.jaut.2007.11.010 (2008).
Kussainova, A. et al. Associations between serum levels of brain-derived neurotrophic factor, corticotropin releasing hormone and mental distress in vitiligo patients. Sci. Rep. 12, 7260. 10.1038/s41598-022-11028-8 (2022).
Di Bartolomeo, L. et al. Vitiligo and Mental Health: Natural Compounds' Usefulness. Antioxid. (Basel Switzerland). 12 10.3390/antiox12010176 (2023).
Wang, Y. et al. Exploring genetic associations between vitiligo and mental disorders using Mendelian randomization. Exp. Dermatol. 33, e14979. 10.1111/exd.14979 (2024).
Mazzieri, A., Montanucci, P., Basta, G. & Calafiore, R. The role behind the scenes of Tregs and Th17s in Hashimoto's thyroiditis: Toward a pivotal role of FOXP3 and BACH2. Front. Immunol. 13, 1098243. 10.3389/fimmu.2022.1098243 (2022).
Afzali, B. et al. BACH2 immunodeficiency illustrates an association between super-enhancers and haploinsufficiency. Nat. Immunol. 18, 813–823. 10.1038/ni.3753 (2017).
Liu, G. & Liu, F. Bach2: A Key Regulator in Th2-Related Immune Cells and Th2 Immune Response. Journal of immunology research 2814510, doi: (2022). 10.1155/2022/2814510 (2022).
Zhang, H. et al. Bach2 Deficiency Leads to Spontaneous Expansion of IL-4-Producing T Follicular Helper Cells and Autoimmunity. Frontiers in immunology 10, doi: (2050). 10.3389/fimmu.2019.02050 (2019).
Winter, S., Hultqvist Hopkins, M., Laulund, F. & Holmdahl, R. A. Reduction in Intracellular Reactive Oxygen Species Due to a Mutation in NCF4 Promotes Autoimmune Arthritis in Mice. Antioxid. Redox. Signal. 25, 983–996. 10.1089/ars.2016.6675 (2016).
Xu, J. et al. NCF4 regulates antigen presentation of cysteine peptides by intracellular oxidative response and restricts activation of autoreactive and arthritogenic T cells. Redox Biol. 72, 103132. 10.1016/j.redox.2024.103132 (2024).
He, C. et al. NCF4 dependent intracellular reactive oxygen species regulate plasma cell formation. Redox Biol. 56, 102422. 10.1016/j.redox.2022.102422 (2022).
He, W. et al. Aberrant Expressions of Co-stimulatory and Co-inhibitory Molecules in Autoimmune Diseases. Front. Immunol. 10 10.3389/fimmu.2019.00261 (2019).
Zhang, L. et al. CD160 Signaling Is Essential for CD8 + T Cell Memory Formation via Upregulation of 4-1BB. Journal of immunology (Baltimore, Md.: 211, 1367–1375, doi: (1950). 10.4049/jimmunol.2200792 (2023).
Agrawal, S. et al. Cutting edge: MHC class I triggering by a novel cell surface ligand costimulates proliferation of activated human T cells. Journal of immunology (Baltimore, Md.: 162, 1223–1226 (1999). (1950).
Tsujimura, K. et al. Characterization of murine CD160 + CD8 + T lymphocytes. Immunol. Lett. 106, 48–56. 10.1016/j.imlet.2006.04.006 (2006).
Cai, G. et al. CD160 inhibits activation of human CD4 + T cells through interaction with herpesvirus entry mediator. Nat. Immunol. 9, 176–185. 10.1038/ni1554 (2008).
El-Far, M. et al. CD160 isoforms and regulation of CD4 and CD8 T-cell responses. J. translational Med. 12, 217. 10.1186/s12967-014-0217-y (2014).
Kaye, J. CD160 and BTLA: LIGHTs out for CD4 + T cells. Nat. Immunol. 9, 122–124. 10.1038/ni0208-122 (2008).
Le Bouteiller, P. et al. CD160: a unique activating NK cell receptor. Immunol. Lett. 138, 93–96. 10.1016/j.imlet.2011.02.003 (2011).
Rabot, M. et al. CD160-activating NK cell effector functions depend on the phosphatidylinositol 3-kinase recruitment. Int. Immunol. 19, 401–409. 10.1093/intimm/dxm005 (2007).
Barakonyi, A. et al. Cutting edge: engagement of CD160 by its HLA-C physiological ligand triggers a unique cytokine profile secretion in the cytotoxic peripheral blood NK cell subset. Journal of immunology (Baltimore, Md.: 173, 5349–5354, doi: (1950). 10.4049/jimmunol.173.9.5349 (2004).
Tu, T. C. et al. CD160 is essential for NK-mediated IFN-γ production. J. Exp. Med. 212, 415–429. 10.1084/jem.20131601 (2015).

No competing interests reported.

Download PDF

Reviewers invited by journal
14 Sep, 2024
Editor assigned by journal
14 Sep, 2024
Editor invited by journal
05 Sep, 2024
Submission checks completed at journal
05 Sep, 2024
First submitted to journal
31 Aug, 2024

You are reading this latest preprint version

Genome-wide association analysis and multi-omic Mendelian randomization study exploring the immune response in vitiligo

Status:

Version 1

Abstract

Figures

Background

Methods

Study design

Data sources of vitiligo data

Datasets of protein, expression, and methylation quantitative trait loci

Meta-analysis of vitiligo GWAS

Annotation of risk variants

Pathway enrichment analysis

Associations of vitiligo risk variants with risk factors for vitiligo

LD score regression

MR analysis

SMR analysis

Single-cell transcriptomic expression analysis

Results

Genome-wide meta-analysis identifies significant vitiligo loci

Gene mapping of significant vitiligo genetic variants

MR and SMR identify immune response genes for vitiligo

Integrating multi-omics and tissue-specific evidence

Phenome-wide association analysis

Cell type-specific expression in skin tissue

Discussion

Conclusion

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1