A multi-omics resource of B cell activation reveals genetic mechanisms for immune-mediated diseases

Abstract

Most genetic variants that confer risk of complex immune-mediated diseases (IMDs) affect gene regulation in specific cell types. Their target genes and focus cell types are often unknown, partially because some effects are hidden in untested cell states. B cells play central roles in IMDs, including autoimmune, allergic, infectious, and cancer-related diseases. Despite this established importance, B cell activation states are underrepresented in functional genomics studies. In this study, we obtained B cells from 26 healthy female donors and stimulated them in vitro with six activation conditions targeting key pathways: the B cell receptor (BCR), Toll-like receptor 7 (TLR7), TLR9, CD40, and a cocktail that promotes differentiation into double negative 2 (DN2) IgD^- CD27^- CD11c⁺ CD21^- B cells, a likely pathogenic subset implicated in autoimmunity and infection. We profiled up to 24 B cell activation states and up to 5 control conditions using RNA-seq, single-cell RNA-seq with surface protein markers (CITE-seq), and ATAC-seq. We characterize how IMD-associated genes respond to stimuli and group into distinct functional programs. High-depth RNA-seq data reveals widespread splicing effects during B cell activation. Using single-cell data, we describe stimulus-dependent B cell fates. Chromatin data reveal transcription factors likely involved in B cell activation, and activation-dependent open chromatin regions that are enriched in IMD genetic risk. We experimentally validate a lupus risk variant in a stimulus-specific open chromatin region that regulates TNFSF4 expression, highlighting the relevance of studying B cell activation to elucidate disease association. These data are shared via an interactive browser that can be used to query the dynamics of gene regulation and B cell differentiation during activation by different stimuli, enhancing further investigation of B cells and their role in IMDs: https://mgalab.shinyapps.io/bcellactivation.

Figure 1:. Experimental design and depiction of B cell activation resource.

Illustration of experimental design, activation conditions used, and assays performed. At the bottom panel, indicated are the number of donors (N) used for each assay, and the activation or control conditions used, which are color coded as indicated in the middle panel. The three donors for which there is ATAC-seq data were profiled for standard bulk RNA-seq as well. “0 hours” indicates untouched isolated B cells not put on media. Gray conditions of 4 and 24 hours indicate B cells in media only (“Unstim”). Created with BioRender.com.

Figure 2:. Low-input RNA-seq of 24 activation and 5 control conditions shows pathway-dependent gene expression and transcriptional dynamics in activated B cells over time.

a, Principal component analysis (PCA) with the top 2,000 most variable genes shows separation of conditions and time points. Inset: PCA plot highlighting conditions that include BCR stimulation. b, Differential gene expression across all pairs of conditions and time points. c, Expression levels of selected B cell activation genes. X-axis represents time and colors represent conditions as in (a). d, Modules of gene expression identified by Weighted Gene Correlation Network Analysis (WGCNA). X-axis: time in hours. Y-axis: module eigengene expression values from WGCNA (eigengene is considered a representative of the gene expression profiles in a module). Main hub genes for each module are shown on the right-hand side of each panel, colored by kIM (WGCNA’s intra-module connectivity). e, Enrichment of gene ontology (GO) biological processes for modules in (d). f, Expression levels of selected IMD-associated genes belonging to modules in (d).

Figure 3:. Single-cell RNA expression coupled with quantification of 137 surface protein markers reveals pathway-dependent B cell fates.

a, Uniform Manifold Approximation and Projection (UMAP) based on gene expression after Harmony correction. b, UMAP as in (a) colored by transcriptome-based clusters. c, UMAP as in (a) colored by protein expression levels of selected surface markers. d, RNA expression levels of selected cluster marker genes. e, Proportion of B cell subsets defined by expression of marker protein/genes (genes in italic) within each condition and time point (colors as in (a)). f, UMAP as in (a) colored by RNA expression of selected IMD-associated genes.

Figure 4:. Open chromatin profiling in activated B cells using bulk ATAC-seq uncovers disease-relevant cell states.

a, Principal component analysis on the 5,000 most variable peaks. b, Pairwise number of differentially accessible peaks across all conditions (1% FDR). In orange, peaks that are more accessible in condition in the rows in comparison to condition in the columns. In purple, peaks that are less accessible in condition in the columns in comparison to conditions in the rows. c, Transcription factor motif enrichment in differentially accessible chromatin regions in stimulated conditions with respect to unstimulated cells at 24 hours. Colors indicate log₂ fold change, and size indicates −log₁₀ p-values from HOMER. d, Enrichment of heritability for IMDs (top panel) and control traits (bottom panel) in differentially accessible chromatin regions in stimulated conditions with respect to unstimulated cells at 24 hours, estimated with LDSC-SEG. Tau*: normalized effect size. Asterisks indicate significance at 1% FDR.

Figure 5:. SLE risk variant overlapping DN2c-specific open chromatin region is in an enhancer that regulates TNFSF4.

a, GWAS data from Langefeld et al. rs2205960 is the most likely causal variant. PIP: posterior inclusion probability from SuSiE. b, Chromatin peaks at the genomic region in (a). c, TNFSF4 gene expression in the same conditions shown in (b). TPM: transcripts per million. d, CRISPRi shows that inhibiting the region containing the putative causal variant leads to down-regulation of TNFSF4 RNA expression.

Figure 6:. Pervasive alternative splicing changes during B cell activation.

a, Principal component analysis on gene expression data from 16 donors. b, Number of differentially spliced genes at 5% FDR between all pairs of conditions. c, Splicing effect size in ΔPSI (percent spliced in) between each activation condition and un-stim by −log₁₀ of the p-value. d, Gene ontology biological processes enriched in differentially spliced genes in (c). P-value and enrichment scores from the fgsea R package. e, Differential splicing event in the CD86 gene between DN2c and un-stim. f, Differential splicing event in the ZBTB38 gene between BCRc and DN2c. Numbers indicate average proportion of junction usage, and linewidths are scaled to this usage.