We analyzed publicly available RNA-seq data for B cell differentiation (GSE148924) and human tissues (PRJEB4337). Reads were aligned to the human genome (GRCh38/hg38) using Salmon (v1.8.0) on default settings. Gene counting was then performed with featureCounts (version 1.6.4) and the resulting count matrix was used as input for DESeq2 (version 1.28.1) (24). Differential gene expression analysis was performed by running DESeq2 with fitType = “local”. Log2 fold change (FC) shrinkage was performed by the adaptive shrinkage algorithm to reduce the effect of low expression levels on fold changes. Significantly differential genes were identified by filtering on Benjamini-Hochberg adjusted P-value (Padj) < 0.001 and absolute log2FC > 0.5 (or < -0.5). For heatmaps, mean normalized expression counts were calculated for each condition. Pairwise correlation of columns was performed using Pearson. Data analysis was performed using Python 3.7.11 with scientific libraries pandas (25), NumPy (26), matplotlib (27), and seaborn. Volcano plots and dot plots were generated using R v. 4.2.1. with libraries tidyverse, ggrepel, dplyr, and ggplot2.

For alternative splicing (AS) analysis reads were aligned to the human genome (GRCh38/hg38) using STAR (version 2.7.5b) on default settings. Differential AS events were determined using rMATS 4.1.2 (28). Sashimi blots were created using the IGV genome browser.

Gene ontology (GO) term enrichment was performed on gene symbols using the enrichGO function in R v. 4.2.1. with libraries clusterProfiler (29, 30), org.Hs.eg.db (31), and AnnotationDbi (32). Only GO enrichments with a false discovery rate (FDR) < 0.05 are presented.

The ShinyGO 0.76 tool (http://bioinformatics.sdstate.edu/go/) was used to identify transcription factor (TF) binding motifs enriched in the promoters of target genes.

The protein 3D structure of the different isoforms was predicted using Robetta (https://robetta.bakerlab.org/) and visualized using PyMol v.2.3.4.

To evaluate differences in gene expression (GE) and alternative splicing (AS) during B cell activation and differentiation, we analyzed publicly available RNA-seq data of preplasmablast (prePB), plasmablast (PB), plasma cells (PC), and memory B cells (MBC) (GSE148924, Figure 1A ). Our principal component analysis (PCA) of the 15.754 gene expression levels showed that each cell subpopulation was segregated according to its developmental stage and that the triplicate samples of the individual cell types clustered well together (as previously demonstrated in 33) ( Figure 1B ). This analysis rules out high variability or heterogeneity of gene expression within one sample group and we therefore proceed using the triplicates in our analysis. We used DeSeq2 to compare changes in gene expression using a log2FC > 0.5 (< -0.5) and Padj <0.001 as cut-offs. Overall, we observed a similar gene expression profile between PB and PC ( Figure 1C ), where only 125 genes were differentially expressed between them. We observed a stronger difference when comparing prePB and PB, and prePB and PC (2383 and 3079 differentially expressed genes, respectively). An even stronger difference in gene expression was observed when comparing prePB, PB, and PC to MBC. We identified 6362 differentially expressed genes between prePB and MBC, 6144 differentially expressed genes between PB and MBC, and 4920 differentially expressed genes between PC and MBC ( Figure 1C ).

Furthermore, we performed a GO term enrichment of differentially expressed genes amongst the different stages of B cell development. Interestingly, the GO terms mRNA processing (GO:0006397) and RNA splicing (GO:0008380) were highly enriched in MBC when compared to PB (FDR < 0.05) ( Figure 1D ). Although the same GO terms were not enriched when comparing MBC to PC, we observed that for “mRNA processing” 88 out of 410 genes were upregulated in MBC, while only 36 were upregulated in PCs. For “RNA splicing”, 75 genes were upregulated in MBC, and 25 were upregulated in PCs (342 total) ( Figures 2A, B ), suggesting that there is also a substantial difference in mRNA processing and splicing in MBC when compared to PC.

It is well known that PCs are responsible for secreting antibodies, and therefore, have a higher rate of protein biosynthesis than MBCs (4). To analyze whether this is accompanied by an adaptation of the secretory pathway at the level of gene expression, we compared the expression of intracellular trafficking related-genes between MBC and PC by analyzing genes included in the following GO terms: intracellular protein transport (GO:0032779), COPII-coated vesicle budding (GO:0090114), retrograde vesicle-mediated transport Golgi to ER (GO:0048220) and protein exit from ER (GO:0032527) ( Figures 2C–F ). We observed that 431 out of 1257 genes related to intracellular protein transport are differentially regulated between MBC and PC without a significant bias towards up- or down-regulation in either of the two cell types. That is about the expected number, given that the expression of about 1/3 of all expressed genes changes between these B cell subtypes (4920 out of 15.754 genes, see above, Figure 2C ). However, when we focus on the early secretory pathway by looking at GO terms COPII-coated vesicle budding, Retrograde vesicle-mediated transport Golgi to ER, Protein exit from ER, and ER to Golgi vesicle-mediated transport, we find strong upregulation of many genes almost exclusively in PC ( Figures 2D–F ; Supplementary Figure 1 ). In contrast, the GO-terms trans-Golgi network membrane and Clathrin-related genes are not specifically regulated between the two cell types ( Supplementary Figure 1 ). These data indicate that especially the expression of components of the early secretory pathway is upregulated in antibody-secreting B cell subtypes to accommodate the increased secretory cargo flux.

To better understand the mechanism that allows PCs to regulate their secretory machinery, we compared the expression level changes of COPII components in PCs and during the activation of CD4+ T cells ( Figure 2G ) (34). Interestingly, the expression of COPII components during T cell activation stays rather constant, whereas PCs show a substantial upregulation of almost all COPII components when compared to non-antibody-secreting prePB. These data suggest that B cells use changes in gene expression to adjust their secretory machinery during activation and differentiation, whereas T cells rely mostly on alternative splicing of key regulators (34, 35). This may also indicate that B cells require stronger adjustments to their secretory system than T cells, which would make B cell differentiation a very interesting model system to study the adaptation of the secretory pathway in the future.

We then performed a global alternative splicing analysis and detected a total of 44.879 AS events, with 3179 events showing significant differences between the differentiation stages we analyzed. All significantly changed alternative splice events are shown in the heatmap in Figure 3A (P-value < 0.001, DPSI > |0.4|). The significantly changed AS events fall into the following four categories: 2683 (80.07%) exon skipping (SE), 53 (2.16%) intron retention (RI), 209 (8.18%) alternative 5 splice site (A5SS), and 234 (9.58%) alternative 3 splice site (A3SS). As observed in our GE analysis, the AS pattern was rather similar amongst prePB, PB, and PC, while MBC showed substantial differences. The splicing PCA looks similar to the GE PCA, with close relation of PB and PC, which are separated from prePB and MBC ( Figure 3B ). The PCA for splicing events shows a high degree of homogeneity in PB and PC and slightly lower in preBP. Only MBCs show a larger extent of heterogeneity, which is still low considering the distance of the MBC cluster to the other cell types. The vastly different global splicing pattern of MBS is consistent with the gene expression analysis showing enrichment of GO terms related to splicing and RNA processing when comparing genes differentially expressed between MBC and PC ( Figures 2A, B ). Notably, ¾ of all significantly changed splicing events are less included in MBCs, which may point to reduced overall splicing efficiency.

Further, we perform a GO term enrichment analysis of differentially spliced genes in MBC when compared to PCs. In addition to the biological processes involved in post−transcriptional regulation of gene expression, we found an enrichment of the GO terms Golgi vesicle transport (GO:0048193), nuclear transport (GO:0051169), and nucleocytoplasmic transport (GO:0006913), suggesting that genes involved in intracellular transport are heavily alternatively spliced in MBC ( Figure 3C ) (full GO term enrichment analysis in Supplementary Figure 2 ). The global changes in alternative splicing in MBC suggests a substantial contribution in controlling the functional changes during the formation of memory B cells and the specific enrichment of GO terms suggests that alternative splicing plays a fundamental role in controlling intracellular trafficking.

To get a more detailed breakdown of intracellular trafficking, we concentrated on the COPII machinery, which controls the first step in protein secretion. In the GO term COPI-coated vesicle budding, we found 21 alternatively spliced events, which mostly show a MBC unique splicing pattern ( Figure 3D ). Important COPII components like SEC13, SEC31A, and SEC24C have multiple splicing events that are differentially spliced in MBCs compared to the other cell types. For example, SEC24C exon 11 is 100% included in prePB, PB, and PC, and only ~50% included in MBC ( Figure 3E ). The length of exon 11 is 116 nucleotides, therefore exclusion of this exon causes a frameshift and a premature stop codon. It is very probable that through nonsense-mediated decay (NMD) this splicing event leads to a lower abundance of SEC24C in MBCs (36), which we detected in our gene expression data in MBCs ( Figure 2D ). SEC24C has a specific IxM cargo binding site and positively regulates the export of the Q-SNARES Syntaxin5, GS27, and Bet1 (11). Therefore, the expression level of SEC24C may contribute to cargo specificity in a cell type-specific manner. SEC24D also contains an alternatively spliced, frameshift-inducing exon (296 nucleotides in length, after exon 9) that is almost exclusively included in MBC ( Figure 3E ). This again correlates with reduced expression of SEC24D in MBC, suggesting that AS-NMD mediated control of gene expression is a more general mechanism to control components of the trafficking machinery in a cell type-specific manner.

While some COPII components are controlled through NMD-inducing isoforms, other splice variants do not induce a frameshift and NMD and thus can control the functionality of the respective protein. For example, the COPII outer coat component SEC31A is heavily alternatively spliced. MBCs include a 50 nucleotides small non-conserved Exon after exon 3. ( Figure 3E ). This alternative exon is likely coupled to an alternative translation start and a shorter N-terminus resulting in the truncation of the first WD40 repeat. The WD40 domain does interact with the WD40 domain of SEC13, the second component of the outer COPII layer (37). Truncation of the SEC31A WD40 domain may thus lead to reduced interaction between SEC31A and SEC13 and may control the efficiency of coat formation or coat structure. Structural analysis of COPII coats formed by splice variants of the individual COPII components will be of high interest to address their cell type-specific functionality in future research.

We have previously shown that alternative splicing of the COPII scaffolding protein SEC16A, responsible for recruiting the COPII machinery, controls the efficiency of ER export during T cell activation (35). This study demonstrates that SEC16A exon 29 is predominantly included in activated T cells, which then increases export efficiency. In B-lymphocytes, we observe a similar splicing pattern. Exon 29 is highly excluded in MBC (PSI = ~15%), whereas in the other differentiation stages it is mostly included (PSI = ~80%) ( Figure 3F ).

KLHL12 is likewise alternatively spliced. Here the MBC exclusive isoform is generated through the exclusion of exons 10 and 11 ( Figure 4A ). While exon 10 exclusion with 99 nucleotides would not shift the reading frame, exon 11 has a length of 187 nucleotides, resulting in a frameshift and an alternative stop codon in the last exon, therefore not inducing NMD. Altogether these splice events, if translated, promote a short C-terminus in the KLHL12 protein ( Figure 4B ). The resulting isoform is missing three out of six kelch repeats, which usually ensure the accurate ubiquitination of SEC31. Without a complete kelch domain, KLHL12 likely is not fully capable of ubiquitinating SEC31, which will decrease the ability to create flexible COPII coats and may result in reduced ability to secrete large cargoes in MBC. Together with the alternative N-terminus of SEC31A discussed above, these splicing-induced changes may lead to substantially altered SEC31A form and function in MBC.

The phosphatidylinositol-binding clathrin assembly protein PICALM regulates clathrin-mediated endocytosis (38) and the size of clathrin-coated vesicles (39). PICALM regulates clathrin assembly through binding clathrin triskelia at its C-terminus, which contains at least two binding sites (40, 41). We find two alternative splice events in PICALM, one of them impacts the clathrin binding site. Exon 13 displays higher inclusion in MBCs, and reduced inclusion in prePB, PB, and PC. On the other hand, exon 17 is basically not included in MBC, has a slight inclusion in prePB, and has a PSI of 0.6 and 0.7 in PB and PC, respectively ( Figure 4C ). The Δ17 isoform codes for the Uniprot isoform 1 (Q13492-1) and Δ13 codes for Uniprot isoform 3 (Q13492-3). Isoform 1 has a shorter C-terminal domain, and therefore, is more efficient in bringing two clathrin triskelia together, thus promoting endocytosis. Through alternative splicing, this effect appears to be enhanced in MBC, weaker in prePB and PB and further reduced in PC ( Figure 4D ).

In our analysis, we found three genes (SEC24B, CUL3, and CSNK1D) in the GO-term COPII-coated vesicle budding, which are upregulated in MBCs and are not regulated through alternative splicing. In order to identify a transcription factor that may be responsible for their MBC-specific upregulation, we performed a Pearson correlation analysis. We used a list of transcription factors (TFs) from the GO terms transcription factor binding (GO:0008134) and DNA-binding transcription factor activity, RNA polymerase II-specific (GO:0000981) and calculated the Pearson correlation index between the three COPII genes and all the TFs from the mentioned GO-terms. Using a threshold of the Pearson correlation coefficient of >0.9, we found 236 TFs that are highly correlated with at least one of the three genes. In the next step we identified 135 TFs that highly correlate with all three, SEC24B, CUL3, and CSNK1D, in MBCs. To validate these findings, we performed the same correlation analysis with the 135 TFs and SEC24B, CUL3, and CSNK1D genes in publicly available RNA-seq data from 17 human tissues (SRA study ERP003613). We then averaged the resulting coefficient indices of SEC24B, CUL3, and CSNK1D and each TF. Additionally, the p-values for the correlation of TFs with each of the three genes in all tissues were calculated and averaged ( Figure 5A ). As a result, we found a small group of TFs, whose expression correlates with the expression of SEC24B, CUL3, and CSNK1D in human tissues. The TF with the highest correlation and the lowest P-value amongst all tissues was NERF/ELF2 ( Figure 5B ). We performed a similar analysis for genes in the GO term COPII-coated vesicle budding, which are upregulated in PC but did not find a strong correlation with a single transcription factor ( Supplementary Table 2 ). Further, we used the ShinyGO 0.76 tool to detect enriched TF binding motifs in the promoters of genes upregulated in PCs. We detected 30 enriched TF binding motifs, but only two of the corresponding TFs (TET1 and E2F6) were also upregulated in PCs when compared to MBCs ( Supplementary Table 3 ). These data suggest that NERF/ELF2 controls the expression of COPII-related genes in MBC, whereas the specific expression of COPII-related genes in PC is not mediated by altered expression of one individual transcription factor and may be dependent on posttranscriptional or posttranslational mechanisms that can control individual targets or the activity of TFs.

Consistent with the hypothesis that NERF/ELF2 controls the expression of trafficking-related genes in MBC, the RNA-Seq data show substantially higher expression of NERF/ELF2 in MBC compared to the other B cell subtypes ( Figure 5C ). In addition to changed gene expression, NERF/ELF2 is also regulated by alternative splicing, as it is expressed in four main protein-coding isoforms. Both NERF-1a and NERF-1b isoforms start with exon 5, resulting in a shorter N-terminus. Additionally, NERF-1a and NERF-1b differ in the use of an alternative 3SS in exon 7 ( Figure 5D ). Similarly, NERF-2a and NERF-2b differ by the use of the exon 7 alternative 3SS but start with exon 1, which adds a 60 amino acid longer N-terminus ( Figure 6A ). The different splice variants are highly regulated in the different B cell subtypes ( Figure 6B ). In MBCs we observe complete skipping of the exon 7 3SS and some usage of the longer N-terminus leading to the expression of NERF-1a and NERF-2a ( Figure 6C ). In contrast, PrePB, PB, and PC have very similar gene expression levels and splice isoform patterns. Around 90% of their isoforms start with exon 5, therefore mostly expressing the short N-terminus isoforms, which are equally distributed between NERF-1a and NERF-1b ( Figure 6C ). This analysis identifies MBC-specific expression of NERF/ELF, regulated at the level of gene expression and alternative splicing, and suggests a role in controlling the MBC-specific expression of intracellular trafficking genes. Future research will address specific functions of different NERF/ELF isoforms and a global role in controlling MBC functionality.