Biological insights into systemic lupus erythematosus through an immune cell-specific transcriptome-wide association study.

OBJECTIVE: Genome-wide association studies (GWAS) have identified >100 risk loci for systemic lupus erythematosus (SLE), but the disease genes at most loci remain unclear, hampering translation of these genetic discoveries. We aimed to prioritise genes underlying the 110 SLE loci that were identified in the latest East Asian GWAS meta-analysis.
METHODS: We built gene expression predictive models in blood B cells, CD4+ and CD8+ T cells, monocytes, natural killer cells and peripheral blood cells of 105 Japanese individuals. We performed a transcriptome-wide association study (TWAS) using data from the latest genome-wide association meta-analysis of 208 370 East Asians and searched for candidate genes using TWAS and three data-driven computational approaches.
RESULTS: TWAS identified 171 genes for SLE (p<1.0×10-5); 114 (66.7%) showed significance only in a single cell type; 127 (74.3%) were in SLE GWAS loci. TWAS identified a strong association between CD83 and SLE (p<7.7×10-8). Meta-analysis of genetic associations in the existing 208 370 East Asian and additional 1498 cases and 3330 controls found a novel single-variant association at rs72836542 (OR=1.11, p=4.5×10-9) around CD83. For the 110 SLE loci, we identified 276 gene candidates, including 104 genes at recently-identified SLE novel loci. We demonstrated in vitro that putative causal variant rs61759532 exhibited an allele-specific regulatory effect on ACAP1, and that presence of the SLE risk allele decreased ACAP1 expression.
CONCLUSIONS: Cell-level TWAS in six types of immune cells complemented SLE gene discovery and guided the identification of novel genetic associations. The gene findings shed biological insights into SLE genetic associations. © Author(s) (or their employer(s)) 2022. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ.

Genome-wide association studies have identified >100 risk loci for systemic lupus erythematosus (SLE), but the disease genes at most of these loci remain unclear.

We built six immune cell-specific gene expression reference panels based on data from East Asians and performed a transcriptome-wide association study for SLE for the first time.

We identified 276 candidate disease genes in 110 SLE loci including 104 genes in novel loci.

We explored allele-specific regulatory mechanisms at ACAP1 that increase SLE risk.

We identified numerous potential drug targets for SLE.

Introduction

Systemic lupus erythematosus (SLE) is an autoimmune disorder with genetic predisposition.1 2 Genome-wide association studies (GWAS) have identified >100 genetic loci robustly associated with SLE.3–8 The latest effort was a genome-wide association meta-analysis in 208 370 East Asians (hereafter, East Asian meta-analysis), which identified 113 SLE loci.9 These GWAS have identified SLE genetic determinants and improved our understanding of disease pathogenesis. However, the disease genes through which genetic associations affect SLE remain unclear at most GWAS loci, hindering translation of genetic discoveries into SLE precision health.10 11

Experiments have been devoted to identify disease genes at SLE loci,12 13 but unsurprisingly are time consuming and expensive. Transcriptome-wide association study (TWAS) is an alternative statistical method for identifying candidate genes at GWAS loci.14 A TWAS usually consists of two steps. In step 1, TWAS learns gene expression predictive models in cohorts with both gene expression and genotype data (hereafter, gene expression reference). In step 2, TWAS uses the predictive models to impute in silico gene expression in cohorts only with genotype or GWAS summary statistics. After that, TWAS tests for associations between imputed gene expression and GWAS traits. TWAS gene expression predictive models use expression quantitative trait loci (eQTL) as predictors. Due to linkage disequilibrium (LD), TWAS can easily implicate hitchhiking genes together with disease genes in GWAS loci. However, methods have been developed to discern disease genes from hitchhikers.15 TWAS has recently nominated candidate genes for many human diseases,16 17 providing biological insights into disease associations. However, no TWAS for SLE has been reported.

TWAS requires gene expression references from populations with the same ancestry as in disease GWAS.18 Most current gene expression references are built from eQTL data sets of European ancestry. Only a few non-European references have been reported,19 limiting the application of TWAS in non-European samples.

Current TWAS mainly use disease-relevant tissue-level expression references that contain composite expression data from multiple distinct cell types in various cellular states. The heterogeneity in expression reference panels can bias TWAS findings and complicate gene association interpretation. In contrast, cell-level expression reference panels have obvious advantages, but few are available. Immune cells, such as B cells, T cells and monocytes, play key roles in SLE pathogenesis.20 TWAS using immune cell-level gene expression references might provide a unique opportunity to further our understanding of SLE pathogenesis.

To that end, we created gene expression predictive models from six types of blood immune cells obtained from East Asians. We performed a TWAS using the latest East Asian meta-analysis findings9 and searched for SLE genes jointly with three other data-driven gene prioritisation approaches. We identified 276 candidate genes, including 104 from recently-identified novel loci. We found that the six cell-level gene expression references complemented SLE gene discovery. We demonstrated that TWAS findings guide the identification of novel genetic associations. Additionally, we explored regulatory mechanisms at ACAP1 in vitro. Our findings provide biological insight into SLE pathogenesis.

Methods

Genome-wide association summary statistics

We previously performed the largest genome-wide association meta-analysis of SLE using data from 208 370 individuals of eight East Asian cohorts and identified 113 loci (including 46 novel loci) at p<5.0×10–8 (online supplemental table 1).9 In the present study, we used the index variants for the 110 autosomal loci and genome-wide single-variant association summary statistics at 11 270 530 genetic markers that were available in at least two member cohorts. We excluded the Human Leukocyte Antigens (HLA) region in further analyses. This study was carried out in compliance with the Helsinki Declaration.

TWAS and Fine-mapping Of CaUsal gene Sets

To infer gene expression changes in SLE, we performed a TWAS in FUSION21 using default parameters and eQTL data sets for six blood immune cell types generated from 105 (21 men and 84 women) healthy Japanese individuals with a mean age of 39 years: B cells, CD4+ T cells, CD8+ T cells, monocytes, natural killer (NK) cells and peripheral blood cells.22 LD was computed from whole-genome sequencing data of 3256 Japanese and 504 East Asians enrolled in the 1000 Genomes Project (1KGP).23 24 We restricted analysis to protein-coding genes. We defined a significant gene association p value threshold after Bonferroni correction for the number of protein-coding genes tested in each cell type (number of genes: B cells: 5055; CD4+ T cells: 5132; CD8+ T cells: 4988; monocytes: 5546; NK cells: 5239; peripheral blood cells: 5614).

To prioritise genes in genomic regions with ≥2 significant genes in TWAS, we implemented a Bayesian fine-mapping analysis using Fine-mapping Of CaUsal gene Sets (FOCUS) as previously described,15 which computed a posterior inclusion probability (PIP) for each gene to quantify the probability for being the true disease gene, and then created a 90% credible gene set that contained the putative disease genes with a probability ≥90%. We estimated gene expression weights in the six eQTL data sets and performed FOCUS for each cell type separately. We regarded TWAS significant genes with PIP ≥0.8 as potential disease genes.

Colocalisation analysis

To evaluate whether SLE GWAS associations share the same causal variants with eQTL, we performed colocalisation analysis between SLE GWAS loci and eQTL for genes with significant TWAS associations around the corresponding SLE GWAS index variants. We created ±100 kilobase (kb) genomic regions centring on SLE GWAS index variants and extracted association summary statistics from SLE East Asian meta-analysis9 and the corresponding eQTL.22 We restricted analysis to genetic variants with sample size N=208 370 in SLE GWAS meta-analysis. We implemented colocalisation analysis in coloc using default parameters.25 We defined significant colocalisation with posterior probability (PPH4) ≥0.8.

Data-driven Expression-Prioritised Integration for Complex Traits

To prioritise SLE genes, we analysed genetic variants with an SLE association p<5.0×10–8 using Data-driven Expression-Prioritised Integration for Complex Traits (DEPICT) V.1 release 194.26 DEPICT clumped the input variants in a 500 kb region at LD r2 >0.1 based on 1KGP East Asian data, yielding 1521 autosomal loci. We identified significant genes using a false discovery rate <5%.

Polygenic Priority Score

To prioritise SLE genes, we applied the Polygenic Priority Score (PoPS) method to East Asian meta-analysis results.27 First, we computed gene-level association statistics and gene–gene correlations from GWAS summary statistics using MAGMA28 and LD estimated from the 1KGP East Asian data. Next, we ran enrichment analysis for gene features listed at https://github.com/FinucaneLab/gene_features using MAGMA. We retained features with p<0.05 in MAGMA. Finally, we computed PoPS for each gene by fitting a joint model for enrichment of all resulting features. After calculating PoPS for a total of 18 383 protein-coding genes, we kept the top 30% of genes and prioritised those with the highest PoPS in a 1 megabase (Mb) window centred on each of the 110 SLE index variants.

Assay for transposase-accessible chromatin using sequencing in blood CD4+ T and CD19+ B cells

To detect open accessible elements of ACAP1, we sorted blood CD4+ T and CD19+ B cells from five healthy Chinese individuals and performed assay for transposase-accessible chromatin using sequencing (ATAC-seq) on the BGISEQ 500 platform as described previously.29 Each participant provided written informed consent.

Luciferase reporter assay

We previously identified rs61759532 as a likely causal variant at the ACAP1 locus.9 To explore the regulatory effects of rs61759532, three identical copies of the 24 bp-element flanking each allele of rs61759532 were subcloned into the luciferase vector, pGL4.26 (luc2/minP/Hydro), between the XhoI and BglII sites upstream of the minimal promoter for the firefly luciferase gene (online supplemental figure 1). The firefly luciferase vector (1 µg) and the normalising Renilla luciferase vector (500 ng) were co-transfected into human leukemia monocytic (THP1) cells for 2 days using Lipofectamine 3000 (Thermo Fisher Scientific). Luciferase activity was measured in five independent biological replicates using the Dual-Luciferase Reporter Assay Kit (Promega) according to the manufacturer’s instructions. Relative fold-change in firefly luciferase activity was normalised by both transfection efficiency, based on Renilla luciferase activity and minimal luciferase activity from the pGL4.26 vector without insert.

Electrophoretic mobility shift assay

Epstein-Barr Virus (EBV)-transformed B or THP1 cells were grown in RPMI 1640 medium including 10% fetal bovine serum and 1% penicillin/streptomycin. Electrophoretic mobility shift assay (EMSA) probes were constructed by annealing biotin-conjugated 30-residue oligonucleotide sequences flanking rs61759532. EMSA was performed using the LightShift Chemiluminescent EMSA Kit (Thermo Fisher Scientific) according to the manufacturer’s instructions.

Results

TWAS in the six immune cell types

To identify SLE genes, we performed a TWAS of East Asian meta-analysis data using gene expression references of six types of immune cells. We identified 57, 51, 48, 46, 44 and 40 significant genes in B, NK, peripheral blood, monocytes, CD4+ and CD8+ T cells, respectively (online supplemental table 2), which together comprised 171 genes. The significant genes were enriched in B cells (χ2 test; p=5.3×10–3). Colocalisation suggested that the same causal variants were shared between 24 SLE loci and 27 of the 171 genes (online supplemental table 3). Notably, only 5 of the 171 genes, namely B3GALT6, ELF1, HEATR3, TPCN2 and UHRF1BP1, attained significance in all six cell types; 114 (66.7%) genes showed significance in only a single cell type (figure 1A). We used PLD4, a newly-identified SLE gene,3 as a positive control. TWAS identified a significant association at PLD4 only in monocytes (p=2.6×10–9), suggesting that monocytes mediate the effects of PLD4 on SLE. A previous study reported that PLD4–/– mice developed more blood monocytes30 and autoimmune phenotypes.3

Figure 1

TWAS of East Asian meta-analysis data. (A) Distribution of significant genes across the six types of immune cells. (B) Number of significant TWAS genes per SLE locus. TWAS, transcriptome-wide association study.

Genes at SLE novel loci

Of the 171 genes, 127 (74.3%) arose within 500 kb from 61 of the 110 SLE index variants.9 For the majority (n=52; 85.2%) of these 61 SLE loci, TWAS identified ≤3 genes (figure 1B). For 33 loci, TWAS identified the closest protein-coding gene (online supplemental table 2).

Of the 127 genes, 35 came from 19 of the 46 recently-identified novel loci.9 For example, we identified a novel association at rs3750996 and prioritised rs3750996 as the putative causal variant for this association.9 TWAS consistently identified STIM1 gene near rs3750996 in B cells (p=1.1×10–7), CD4+ T cells (p=4.7×10–9), monocytes (p=5.8×10–9) and peripheral blood cells (p=3.38×10–8). Colocalisation suggested that the same causal variant was associated with both SLE risk and STIM1 expression in all the four types of cells (PPH4 >0.97; online supplemental table 3). STIM1 encodes a calcium channel sensor that regulates type I interferon response31 and plays an essential role in effector functions of T and B cells.32–34 Mutations in STIM1 cause severe immune deficiency in humans.35 STIM1 is a potential lupus therapeutic target.36

As another example, we previously identified a novel association at rs58107865 as an East Asian-specific SLE locus.9 TWAS identified at this locus the LEF1 gene (p=1.3×10–10), which encodes a transcription factor that binds to the T-cell receptor-α enhancer site. LEF1 controls the maintenance and functional specification of Treg subsets to prevent autoimmunity.37 LEF1 antagonist demonstrated tumour inhibition for B-cell chronic lymphocytic leukaemia,38 suggesting the potential of LEF1 as a drug target.

TWAS-guided identification of novel GWAS association

TWAS identified 44 genes (44=171–127) outside the 110 SLE loci, suggesting that future studies with larger sample sizes might detect novel GWAS associations with SLE around these 44 genes.

TWAS identified a significant association of CD83 with SLE in B cells (p=7.7×10–8). East Asian meta-analysis only found a borderline significant single-variant association (rs12530098 with the lowest p value; OR=1.10, p=6.9×10–8). We recruited two additional cohorts of 1498 SLE cases and 3330 controls in China39 and meta-analysed their summary-level associations with East Asian findings in this region. We identified a genome-wide significant association at rs72836542 around CD83 for the first time (OR=1.11, p=4.5×10–9; LD r2=0.93 with rs12530098; figure 2). SNP rs72836542 regulates CD83 expression in blood B cells (β=−1.51, p=4.2×10–20),22 suggesting that CD83 might mediate the association with SLE.

Figure 2

Locuszoom plot for a new single-variant association at the CD83 gene. Mb, megabase.

CD83 acts as an essential factor during the differentiation of T and B lymphocytes.40 Soluble human CD83-treated mice showed lower concentrations of anti-histone IgG autoantibodies and significantly delayed onset of anti-dsDNA autoantibody production.41 These reports suggest that CD83 is a promising drug target for SLE.41

Fine-mapping of TWAS genes

Of the 171 genes, 53 arose in TWAS of the same cell types. The flanking regions (±500 kb) for these 53 genes overlapped, suggesting they might arise at the same loci. These 53 genes comprised 17 genomic loci. To prioritise disease genes, we implemented a FOCUS analysis.15 We identified these 53 genes in the 90% gene credible sets and suggested 23 (43.4%) as likely disease genes (PIP ≥0.8; online supplemental table 4). Among them, FNIP1, HEATR3, and CD37 arose at SLE novel loci.

We reported a genome-wide association between SLE and rs11288784 for the first time in our East Asian meta-analysis.9 At this locus, TWAS identified three genes, ADCY7, BRD7 and HEATR3 (p<2×10–8; online supplemental table 2). Fine-mapping analysis suggested that HEATR3, the closest gene to the association, is most likely the disease gene (PIP >0.998). The eQTL for HEATR3 was colocalised with the SLE association (PPH4 >0.89; online supplemental table 3). HEATR3 plays a role in NOD2-mediated NF-κB signalling and has been implicated in Crohn’s disease.42

Complementary gene identification

To complement gene identification at the 110 SLE loci,9 we implemented three additional gene prioritisation approaches: (1) the nearest protein-coding gene; (2) DEPICT; and (3) PoPS. DEPICT and PoPS identified 54 and 107 protein-coding genes, respectively (online supplemental tables 5-7); 24 (44.4%) and 41 (38.3%) are the closest protein-coding genes to the corresponding SLE associations; 12 and 10 genes arose at SLE novel loci.9 TWAS and these three gene prioritisation approaches together identified 276 genes within the 110 SLE loci, including 104 genes at novel loci (online supplemental table 7). Notably, only seven genes (BANK1, IRF5, BLK, NCOA2, WDFY4, SLC15A4 and RASGRP1) were identified by all four methods, of which NCOA2 arises in the novel SLE locus at rs142937720.9 Colocalisation analysis using genetic associations with SLE susceptibility and NCOA2 expression revealed the sharing of causal variant (PPH4=0.93; online supplemental table 3). NCOA2 encodes a transcriptional co-activator of interferon regulatory factor 143 that plays a role in SLE.44

Regulatory mechanisms at ACAP1

One hundred and eighty-six (67.4%) of the 276 genes were identified in only one approach (figure 3). For example, DEPICT identified ACAP1 at the novel SLE locus around rs61759532.9 ACAP1 encodes a key regulator of integrin traffic for cell adhesion and migration.45 Fine-mapping analysis previously prioritised rs61759532, an intronic variant of ACAP1, as the likely causal variant (posterior probability of being causal =0.999; figure 4A). We found that rs61759532 overlaps with an accessible open chromatin region in blood B and T cells (figure 4B). GeneHancer46 suggested that rs61759532 resides in an enhancer/promoter element of ACAP1. Transcriptional reporter assays showed significant allelic differences in the enhancer activity of rs61759532 in THP1 monocyte cell lines (two-sided t-test p=8.1×10–3; figure 4C), consistent with the regulatory effect of the risk allele, T, in reducing ACAP1 expression in whole blood (p=1.7×10–47; figure 4D).47 EMSA revealed that allele-specific biotin-labelled probes containing T (risk allele) formed fewer nuclear protein-probe complexes than probes with C (non-risk allele) in THP1 and EBV-transformed B cell lines (figure 4E).

Figure 3

Venn diagram of candidate disease genes at the 110 SLE loci identified using four gene discovery approaches. DEPICT, Data-driven Expression-Prioritised Integration for Complex Traits; PoPS, Polygenic Priority Score; TWAS, transcriptome-wide association study.

Figure 4

Allele-specific regulatory effect of rs61759532 on ACAP1. (A) Regional association plot for the ACAP1 locus. The lead variant rs61759532 is labelled as a purple diamond. Linkage disequilibrium was estimated using data from 7021 Chinese individuals. (B) Location of rs61759532 within an assay for transposase-accessible chromatin using sequencing open chromatin accessible region in CD19+ B and CD4+ T cells (green tracks) and within active ChromHMM chromatin states (bars on the bottom panel) in primary CD8+ T naive cells (CD8.NPC), T helper naive cells (CD4.NPC) and primary B cells (BLD.CD19.PPC). Chromatin states are coloured red (active transcription start site), orange red (flanking active transcription start site), or yellow (enhancers). (C) Allelic differential enhancing activity of rs61759532 in THP1 cells. None, 3×C, and 3×T denote an empty vector containing a minimal promoter, and vectors with the C and T alleles of rs61759532, respectively. Relative luciferase activities, measured in five independent biological replicates, were significantly higher for inserts with the C allele (two-tailed t-test p=8.1×10–3). Error bars indicate SEMs of five independent biological replicates. (D) Association between the risk allele (T) of rs61759532 and decreased expression of ACAP1 in GTEx v8 whole blood (p=1.7×10–47). The white line in the centre of each box indicates the median expression value, while the box for each genotype represents the IQR of ACACP1 expression. (E) Allelic differential protein-DNA binding by rs61759532 in EMSAs. Biotin-conjugated 30-nucleotide probes flanking rs61759532 (denoted as C or T, according to the allele) were incubated with nuclear extracts (10 µg) from EBV-transformed B cells or THP1 cells in EMSAs. Shifted bands (indicated by red arrows) had stronger intensities with the biotin-conjugated C allele probes than the T allele probes and were not detected in the presence of excess non-conjugated probes. EBV: Epstein-Barr Virus; EMSA, electrophoretic mobility shift assay; Mb, megabase; THP1: human leukemia monocytic cell line.

Discussion

Here, we performed a TWAS for SLE for the first time and identified 171 genes associated with SLE risk. We nominated 276 genes at 110 SLE loci through TWAS and three computational approaches. One hundred and four genes arise at SLE novel loci; multiple show therapeutic potential. These findings provide insights into SLE biology and can guide future functional experiments.

SLE GWAS have identified >100 risk loci, but the disease genes and underlying molecular mechanisms remain largely unknown.3–8 TWAS is widely used to identify disease genes and determine disease mechanisms.14 In TWAS, population ancestry, tissue/cell relevance and cell sources of gene expression references are critical.48 49 Here, we created cell-level gene expression references from six types of immune cells in East Asians, ensuring that the reference panels were constructed from individuals with the same ancestry as the SLE GWAS. Various immune cells play a role in SLE pathogenesis.20 Studies suggested that loci identified in SLE GWAS could contribute to the risk of SLE through their effects on immune cells.50 Our study showed that 66.7% of these significant genes attained significance only in TWAS of one of the six immune cell types. This finding highlights the value of evaluating diverse immune cells in TWAS of SLE.

TWAS identified 44 genes associated with SLE in regions without prior GWAS associations. Among these 44 genes, we identified a genome-wide single-variant association at CD83 for the first time. CD83 modulates the production of autoantibodies and might have therapeutic effects in SLE.41 This result demonstrates that TWAS can help guide the identification of novel GWAS associations.

For the 110 SLE loci that we recently identified in our latest East Asian meta-analysis,9 TWAS and the three data-driven approaches identified a pool of 276 gene candidates, 186 of which were identified using a single approach. These gene findings warrant careful interpretation. We previously identified rs61759532 as a putative causal variant of SLE.9 In the present study, we demonstrated in vitro the molecular effects of the different alleles of rs61759532 on ACAP1 expression levels. We showed that rs61759532 resides in an open chromatin region and exhibited enhancing activity on ACAP1. The risk allele T of rs61759532 reduces the expression of ACAP1 in whole blood.

This study has several limitations. The modest study sample size in the cell-level gene expression references likely limited the power and precision of TWAS. SLE has various systemic manifestations, suggesting that many tissues/cells contribute to disease pathogenesis in addition to the immune cells that we studied.20 Increasing the breadth of cell types and cell state resources in gene expression references would increase the precision of TWAS. We only experimentally explored functional mechanisms for one significant SNP (rs61759532) in one gene, ACAP1. The role of ACAP1 and the biological pathways mediating the effects of ACAP1 on SLE are worthy of further investigation.

In summary, we performed a TWAS for SLE for the first time and identified 276 gene candidates at SLE loci. These findings help elucidate the genetic mechanisms underlying SLE and provide potential SLE therapeutic targets.