Abstract
Background Systemic lupus erythematosus (SLE) is a chronic autoimmune disease with a strong genetic component (Block, Winfield et al. 1975, Deapen, Escalante et al. 1992, Alarcon- Segovia, Alarcon-Riquelme et al. 2005) and most of the genes involved have modest impact, complicates genetic dissection of disease etiology (Eberle, Ng et al. 2007, Stahl, Wegmann et al. 2012). To overcome this, a sensitive approach to identify risk alleles with small gene effects called genome-wide association studies (GWAS) has been applied for SLE studies and more than one hundred loci have been linked with disease risk (Gateva, Sandling et al. 2009, Bentham, Morris et al. 2015, Morris, Sheng et al. 2016, Langefeld, Ainsworth et al. 2017, Zhang, Zhang et al. 2018). Therefore, decoding GWAS is a promising strategy to identify novel drug targets in SLE. However, translating GWAS hits into disease pathogenic mechanisms has proven difficult as most of the identified disease-associated hits are noncoding single-nucleotide polymorphisms (SNPs), suggesting a regulatory role, and cannot be distinguished from others that reside incidentally within risk loci due to Linkage Disequilibrium (LD), which significantly increase the number of SNPs need to be investigated. A high-throughput screen method is therefore ideal for this task. Massively parallel reporter assays (MPRAs) that test the translational impact of candidate SNPs have been applied for variants library associated with different traits, including one lupus study (Tewhey, Kotliar et al. 2016, Ulirsch, Nandakumar et al. 2016, Lu, Chen et al. 2021). Here we show another lupus MPRA using a different B cell line.
Methods A library of 151 bp oligo pool including each alleles of 2342 SNPs (2176 candidate lupus associated SNPs covering 87 linked loci, 80 positive control SNPs and 86 negative control SNPs) was cloned into a luciferase report gene expressing plasmid, upstream of the minimal promoter to test the enhancer capability of each allele. Each allele was also barcoded with 20 different short oligo sequences downstream of luciferase gene to avoid sequencing preference and have higher statistical power. In the end, a library of 93,840 plasmid constructs was nucleofected into Daudi B cells (three biological repeats). Total RNAs were extracted 24 hours after nucleofection and sent for next generation sequencing (NGS) along with the plasmid library. Alleles with better enhancer role will have higher total barcode counts. Oligos (alleles) with all 20 unique barcodes from the sequencing data were included for analysis. Barcode count totals for each oligo, including all SLE variants and the control variants, were passed into DESeq2 in R to estimate the fold change and difference between each allele of the same SNP. A Benjamini–Hochberg FDR adjusted p-value smaller than 0.05 were required to show significant allelic difference. Allelic enhancer variants were then further analyzed using Regulatory Element Locus Intersection (RELI) and Measurement of Allelic Ratios Informatics Operator (MARIO) bioinformatic analysis to identify transcriptional factor binding to the SNP.
Results Principle component analysis and correlation analysis showed good quality of NGS sequencing data (figure 1A and 1B), from which 34 SNPs are found to be allelic enhancer variants from a total of 2342, including 5 from the positive control, 29 from the SLE library and none from the negative control (figure 1C). Furthermore, with RELI, we observed significant enrichment for overlap between the 29 identified lupus allelic enhancer variants and histone modification marks, including H3K27ac and H3K4me1 from 576 GM12878 ChIP-seq datasets available in the NCBI Gene Expression Omnibus (GEO) database. As expected, we did not identify enrichment for repressive marks such as H3K9me3 or H3K27me3 (figure 2A). Altogether, the genomic features present confirm these SNPs likely alter transcriptional regulation. As to which transcriptional factor binds to those SNPs and regulate gene expression, MARIO pipeline searches for allelic binding events at these SLE variants within 576 GM12878 ChIP-seq datasets, rs205286 has incredibly strong and reproducible allele-dependent binding of CCCTC- binding factor (CTCF) and Cohesin subunit SA-1 (STAG1) to the C allele (figure 2B).
Therefore, STAG1 and CTCF might work together to bind to rs205286 to change the formation of DNA loop while transcription, hence regulate associated gene expression.
Conclusion Using an unbiased experimental high-throughput screen method massively parallel reporter assay, we successfully identified 29 SNPs with allelic regulatory activity, represent 1.3% from initial screening library. The integrative bioinformatic analyses further reveal one possible molecular mechanism underlying genotype-dependent transcriptional regulation. Therefore, we conclude that the MPRA is a robust tool for dissecting the genetic etiology of SLE.
