We first collected all SLE-associated risk loci reaching genome-wide association significance (p < 5 × 10⁻⁸) published through March 2018 (Supplementary Data 1). Studies of all ancestral groups were included, and independent risk loci were defined as loci with lead (tag) variants at r² < 0.2. For each of these 91 risk loci, we performed linkage disequilibrium (LD) expansion (r² > 0.8) in each ancestry of the initial genetic association(s), to include all possible disease-relevant variants (Supplementary Data 2). In total, this procedure identified 3073 genetic variants. All established alleles of these variants were included, with 36 variants having three or more alleles. We also included 20 additional genetic variants from a previously published study¹⁹ as positive and negative controls to assess the library’s performance (Supplementary Data 3).

For each variant, we generated a pair of 170 base pair (bp) DNA oligonucleotides (subsequently referred to as “oligos”) for each allele, with the variant located in the center and identical flanking genomic sequence across the alleles (Supplementary Data 4). A total of 12,478 oligos (3093 variants with 6239 alleles) were synthesized. For barcoding, a random 20mers were added to each oligo through PCR. Each unique barcode was matched with perfectly synthesized oligos. The number of unique barcodes per oligo had an approximately normal distribution with a median of 729 barcodes per oligo (Supplementary Fig. 2a and Supplementary Data 5). Only oligos with at least 30 unique barcodes were used for downstream analyses. A fragment containing a minimal promoter and an eGFP gene was inserted between the oligo and barcode to generate the MPRA transfection library. We note that the use of a minimal promoter allows us to effectively measure the ability of alleles to enhance, but not reduce, transcriptional activity. Three aliquots of the library were independently transfected into the EBV-transformed B cell line GM12878. We then used nucleic acid capture to enrich for eGFP mRNA and sequenced the barcode region. The normalized barcode ratio between the eGFP mRNA and the plasmid DNA was used to quantify the amount of enhancer activity driven by each oligo (Supplementary Note 1 and Supplementary Fig. 5f, g). This mRNA to DNA ratio measures the enhancing effect of an allele on eGFP expression under the control of a minimal promoter (Fig. 1 and Supplementary Fig. 1). We observed strong correlation of enhancer activity between experimental replicates (mean pairwise Pearson correlation of 0.99) (Supplementary Fig. 2b-d). Likewise, calibration variants showed high accuracy, with 17 of the 20 variants matching the results of a previous study¹⁹ (87.5% sensitivity and 75% specificity), collectively demonstrating a robust experimental system (Supplementary Data 3).

Using the SLE MPRA library, we next identified genetic variants capable of driving enhancer activity in the EBV-transformed B cell line GM12878. An SLE-risk variant was considered a candidate for enhancer activity if an oligo corresponding to any allele had significantly increased transcriptional regulatory activity compared to controls (see “Methods”). Not all statistically significant changes in transcriptional activity are necessarily biologically relevant—a highly consistent, but slight change in expression levels is statistically, but not biologically, meaningful. We therefore considered an oligo to have enhancer activity only when (1) the oligo had statistically significant enhancer activity (p_adj < 0.05) and (2) we observed at least a 50% increase in transcriptional activity compared to the corresponding barcode counts in the plasmid control. Based on these criteria, 16% of SLE-risk variants (482 variants, 853 alleles) demonstrated enhancer activity, henceforth referred to as “enhancer variants” (enVars) and “enhancer alleles” (enAlleles), respectively (Fig. 2a and Supplementary Data 6).Fig. 2

a Distribution of MPRA regulatory activity. The normalized fold change of MPRA activity relative to plasmid control (X-axis) was calculated using DESeq2 (n = 3 biological replicates). Enhancer alleles (enAlleles) (blue) were identified as those alleles with significant activity relative to control (p_adj < 0.05) and at least a 50% increase in activity (see “Methods”). The p-values were generated by two-sided Wald tests with Benjamini–Hochberg multiple testing correction. Full results are provided in Supplementary Data 6. b Enrichment of histone marks in GM12878 cells at enVars compared to non-enVars. p-values were estimated by one-sided z-test with Bonferroni multiple testing correction using RELI (see “Methods”). Full results are provided in Supplementary Data 9. c Enrichment of regulatory protein and transcription factor (TF) binding at enVars compared to non-enVars. p-values were estimated by one-sided z-test with Bonferroni multiple testing correction using RELI (see “Methods”). The top 15 TFs (based on RELI p-values) that overlap at least 10% of enVars are shown. Full results are provided in Supplementary Data 9. d TF binding site motif enrichment for enVars compared to non-enVars. p-values were estimated by one-sided hypergeometric test with Benjamini–Hochberg multiple testing correction by HOMER using the full oligo sequences of enVars and non-enVars (see “Methods”). The top 15 enriched TF motif families are shown. Full results are provided in Supplementary Data 10.

We next explored the potential effects of enVars on gene expression. We connected each enVar to one or more genes using an approach that takes into account chromatin looping interactions, expression quantitative trait loci (eQTLs), and gene proximity (Supplementary Data 2 and 7) (see “Methods”). This approach identified 1006 genes in total, which are enriched for expected SLE-related processes such as the interferon pathway, the antigen processing and presentation pathway, and cytokine-related pathways (Supplementary Fig. 3 and Supplementary Data 8), providing functional support for the enVars we identified.

Next, we searched for functional genomic features enriched within enVars relative to non-enVars using the RELI algorithm¹⁰. In brief, RELI estimates the significance of the intersection between an input set of genomic regions (e.g., enVars) and each member of a collection of functional genomics datasets (e.g., ChIP-seq for a particular histone mark or TF). For this analysis, we identified, curated, and systematically processed the 576 GM12878 ChIP-seq datasets available in the NCBI Gene Expression Omnibus (GEO) database (see “Methods”). Using RELI, we observed significant enrichment for overlap between enVars and multiple histone modification marks, including H3K4me3 (5.8-fold, p_corrected < 10⁻²¹) and H3K27ac (2.0-fold, p_corrected < 10⁻¹³) (Fig. 2b and Supplementary Data 9). As expected, we did not identify enrichment for repressive marks such as H3K9me3 or H3K27me3²⁸ (Fig. 2b and Supplementary Data 9). Altogether, the genomic features present within enVars confirm that many SLE genetic risk loci likely alter transcriptional regulation in EBV-transformed B cells.

We next asked if the genomic binding sites of particular TFs were enriched within our enVars using RELI and the TF ChIP-seq datasets from GM12878. As expected, the enVars are highly enriched for ChIP-seq signal of TFs involved in regulation of the immune response, relative to variants lacking enhancer activity (Fig. 2c and Supplementary Data 9). In particular, we found significant enrichment for all members of the NFκB TF family: REL/C-Rel (6.4-fold, p_corrected < 10⁻²⁶), NFKB1/p50 (3.0-fold, p_corrected < 10⁻¹⁸), RELA/p65 (3.1-fold, p_corrected < 10⁻¹⁶), RELB (2.7-fold, p_corrected < 10⁻¹⁰), and NFKB2/p52 (2.2-fold, p_corrected < 10⁻⁷). These results are consistent with our previous findings that altered binding of NFκB TFs is likely an important mechanism conferring SLE risk¹⁰. We also found significant enrichment for other TFs that have been previously implicated in SLE pathogenesis, such as PAX5²⁹, MED1³⁰, IKZF1³¹, ELF1³², and the EBV-encoded EBNA2 transactivator¹⁰ (Fig. 2c and Supplementary Data 9). As a complementary approach, we next assessed enrichment for TF binding site motifs in the enAllele DNA sequences using HOMER³³ and motifs contained in the Cis-BP database³⁴ (see “Methods”). This analysis also revealed enrichment of multiple TF families with known roles in SLE, including ETS, NFκB, and IRF³ (Fig. 2d and Supplementary Data 10). Many of these same TFs also have enriched ChIP-seq peaks at SLE-risk loci¹⁰. Collectively, these results indicate that particular TFs tend to not only concentrate at SLE-risk loci¹⁰, but also concentrate at alleles capable of driving gene expression in EBV-transformed B cells.

We next used our MPRA library to identify SLE genetic risk variants that drive allele-dependent (allelic) enhancer activity. Allelic activity was assessed for each enVar by comparing enhancer activity between each pairs of alleles. We considered a SLE variant allelic if (1) at least one of its alleles is an enAllele; (2) we observed significant genotype-dependent activity using Student’s t-test^19,35 (Supplementary Note 2 and Supplementary Fig. 6); and (3) the oligos had more than a 25% change between any pair of alleles. Using these criteria, we identified 51 SLE-risk variants (11% of enVars, 1.7% of all SLE-risk variants) as allelic enVars in GM12878 (Fig. 3a and Supplementary Data 11). For 31 of these 51 allelic enVars, the risk allele decreased enhancer activity relative to the non-risk allele, which is statistically indistinguishable from the 20 variants with increased risk allele activity (p = 0.1). Three of the allelic enVars can also alter the amino acid sequence of proteins—rs1059702 (IRAK1), rs1804182 (PLAT), and rs3027878 (HCFC1), consistent with previous studies identifying dual-use codons in the human genome³⁶. Collectively, these 51 variants represent causal variant candidates for 27 SLE-risk loci (30% of all tested loci) (Supplementary Data 12). For these 27 risk loci, our approach reduced the number of potential causal variants with allelic activity in GM12878 from an average of 84 variants to an average of two variants per risk locus (Fig. 3b). For example, at 17q12 (marked by rs8079075), we reduce the candidate causal variant set from 249 to one, with the rs112569955 “G” risk allele showing a 36% increase in enhancer activity compared to the “A” non-risk allele.Fig. 3

a Identification of allelic enVars. Genotype dependence (Y-axis) is defined as the normalized fold change of MPRA activity between the non-reference and reference alleles (n = 3 biological replicates, see “Methods”). MPRA activity (X-axis) is presented as the maximum normalized fold change of MPRA activity for any allele of the variant. Allelic enVars (red) were defined as variants with a significant difference in MPRA activity (p_adj < 0.05) between any pair of alleles and at least a 25% change in activity difference (see “Methods”). The p-values were generated by two-sided Student’s t-test with Benjamini–Hochberg multiple testing correction. Full results are provided in Supplementary Data 11. b MPRA enhancer activity at the 27 risk loci with at least one allelic enVar. Bar plots indicate the total number of variants at each locus. Variants with allelic enhancer activity (allelic enVars) are shown in red. Variants lacking allelic enhancer activity are shown in gray.

To identify candidate regulatory proteins that might participate in allelic SLE mechanisms, we next used RELI to identify GM12878 ChIP-seq datasets that significantly overlap allelic enVars (Supplementary Data 13). Many of the top results are consistent with our previous study¹⁰, including the enriched presence of general enhancer features such as the H3K27ac histone mark (17 of 51 allelic enVars, 13.6-fold enriched, p_corrected < 10⁻³⁸), mediator complex subunit MED1 (17 of 51 allelic enVars, 13.0-fold enriched, p_corrected < 10⁻³⁴), and the histone acetyltransferase p300 (16 of 51 allelic enVars, 12.4-fold enriched, p_corrected < 10⁻³²), along with particular regulatory proteins that participate in “EBV super enhancers”³⁷ and play key roles in B cells such as ATF7 (15 of 51 allelic enVars, 11.3-fold enriched, p_corrected < 10⁻²⁶), Ikaros/IKZF1 (19 of 51 allelic enVars, 9.7-fold enriched, p_corrected < 10⁻²⁵), and the NFκB subunit RELA (13 of 51 allelic enVars, 12.4-fold enriched, p_corrected < 10⁻²⁴). Also consistent with our previous study¹⁰, we observe strong enrichment for the EBV-encoded EBNA2 protein (7 of 51 allelic enVars, 17.7-fold enriched, p_corrected < 10⁻¹⁹). Collectively, these data reveal particular regulatory proteins that might participate in the mechanisms contributing to SLE at multiple risk loci by driving allelic enhancer activity.

We next used the MARIO pipeline¹⁰ to search for allelic binding events (i.e., allelic imbalance between sequencing read counts) at SLE variants within 1058 LCL ChIP-seq datasets (576 from GM12878). By necessity, this approach is limited to the 47 allelic enVars that are heterozygous in at least one of these cell lines. In total, this procedure identified 11 variants with strong allelic imbalance (MARIO ARS value >0.4) in at least one ChIP-seq dataset (Supplementary Data 14), revealing groups of TFs and transcriptional regulators that allelically bind SLE-risk variants with genotype-dependent MPRA enhancer activity. For example, the rs3101018 variant, which is associated with SLE³⁸ and rheumatoid arthritis³⁹ in Europeans, shows 1.7-fold stronger enhancer activity for the reference/non-risk ‘C’ allele compared to the non-reference/risk ‘T’ allele (Fig. 4a). These results are consistent with a previously established eQTL obtained from GTEx⁴⁰, which demonstrates higher Complement C4A (C4A) expression in EBV-transformed B cell lines for the rs3101018 ‘C’ allele than the ‘T’ allele (Fig. 4b). Our MARIO allelic ChIP-seq analysis reveals 15 regulatory proteins that prefer the ‘C’ allele and 2 that prefer the ‘T’ allele (Fig. 4c and Supplementary Fig. 4a). Among these, particularly robust signal is obtained for ATF7, with one experimental replicate in GM12878 displaying 77 vs. 18 reads (‘C’ vs. ‘T’) and another showing 66 vs. 23 reads (‘C’ vs. ‘T’) (Supplementary Data 14). Moreover, CREB1 and CREM strongly favor the ‘C’ allele as well (Supplementary Data 14). In agreement with these data, computational analysis of the DNA sequences surrounding this variant predicts that ATF7, CREB1, and CREM will all bind more strongly to the ‘C’ than the ‘T’ allele (Fig. 4d). Intriguingly, ten additional proteins (FOXK2, PKNOX1, ARID3A, ZBTB40, ZNF217, ARNT, ELF1, IKZF2, MEF2B, and FOXM1) also bind allelically and have known DNA-binding motifs⁴¹, but none of them have binding sites altered by the variant. Further, we do not observe allelic chromatin accessibility in an available GM12878 ATAC-seq dataset (9 vs. 7 unique reads). Together, these results reveal a potentially causative SLE regulatory mechanism involving weaker direct binding of ATF7/CREB1/CREM to the ‘T’ risk allele, altering the recruitment of additional proteins to the locus and lowering the expression of C4A.Fig. 4

a, e Normalized MPRA enhancer activity of each experimental replicate for rs3101018 and rs26069235. b, f Expression trait quantitative loci (eQTLs) revealing genotype-dependent expression of C4A and SNYGR1 for rs3101018 (CC, n = 127 biologically independent samples; CT, n = 17; TT, n = 3) and rs26069235 (GG, n = 72 biologically independent samples; GA, n = 66; AA, n = 9) in EBV-transformed B cell lines (GTEx). c, g Genotype-dependent activity of transcription factors, transcriptional regulators, and histone marks in EBV-transformed B cell lines for rs3101018 and rs26069235. Results with MARIO ARS value >0.4 and consistent allelic imbalance across ChIP-seq datasets are included (see “Methods”). The X-axis indicates the preferred allele, along with a value indicating the strength of the allelic behavior, calculated as one minus the ratio of the weak to strong read counts (e.g., 0.5 indicates the strong allele has twice the reads of the weak allele). The median value is plotted when data from multiple cell lines are available, with full results provided in Supplementary Fig. 4. The numbers in parentheses represent the number of ChIP-seq datasets with significant allelic activity (i.e., MARIO ARS value >0.4) out of the number of datasets where the given variant is inside a ChIP-seq peak and is also heterozygous in the given cell line. Variant overlapping TFs are indicated in black. Variant adjacent TFs are shown in green (see definition in Fig. 5a). d, h DNA-binding motif logos are shown for the ATF/CREB/CREM family, and ELF1 in the context of the DNA sequence surrounding rs3101018 and rs2069235, respectively. Tall nucleotides above the X-axis indicate preferred DNA bases. Bases below the X-axis are disfavored. In (b) and (f), data are represented as a violin plot where the middle line is the median, the lower and upper hinges correspond to the first and third quartiles, with the rotated kernel density plot shown on each side. The data used for the analyses were obtained from the Genotype-Tissue Expression (GTEx) Portal on 11/12/2020. The GTEx Project was supported by the Common Fund of the Office of the Director of the National Institutes of Health, and by NCI, NHGRI, NHLBI, NIDA, NIMH, and NINDS.

We observe a similar phenomenon for the rs2069235 variant, which is associated with SLE in Asian ancestries⁴² and rheumatoid arthritis in Europeans⁴³. rs2069235 displays much stronger enhancer activity for the ‘A’ (non-reference/risk) allele compared to the ‘G’ (reference/non-risk) allele (Fig. 4e), consistent with the established synaptogyrin 1 (SYNGR1) eQTL in EBV-transformed B cell lines⁴⁰ (Fig. 4f). Inspection of our allelic ChIP-seq data reveals 14 proteins preferentially binding the ‘A’ allele, with none preferring the ‘G’ allele (Fig. 4g and Supplementary Fig. 4b). Among these 14 proteins, only ELF1 has its binding site directly altered by the variant (Fig. 4h). Strikingly, 55 of the 78 available H3K27ac datasets are allelic at this variant, with all 55 preferring the ‘A’ allele. Likewise, 24 of 46 H3K4me1 datasets are allelic, with all of them also preferring the ‘A’ allele. Both of these histone marks are indicative of active chromatin²⁸. Only a single histone mark dataset prefers the ‘G’ allele—the H3K27me3 mark, which is indicative of silenced chromatin²⁸ (Fig. 4g and Supplementary Fig. 4b). Together, these data are consistent with a potentially causative SLE molecular mechanism involving an allele-dependent enhancer consisting of stronger direct binding of ELF1 to the ‘A’ risk allele, along with indirectly altered binding of multiple additional TFs to this locus.

As illustrated by the above examples, a particular TF can be involved in allelic mechanisms that are either directly impacted by a given variant (i.e., the variant directly alters the TF’s DNA-binding site) or indirectly impacted by the variant (i.e., the variant alters the DNA binding of the TF’s physical interaction partner, modulates chromatin accessibility, or affects another mechanism). At a given locus, we designate such TFs as variant overlapping and variant adjacent TFs, respectively (Fig. 5a). We next sought to discover such TFs at the 51 allelic enVars. At each allelic enVar locus, we identified variant overlapping TFs as those TFs predicted to have strong binding to one allele and weak binding to the other allele. Likewise, we identified variant adjacent TFs as those TFs with proximal strong predicted binding sites that do not directly overlap the variant (see “Methods”). We then searched for particular TFs that tend to act as variant overlapping TFs or as variant adjacent TFs at the 51 allelic eVars using a proportion test (see “Methods”) and confirmed that their binding site locations are distributed relative to the variant as expected (Fig. 5b). Consistent with our results at the C4A and SYNGR1 loci, variant overlapping TFs include members of the ETS (e.g., ELF1) and ATF-like (e.g., ATF7) families, along with other TFs whose genetic loci are associated with SLE, including IRF5⁴⁴ (Fig. 5c and Supplementary Data 15). Variant adjacent TFs represent a distinct class, but also include several TFs with SLE genetic associations, including NFκB⁴⁵, the Ikaros (IKZF) family⁴⁶, and HMGA family members⁴⁷ (Fig. 5d and Supplementary Data 15). Collectively, these analyses reveal two distinct classes of TFs at a given SLE-associated locus that both likely play key roles in SLE mechanisms, along with particular TFs that tend to participate in one class or the other.Fig. 5

a Model of variant overlapping and variant adjacent transcription factors (TFs). Variant overlapping TFs (blue) allelically bind on top of variants, while variant adjacent TFs (orange) allelically bind near variants. b TF binding site location distribution for variant overlapping (blue) and variant adjacent (orange) TFs, relative to allelic enVars. c TF motif families enriched for participating as variant overlapping TFs at allelic enVars. Motif disruption p-values were estimated by a two-sided proportions test by comparing the fraction of motif disruption events at allelic enVars to the fraction observed at non-allelic enVars (see “Methods”). d TF motif families enriched for participating as variant adjacent TFs at allelic enVars. Motif enrichment p-values were estimated by a two-sided proportions test by comparing the fraction of predicted TF binding sites in allelic enVars to random expectation (see “Methods”). For both the variant overlapping and variant adjacent analyses, motif families are shown with p_adj < 0.0001 and three or more allelic events at allelic enVar loci, or five or more predicted binding sites at allelic enVar loci, respectively.

To explore the cell-type specificity of allelic enhancer activity, we transfected our SLE MPRA library into Jurkat cells, a T cell line, as T cells are another important cell type in SLE². We identified 92 SLE-risk variants as allelic enVars in Jurkat cells, 25% of which were also found in GM12878 (Supplementary Fig. 5a, b and Supplementary Data 11). We then repeated the experiment in Jurkat cells stimulated with the inflammatory cytokine TNFα, a key cytokine in SLE development⁴⁸, to identify stimulation-dependent allelic enVars. This resulted in the identification of 102 allelic enVars, 28 of which were specific to the stimulated Jurkat cells (Supplementary Fig. 5c, d and Supplementary Data 11). Altogether, our study identified a total of 145 allelic enVars across 50 independent SLE-risk loci (Supplementary Fig. 5e). These results highlight allelic transcriptional regulatory mechanisms that are both cell-type and inflammatory signaling-dependent.

In summary, through the application of an allelic MPRA library to the EBV-transformed B cell line GM12878, we identified global transcriptional enhancer activity at 16% of SLE-associated genetic variants (enVars), with particular transcriptional regulatory proteins concentrated at these genomic locations. We further identified 51 SLE-risk variants with allelic enhancer activity (allelic enVars) that we now nominate as plausibly causal by acting through genotype-dependent changes to enhancer activity in GM12878. Upon comparison to allelic enhancer activity in the Jurkat T cell line, we identified shared and unique allelic transcriptional regulatory mechanisms at SLE-risk loci. Using experimental TF ChIP-seq data and TF binding site motif scanning, we propose a model where the collective action of the genotype-dependent binding of particular variant overlapping and variant adjacent TFs leads to genotype-dependent transcriptional activity at SLE-risk loci.

All SLE-associated genetic risk loci reaching genome-wide significance published through March 2018 were included in this study^{31,38,44,58,69–78}. A total of 91 genetic risk loci were used for linkage disequilibrium (LD) expansion (r² > 0.8) based on 1000 Genomes Data⁷⁹ in the ancestry(ies) of the initial genetic association using PLINK(v1.90b)⁸⁰ (Supplementary Data 1). All expanded variants were updated to the dbSNP 151 table⁸¹ from the UCSC table browser⁸² based on either variant name or genomic location. Unmappable variants were discarded. We also included 20 genetic variants from the Tewhey et al.¹⁹ study as positive and negative controls.

For single nucleotide polymorphisms, we pulled 170 base pairs (bps) of hg19-flanking DNA sequences for every allele, with the variant located in the center (84 bps upstream and 85 bps downstream of the variant). For the other types of variants (indels), we designed the flanking sequences to ensure that the longest allele has 170 bps. Adapters (15 bps) were added to each sequence at either end (5′-ACTGGCCGCTTGACG - [170 bp oligo] - CACTGCGGCTCCTGC-3′) to make a 200 bp DNA sequence (Supplementary Data 4). For all resulting sequences, we created a forward and reverse complement sequence to compensate for possible DNA synthesis errors. A total of 12,478 oligos (3093 variants, 6239 alleles) were obtained from Twist Bioscience.

For assembly of the MPRA library, we followed the procedure described by Tewhey et al.¹⁹ with minor modifications. In brief, we first created the empty vector pGL4.23∆xba∆luc from pGL4.23[luc2/minP] using primer Q5_deletion_rev and Q5_deletion_fwd following the manufacturer’s instruction of the Q5 Site-Directed Mutagenesis Kit. Then, 20 bps barcodes were added to the synthesized oligos through 24X PCR with 50 μL system, each containing 1.86 ng oligo, 25 μL NEBNext® Ultra™ II Q5® Master Mix, 1 μM MPRA_v3_F, and MPRA_v3_201_R. PCR was performed under the following conditions: 98 °C for 2 min, 12 cycles of (98 °C for 10 s, 60 °C for 15 s, 72 °C for 45 s), 72 °C for 5 min. Amplified product was purified and cloned into SfiI digested pGL4.23∆xba∆luc by Gibson assembly at 50 °C for 1 h. The assembled backbone library was purified and then transformed into Escherichia coli (E. coli) through electroporation (2 kV, 200 ohm, 25 μF). Electroporated E. coli was expanded in 200 mL of LB Broth buffer supplemented with 100 μg/mL of carbenicillin at 37 °C for 12 to 16 h. Plasmid was then extracted using the QIAGEN Plasmid Maxi Kit.

We next created the pGL4.23[eGFP/miniP] plasmid. An eGFP fragment was amplified from MS2-P65-HSF1_GFP (Addgene #61423) through PCR with a 50 μl system containing 1 ng plasmid, 25 μL NEBNext® Ultra™ II Q5® Master Mix, 0.5 μM GFP_seq_MS2-P65-HSF1_GFP_FWD, and GFP_seq_MS2-P65-HSF1_GFP_REV. PCR was performed under the following conditions: 98 °C for 2 min, 20 cycles of (98 °C for 10 s, 60 °C for 15 s, 72 °C for 30 s), 72 °C for 5 min. The amplified fragment was purified and then inserted into XbaI and NcoI digested pGL4.23[luc2/minP] through Gibson assembly at 50 °C for 1 h. The assembled plasmid was purified and then transformed into E. coli through chemical transformation. Transformed E. coli was expanded in 100 mL of LB Broth buffer supplemented with 100 μg/mL of carbenicillin at 37 °C for 12–16 h. Plasmid was then extracted using the QIAGEN Plasmid Maxi Kit.

A miniP + eGFP fragment was amplified from pGL4.23[eGFP/miniP] through 8X PCR with 50 μL system, each containing 1 ng plasmid, 25 μL NEBNext® Ultra™ II Q5® Master Mix, 0.5 μM 200-MPRA_v3_GFP_Fusion_v2_F, and 201-MPRA_v3_GFP_Fusion_v2_R. PCR was performed under the following conditions: 98 °C for 2 min, 20 cycles of (98 °C for 10 s, 60 °C for 15 s, 72 °C for 45 s), 72 °C for 5 min. The amplified product was purified and then inserted into AsiSI digested backbone library through Gibson assembly at 50 °C for 1.5 h to create the transfection library. The resulting library was re-digested by RecBCD and AsiSI, purified, and then transformed into E. coli through electroporation (2 kV, 200 ohm, 25 μF). Transformed E. coli was cultured in 5 L of LB Broth buffer supplemented with 100 μg/mL of carbenicillin at 37 °C for 12–16 h. The plasmid was then extracted using the QIAGEN Endo-free Plasmid Giga Kit.

The oligo and barcode regions were amplified from the backbone library through 4X PCR with a 100 μL system containing 200 ng plasmid, 50 μL NEBNext® Ultra™ II Q5® Master Mix, 0.5 μM TruSeq_Universal_Adapter_P5, and MPRA_v3_TruSeq_Amp2Sa_F_P7. PCR was performed under the following conditions: 95 °C for 20 s, 6 cycles of (95 °C for 20 s, 62 °C for 15 s, 72 °C for 30 s), 72 °C for 2 min. The product was then purified, and indices were added through a 100 μl system containing all purified product, 50 μl NEBNext® Ultra™ II Q5® Master Mix, 0.5 μM TruSeq_Universal_Adapter_P5, and index primer. PCR was performed as above, except for only five cycles. Samples were purified, molar pooled, and sequenced using 2 × 125 bp on Illumina NextSeq 500.

The GM12878 cell line was grown in RPMI medium supplemented with 10% FBS, 100 units/mL of penicillin, and 100 µg/mL of streptomycin. Cells were seeded at a density of 5 × 10⁵ cells/mL the day before transfection. For triplicate transfections, we collected a total of 5 × 10⁷ cells per replicate. Cells were then suspended with 50 μg transfection library plasmid in 400 μL Buffer R. Electroporation was performed with the Neon transfection system in 100 μl tips with 3 pulses of 1200 V, 20 ms each. After transfection, cells were recovered in 50 mL pre-warmed RPMI medium supplemented only with 10% FBS for 24 h. Cells were then collected for preparation of the sequencing library for barcode counting.

The Jurkat cell line was grown in RPMI medium supplemented with 10% FBS, 100 units/mL of penicillin, and 100 µg/mL of streptomycin. Cells were seeded at a density of 5 × 10⁵ cells/mL the day before transfection. For each experimental group, we collected a total of 5 × 10⁷ cells per replicate for 5 replicates. Cells were then resuspended with 50 μg transfection library plasmid in 400 μL Buffer R. Electroporation was performed with the Neon transfection system in 100 μl tips with 3 pulses of 1350 V, 10 ms each. After transfection, cells were recovered in 50 mL pre-warmed RPMI medium supplemented only with 10% FBS for 24 h. After recovery, cells were supplemented with or without 100 ng/ml TNFα for 24 h. Cells were then collected for preparation of the sequencing library for barcode counting.

For samples from GM12878 cells, total RNA of transfected cells was extracted by the RNeasy Midi Kit following the manufacturer’s instruction. Extracted RNA was subjected to DNase treatment in a 375 μL system with 2.5 μL Turbo DNase and 37.5 μL Turbo DNase Buffer at 37 °C for 1 h. 3.75 μL 10% SDS and 37.5 μL 0.5 M EDTA were added to stop DNase with 5 min of incubation at 75 °C. The whole volume was used for eGFP probe hybridization in an 1800 μL system, with 450 μl 20X SSC buffer, 900 μL Formamide and 1 μL of each 100 μM Biotin-labeled GFP probe One to Three. The probe hybridization was performed through incubation at 65 °C for 2.5 h. 200 μL Dynabeads™ MyOne™ Streptavidin C1 was prepared according to the manufacturer’s instruction. The beads were suspended in 250 μL 20X SSC Buffer and incubated with the above probe hybridization reaction at room temperature for 15 min. Beads were then collected on a magnet and washed with 1X SSC Buffer once, and 0.1X SSC Buffer twice. eGFP mRNA was eluted first through adding 12.5 μL ddH₂O, heating at 70 °C for 2 min and collecting on a magnet, then adding another 12.5 μL ddH₂O, heating at 80 °C for 2 min and collecting on a magnet. All collected elution was performed with another DNase treatment in a 30 μL system containing 0.5 μL Turbo DNase and 3 μL Turbo DNase Buffer at 37 °C for 1 h. 0.5 μL 10% SDS was added to halt DNase reaction. Eluted mRNA was purified through RNA Clean SPRI Beads. mRNA was reverse transcribed to cDNA using SuperScript™ IV First-Strand Synthesis System with gene specific primer MPRA_v3_Amp2Sc_R, following the manufacturer’s instruction. cDNA and plasmid control were then used for building sequencing libraries following the Tag-seq Library Construction section in the paper of Tewhey et al.¹⁹. In brief, 1 μL of cDNA and plasmid control samples were used to estimate the relative concentration of eGFP in the 10 μL system containing 5 μL NEBNext® Ultra™ II Q5® Master Mix, 0.6 μL SYBR green I diluted 1:1000 (Life Technologies, S7563), and 0.5 μM TruSeq_Universal_Adapter_P5 and MPRA_V3_Illumina_GFP_F. PCR was performed under the following conditions: 95 °C for 20 s, 40 cycles of (95 °C for 20 s, 65 °C for 20 s, 72 °C for 30 s), 72 °C for 2 min. According to the cycle threshold, all cDNA and plasmid control samples were diluted to match the sample with the lowest concentration. A total of two PCRs were needed for building the sequencing library. The first PCR was performed with 10 μL of normalized samples in the 50 μL system containing 25 μL NEBNext® Ultra™ II Q5® Master Mix, 0.5 μM TruSeq_Universal_Adapter_P5, and MPRA_V3_Illumina_GFP_F. PCR was performed under the following conditions: 95 °C for 20 s, corresponding cycles of (95 °C for 20 s, 65 °C for 20 s, 72 °C for 30 s), 72 °C for 2 min. The product was then purified, and indices were added through a 100 μl system containing all purified product, 50 μl NEBNext® Ultra™ II Q5® Master Mix, 0.5 μM TruSeq_Universal_Adapter_P5, and index primer. PCR was performed as above, except for only 6 cycles. Samples were purified, molar pooled, and sequenced using 1 × 75 bp on Illumina NextSeq 500.

For samples from Jurkat cells, total DNA and RNA of transfected cells were extracted by the Qiagen ALLPrep DNA/RNA Mini Kit following the manufacturer’s instruction²³. Extracted RNA was processed the same as above to obtain cDNA. cDNA, extracted DNA, and plasmid control were then used for building sequencing libraries with the same protocol described above. Samples were purified, molar pooled, and sequenced using 1 × 100 bp on Illumina NovaSeq 6000.

All primers used in this study are provided in Supplementary Table 1.

Paired-end, 125 bp reads were first quality filtered using Trimmomatic (v0.38)⁸³ (flags: PE -phred33, LEADING:25, TRAILING:25, MINLEN:80). Read 1 was then separated into the 20 bp barcode region and the oligo-matching region. The trimmed oligo-matching regions of Read 1 and Read 2 were mapped back to the synthesized oligo sequences using Bowtie2 (v2.3.4.1)⁸⁴ (flags: -X 250, -very-sensitive, -p 16). Barcodes were then associated with the oligo sequences using the read ID. Only uniquely mapped barcodes were used for downstream analysis.

Single-end 75/100 bp reads were first quality filtered using Trimmomatic (v0.38)⁸³ (flags: PE -phred33, LEADING:3, TRAILING:3, MINLEN:70). Each read was then separated into the 20 bp barcode region and the constant region. The trimmed constant regions of the reads were mapped back to the constant region within the eGFP 3′ UTR using Bowtie2 (v2.3.4.1)⁸⁴ (flags: -very-sensitive, -p 16). Only reads with Levenshtein distance of 4 or less within the constant region and perfect matches to the two bases directly adjacent to the barcode were kept. Barcodes were then associated with the retained reads using the read ID. Only barcodes that met our quality threshold requirements described above in the Methods section “Oligo and barcode association” were used for downstream analysis.

We followed the procedures described in the “Identification of Regulatory Oligos” section of Tewhey et al.¹⁹ with minor modifications. In brief, oligos (alleles) with 30 or more unique barcodes from the plasmid control were included for analysis. All barcodes were summarized at the oligo level. Barcode count totals for each oligo, including all SLE variants and the 20 control variants, were passed into DESeq2 (v1.28.1)⁸⁵ in R (v3.5.3) to estimate the fold change and significance between plasmid controls (Supplementary Note 1 and Supplementary Fig. 5f, g) and the experimental replicates. A Benjamini–Hochberg FDR adjusted p-value of <0.05 was required for significance. Only significant alleles with greater than or equal to a 1.5x fold change were identified as enhancer alleles (enAlleles). A variant was identified as an enhancer variant (enVar) if any allele of this variant was an enAllele. Results for the 20 control variants were compared to data from Tewhey et al.¹⁹ to estimate accuracy, sensitivity, and specificity.

Only enVars were considered for allelic analysis. The barcode counts from every allele of each enVar were used for calculating p-values by comparing the log2 ratios of the non-reference allele vs the reference allele, normalized by plasmid controls, using Student’s t-test^19,35. p-values were adjusted with the Benjamini–Hochberg FDR-based procedure. A corrected p-value of <0.05 was required for significance. Only significant alleles with 25%-fold changes or greater were identified as allelic enVars (Supplementary Note 2 and Supplementary Fig. 6). We have created an R package (mpraprofiler) for performing this analysis, which is available on the Weirauch lab GitHub page (https://github.com/WeirauchLab/mpraprofiler).

We annotated each SLE genetic variant with its nearest gene using the NCBI RefSeq table⁸⁶ downloaded from the UCSC table browser⁸². enVars were annotated using a combination of DNA looping interactions (GM12878 Capture Hi-C data^87,88) and eQTL data obtained from the eQTL Catalog, a resource that contains quality-controlled, uniformly re-computed eQTLs from 19 eQTL publications^61,89–108, EBV-transformed B cell lines (GTEx Analysis V7 (dbGaP Accession phs000424.v7.p2))⁴⁰ and other individual studies^109–112. For all variants, the target genes were annotated (Supplementary Data 2) using the union of promoter interacting genes and eQTL genes from B cells with and without EBV transformation, when available. Otherwise, target genes were annotated as the nearest gene. Allelic enVar gene targets were classified into four tiers: a Tier (1) variant is both an eQTL and also loops to the promoter of the same gene; a Tier (2) variant has an eQTL for at least one gene; a Tier (3) variant only loops to the promoter of at least one gene; a Tier (4) variant is neither an eQTL nor loops to the promoter of any gene (Supplementary Data 12).

To identify specific TFs whose binding might contribute to the enhancer activity observed in our MPRA experiments, we performed HOMER (v4.9)³³ TF binding site motif enrichment analysis. Specifically, we used HOMER to calculate the enrichment of each motif in the sequence of enAlleles compared to the sequences of non-enAlleles. HOMER was modified to use the large library of human position weight matrix (PWM) binding site models contained in build 2.0 of the Cis-BP database³⁴ and a log base 2 likelihood scoring system.

Enrichr^113,114 was used for GO enrichment analysis. In short, the target genes of enVars were passed to Enrichr for analysis. Results from the GO biological process (2018) category were used (Supplementary Data 8, Supplementary Fig. 3).

1058 ChIP-seq datasets were obtained from the Gene Expression Omnibus (GEO)¹¹⁵ using custom scripts that searched for ChIP-seq experiments performed in EBV-transformed lymphoblastoid cell lines (LCLs). The annotations for every dataset (assay type, cell line, assayed molecule) were manually checked by two authors (MTW and LCK) to ensure accuracy. The Sequence Read Archive (SRA) files obtained from GEO were analyzed using an automated pipeline. Briefly, the pipeline first runs QC on the FastQ files containing the sequencing reads using FastQC (v0.11.2). If FastQC detects adapter sequences, the pipeline runs the FastQ files through Trim Galore (v0.4.2)¹¹⁶, a wrapper script that runs cutadapt (v1.9.1)¹¹⁷ to remove the detected adapter sequence from the reads. The quality-controlled reads are then aligned to the reference human genome (hg19/GRCh37) using bowtie2 (v2.3.4.1)⁸⁴. The aligned reads (in .BAM format) are then sorted using samtools (v1.8.0)¹¹⁸ and duplicate reads are removed using picard (v1.89)¹¹⁹. Finally, peaks are called using MACS2 (v2.1.2) (flags: callpeak -g hs -q 0.01 -f BAM)⁴⁴. ENCODE blacklist regions¹²⁰ were removed from the peak sets using the hg19-blacklist.bed.gz file available at https://github.com/Boyle-Lab/Blacklist/tree/master/lists/Blacklist_v1. ChIP-seq datasets GSM1666207, GSM2748907, and GSM1599157 were removed due to the low number of cells used in the experiments.

We used the RELI(v0.9)¹⁰ algorithm to identify genomic features (TF binding events, histone marks, etc.) that coincide with enVars. As input, RELI takes the genomic coordinates of enVars. RELI then systematically intersects these coordinates with one of the GM12878 ChIP-seq datasets, and the number of input regions overlapping the peaks of this dataset (by at least one base) is counted. Next, a p-value describing the significance of this overlap is estimated using a simulation-based procedure. To this end, a ‘negative set’ is created for comparison to the input set, which in this study contains the set of non-enVars (i.e., variants with no allele having an adjusted p-value of <0.05 and more than 10%-fold change in the DESeq2 result). A distribution of expected overlap values is then created from 2000 iterations of randomly sampling from the negative set, each time choosing a set of negative examples that match the input set in terms of the total number of genomic loci. The distribution of the expected overlap values from the randomized data resembles a normal distribution and can thus be used to generate a Z-score and corresponding p-value estimating the significance of the observed number of input regions that overlap each ChIP-seq dataset.

We performed similar RELI analysis for allelic enVars. As input, we used the allelic enVar sites. For the ‘negative set’, we used the set of common SNPs taken from the dbSNP142 database downloaded from the UCSC table browser⁸².

To identify possible mechanisms underlying our allelic enVars, we applied our MARIO(v3.93)¹⁰ method to the LCL ChIP-seq dataset collection described above. In brief, MARIO identifies common genetic variants that are (1) heterozygous in the assayed cell line and (2) located within a peak in a given ChIP-seq dataset. It then examines the sequencing reads that map to each heterozygote in each peak for imbalance between the two alleles. Results are combined across experimental replicates to produce a robust Allelic Reproducibility Score (ARS). Results with MARIO ARS value >0.4 that also pass the following three post-processing filters were considered allelic. (1) The variant must be significantly allelic for a given protein/histone mark (ARS > 0.4) in at least 50% of the datasets in which that variant was heterozygous; (2) The same allele must be significantly preferred (ARS > 0.4) in at least 75% of the datasets where that variant shows significant allelic behavior; and (3) The replicates of a given experiment must all prefer the same strong allele. These post-processing filters were applied to remove results with inconsistent allelic imbalance, extending the procedures of our previous study¹⁰.

Variant overlapping TFs were identified using an algorithm that compares predicted TF binding motif scores between the different alleles of each allelic enVar. First, we padded each allele of a given allelic enVar with 25 bps of upstream and downstream DNA sequence (a sufficient length to account for any known human TF binding sites¹²¹). The algorithm consists of two major components: (1) individually scoring the two alleles of a given variant with a given TF model; and (2) quantifying the difference in the binding intensity between these two alleles. DNA sequences are scored using the large collection of human TF position  weight matrix (PWM) models contained in the Cis-BP database³⁴ and the log-likelihood PWM scoring system¹²². Since log-likelihood score distributions vary substantially (depending on the information content of a given motif), we employ a simple scaled scoring system that maps a given log-likelihood score to the percentage of the maximum achievable log-likelihood score of the given motif—we refer to this value as the “relative PWM score”. We identify binding site altering events (i.e., “creating” or “breaking” a predicted binding site for a given TF motif) as cases where one allele has a relative PWM score of 70% or higher, and the other allele has a score of <40%. For a given variant, any TF with allelic ChIP-seq sequencing reads (see above) and a binding site altering event for any of its motifs was deemed a variant overlapping TF. Any TF with allelic ChIP-seq sequencing reads and a lack of a binding site altering event for any of its motifs was deemed a variant adjacent TF.

We next sought to identify particular TFs that tend to be variant overlapping TFs at SLE allelic enVars. To this end, we calculated the fraction of times each TF motif has a binding site altering event (as defined above) at SLE allelic enVars. As background, we calculated the fraction of times each TF motif has a binding site altering event at non-allelic enVars. The significance of the difference between these two fractions was then calculated using a proportions test. Results are provided in Supplementary Data 15.

We used a similar procedure to identify particular TFs that tend to be variant adjacent TFs at SLE allelic enVars. We performed HOMER(v4.9) motif enrichment analysis using the full 170 bp allelic enVar DNA sequences as input. The dinucleotide scrambled version of input sequences were used as background. The fractions of motif “hits” obtained in the foreground vs. background set were then compared, and significance was again calculated using a proportions test. Results are provided in Supplementary Data 15.

Further information on research design is available in the Nature Research Reporting Summary linked to this article.