Genome-wide association study identifies Sjögren's risk loci with functional implications in immune and glandular cells

Abstract

Sjögren's disease is a complex autoimmune disease with twelve established susceptibility loci. This genome-wide association study (GWAS) identifies ten novel genome-wide significant (GWS) regions in Sjögren's cases of European ancestry: CD247, NAB1, PTTG1-MIR146A, PRDM1-ATG5, TNFAIP3, XKR6, MAPT-CRHR1, RPTOR-CHMP6-BAIAP6, TYK2, SYNGR1. Polygenic risk scores yield predictability (AUROC = 0.71) and relative risk of 12.08. Interrogation of bioinformatics databases refine the associations, define local regulatory networks of GWS SNPs from the 95% credible set, and expand the implicated gene list to >40. Many GWS SNPs are eQTLs for genes within topologically associated domains in immune cells and/or eQTLs in the main target tissue, salivary glands.

Fig. 1. Schematic of the Sjögren’s GWAS and Bioinformatic Workflow.

a Workflow of the six genotyped datasets (DS1–6) and one ImmunoChip dataset (DS7) used in this study, including the number of cases, controls, and SNPs included in each dataset pre- and post-quality control, and after whole-genome imputation. The post-imputation merged dataset (PI1) containing DS1–6 was used to perform the SNP-Sjögren’s single marker trait analysis (orange), the polygenic risk score (PRS) analysis (yellow), the genetic correlation analyses (blue), and the epigenetic enrichment analyses (blue). Meta-analysis was performed using the genotyped PI1 merged dataset and DS7 ImmunoChip dataset (green) merged using the DS7 genotyping platform. See Supplementary Data 1 for detailed information for each dataset. b Statistical and bioinformatic analysis workflow applied to each novel risk locus to identify and predict functionality of likely functional SNPs.

Fig. 2. Summary of the SNP-Sjögren associations in a European population.

a Manhattan plot shows the summary data from the meta-analysis of the 7.3 × 10⁵ SNPs shared between the GWAS and ImmunoChip dataset (Supplementary Data 1) after imputation. The −log₁₀(P) for each variant is plotted according to chromosome and base pair position. A total of seven novel loci (indicated in red) exceeded genome-wide significance (GWS) of P_GWAS < 5 × 10⁻⁸ (red dashed line). Three additional novel loci (indicated in royal blue) exceed GWS after meta-analysis (P_META < 5 × 10⁻⁸). Several previously established loci were replicated (indicated in light blue). The suggestive GWAS and meta-analysis threshold (P_Suggestive < 5 × 10⁻⁵) is indicated by the blue dashed line. b–h Logistical regression analysis was performed on GWAS dataset PI1 after imputation, identifying the top SNP associations (e.g., index SNPs) of the novel GWS regions of association: NAB1 (b), PTTG1-MIR146A (c), XKR6 (d), MAPT-CRHR1 (e), RPTOR-CHMP6-BAIAP6 (f), TYK2 (g), SYNGR1 (h). i–k Logistical regression analysis was performed after meta-analysis of the GWAS and ImmunoChip data, identifying the top SNP associations (e.g., index SNPs) of the novel GWS regions of association: CD247 (i), PRDM1-ATG5 (j), TNFAIP3 (k).

Fig. 3. Polygenic risk score analysis of the Sjögren-SNPs in all genotype Sjögren’s cases and Ro⁺ Sjögren’s cases with or without the HLA region.

a Heatmap of LDSC-estimated genetic correlations between Sjögren’s and 19 other immune-mediated diseases and other common traits using European GWAS summary data from the 1000 Genomes Project. Box color indicates magnitude of correlation; * indicates significant P-value after Bonferroni correction. b–m Polygenic risk scores (PRS) were calculated for all genotyped individuals from the Sjögren-All (b–g) or Sjögren-Ro⁺ (h–m) subsets, divided 2/3 into training and 1/3 into testing datasets, using LD-pruned genotyped SNPs (r² > 0.2) including (b–d, h–j) or excluding (e–g, k–m) SNPs from the HLA region (6p21.3–22.3.). (b, e, h, k) Bar plot of multiple P-value thresholds (P_T) for PRS prediction of Sjögren’s. c, f, i, l Histogram of the PRS distribution in Sjögren’s cases (orange) and controls (teal). d, g, j, m Strata plot of the odds ratio (OR) when comparing PRS from different quantile ranges. Bars indicate the 95% confidence intervals of the OR from each quantile range. n Area under the receiver operating characteristic curve (AUROC) values of PRS using LD-pruned genotyped Sjögren-SNPs including SNPs from 6p21.3–22.3 in all genotyped Sjögren’s cases (dark blue) or Ro⁺ Sjögren’s cases (green) relative to population controls, and LD-pruned genotyped Sjögren-SNPs excluding SNPs from 6p21.3–22.3 in all genotyped Sjögren’s cases (red) or Ro⁺ Sjögren’s cases (light blue) relative to population controls.

Fig. 4. Epigenetic and transcription factor enrichment analysis of Sjögren’s risk loci in specific cell types and tissues.

a Enrichment analysis of 30 histone marks in 127 different cell types and tissues from the Roadmap Epigenome Consortium Project was performed using GREGOR software. Heatmap displays the enrichment P-value for each histone mark plotted relative to specific immune cell types from the blood (red), spleen (yellow), and thymus (green). Black boxes indicate missing data. Complete analysis of all 127 cell types and tissues from the Roadmap Epigenome Consortium Project are reported in Supplemental Fig. 8 and Supplemental Data 3. b Enrichment of heritability in cell type-specific enhancers reported in the EnhancerAtlas2.0 database. Enrichment was calculated as partitioned heritability divided by the proportion of GWS Sjögren-SNPs that intersect with enhancer regions in each cell type. Cell types exhibiting significant GWS Sjögren-SNPs in enhancers are indicated in red. Blue dashed line is the threshold of significance after correction for multiple testing of P ≲ 0.05. c Enrichment of cell type-specific transcription factor binding sites. Transcription factors are indicated by a blue box if at least one Bonferroni-corrected (P ≲ 0.05) intersection between an indicated transcription factor and Sjögren-associated risk locus was detected. Red font indicates transcription factors that are also associated with the EBNA2 super enhancer. The cell type with the most significant interaction is listed in parentheses for each transcription factor.

Fig. 5. Chromatin Interactions and eQTLs of the ten novel Sjögren-associated genetic risk loci.

a Circos plot shows the zoom regional Manhattan plots for each genetic risk locus (outer most layer); SNPs with P-value <0.05 (black); r² > 0.08 (red); r² > 0.06 (orange). Index SNP rsIDs are indicated in red. Black rsIDs are prioritized SNPs from the 95% credible set that are also eQTLs that exhibit chromatin-chromatin interactions and are shown in (b). Outer circle displays the chromosome coordinate with the genomic risk loci highlighted in blue. Genes that are eQTLs (green) or exhibit chromatin interaction by Hi-C in Epstein–Barr virus (EBV)-transformed B lymphocytes (orange) are reported on the inner circles as text or interaction links. Each index gene is colored blue. Genes that are eQTLs and engage in chromatin interactions are reported in purple. b Cell type-specific functional annotations (horizontal rectangles), select eQTLs (top triangles), and chromatin-chromatin interactions (bottom triangle) are shown for the indicated prioritized SNPs from each 95% credible set. MIR146A was omitted because mined eQTL databases did not test MIR146A. Complex linkage disequilibrium of the CRHR1 association impaired refinement and fine-mapping of the region. c IMPACT annotation of the most likely functional Sjögren-SNPs to quantify SNP position in 700 cell-type-specific active transcription factor binding sites. Top panel depicts SNP position (blue lines) relative to genomic coordinates (Mb) of each indicated locus. Bottom panel shows the total number of active transcription factor binding sites detected at each indicated SNP.