The SLE blood samples were obtained from subjects that met the American College of Rheumatology revised criteria for the classification of SLE and were followed under the Studies of the Pathogenesis and Natural History of Systemic Lupus Erythematosus (SLE) at the National Institute of Arthritis and Musculoskeletal and Skin Diseases of the National Institutes of Health, Bethesda, MD (30). The studies were approved by the Institutional Review Board of the National Institute of Arthritis and Musculoskeletal and Skin Diseases (protocol 94-AR-0066, and 00-AR-0222, respectively). The demographics and clinical characteristics of these SLE blood samples are listed in Supplementary Table 1 (transcriptomic analysis) (1) and Supplementary Table 2 (IgH repertoire analysis). B cell sorting phenotypes are shown in Supplementary Figure 1 have been described previously (1); the purity of the sorted populations was routinely >90%. Peripheral blood B cells from SLE patients (n = 4 independent samples) were isolated using a human B cell enrichment kit (StemCell Technologies) and stained as described above. Using a FACS Aria Fusion (BD Biosciences), B cells from SLE patients were sorted as: CD19⁺ CD11c⁻ CD27⁻ IgD⁺ naïve B cells, CD19⁺ CD11c⁻ CD27⁺ IgD⁻ memory B cells, CD19⁺ CD11c^hi IgD⁺ B cells or CD19⁺ CD11c^hi IgD⁻ B cells.

For data sets derived from Affymetrix platforms, GCRMA normalized expression values were variance corrected using local empirical Bayesian shrinkage before calculation of differential expression (DE) using the ebayes function in the open source BioConductor LIMMA package (31). Resulting p-values were adjusted for multiple hypothesis testing and filtered to retain DE probes with a False Discovery Rate (FDR) < 0.05 (32). For RNAseq datasets, FASTQC files were obtained from GEO. Trimmomatic was performed to cut adapter sequences, low quality reads, and the first six reads for each sequence because of non-random primer bias. Reads were aligned to human reference genome hg38 in STAR with default parameters, and then SAM files were converted to BAM files using Sambamba. Relative differential expression counts were generated using featureCounts. After careful examination with PCA, one outlier RA patient from GSE110999 was removed. Data were normalized using DeSeq2 bioconductor R package. Differentially expressed gene (DEG) comparisons between various cell populations were done using generalized linear modeling from the DeSeq2 package. Resulting p-values were adjusted for multiple hypothesis testing and filtered to retain DEG probes with an FDR < 0.05. FastQC, Trimmomatic, STAR, Sambamba, and the featureCounts programs are all free, open source programs available at the following web addresses:

Exploratory analysis like Principal Component Analysis (PCA) was conducted on the DeSeq2 normalized counts using R Bioconductor package factoextra (v1.0.7). The arguments addEllipses, ellipse.type, ellipse.level, were used to draw stat ellipses around each cell type at 95% confidence level in order to determine any outlier samples. The top 1,000 most differentially expressed genes were used for PCA analysis.

The database for annotation, visualization and integrated discovery (DAVID) (http://david.abcc.ncifcrf.gov/) was used to determine enriched gene ontology (GO) biological pathways (BP) for increased or decreased human or mouse gene symbols. BIG-C, a functional clustering tool that sorts increased and decreased genes into 52 categories based on their most likely biological function and/or cellular localization, was also employed. The BIG-C is based on information from multiple online tools and databases including UniProtKB/Swiss-Prot, GO Terms, MGI database, KEGG pathways, NCBI PubMed, and the Interactome. Each gene is placed into only one category based on its most likely function to eliminate the redundancy in enrichment, sometimes found in GO BP annotation (33).

GSVA (V1.25.0) was used as a non-parametric, unsupervised method for estimating the variation of pre-defined gene sets in samples of microarray expression data sets (www.bioconductor.org/packages/release/bioc/html/GSVA.html) (34). The inputs for the GSVA algorithm were a gene expression matrix of log2 expression values for pre-defined gene sets co-expressed in SLE datasets (Supplementary Data File 1). Enrichment scores (GSVA scores) were calculated non-parametrically using a Kolmogorov Smirnoff (KS)-like random walk statistic and a negative value for a particular sample and gene set, meaning that the gene set has a lower expression than the same gene set with a positive value. The enrichment scores (ES) were the largest positive and negative random walk deviations from zero, respectively, for a particular sample and gene set. The positive and negative ES for a particular gene set depend on the expression levels of the genes that form the pre-defined gene set. GSVA calculates enrichment scores using the log2 expression values for a group of genes in each SLE patient and healthy control and normalizes these scores between −1 (no enrichment) and +1 (enriched).

Morpheus (Morpheus, https://software.broadinstitute.org/morpheus) was used for unsupervised hierarchical clustering analysis. The inputs for Morpheus were the DEGs (log fold change) for each B cell population compared to either naïve or memory B cells. The distance measure was one minus Pearson's correlation coefficient with average linkage and clustering set to both rows and columns.

The Bioconductor package entitled homologene (v1.1.68, https://www.rdocumentation.org/packages/homologene/versions/1.1.68) was used for conversion of mouse gene symbols to their human ortholog gene symbols. Mouse gene symbols which did not convert to a human ortholog equivalent were entered manually into the Mouse Genomic Informatics database (informatics.jax.org) to check if there was a human homolog.

RNA from eight SLE patients (Supplementary Table 2) was isolated by resuspending cell pellets in 100 μL PBS and adding 1 mL TRIzol reagent (Thermofisher Scientific). 200 μL chloroform was added to each sample, mixed, and centrifuged at 12,000 rfc to separate RNA (aqueous phase) from DNA (interphase). RNA was precipitated with 2 μL polyacyl carrier and 500 μL isopropanol, followed by centrifugation and resuspension in 25 μL RNAase-free water. cDNA was generated using Superscript III enzyme (Thermofisher Scientific) per manufacturer's recommendation. IgH transcripts were amplified with Herculase II Fusion DNA Polymerase (Agilent) and leader (L) primers for each V_H family and primers for each C_H gene: 5′L-VH1; ACA GGT GCC CAC TCC CAG GTG CAG, 5′L-VH3; AAG GTG TCC AGT GTG ARG TGC AG, 5′L-VH4/6; CCC AGA TGG GTC CTG TCC CAG GTG CAG, and 5′L-VH5; CAA GGA GTC TGT TCC GAG GTG CAG (35) that were paired with constant gene primers: hIgM; GGC CAC GCT GCT CGT ATC C, hIgD; GAC CAC AGG GCT GTT ATC CTT TGG, hIgG; CGC CTG AGT TCC ACG ACA CC, and hIgA; GGA AGA AGC CCT GGA CCA GGC. PCR products were cloned into pSC-A-amp/kan cloning vector (Agilent) and Sanger sequenced. IgH VDJ sequence identity and mutational abundance were identified using IMGT database (http://www.imgt.org/). Clonal sequences were identified as having the same V- and J-gene segments with identical CDR3 nucleotide sequences. The use of family specific V-gene primers was chosen to allow for Sanger sequencing to minimize errors generated during other forms of IgH amplification. This technique will introduce minor amplification bias into the analysis; however, the amplification bias will be identical for all populations analyzed.

GraphPad PRISM 8 version 8.2.1 was used to perform mean, median, mode, standard deviation, and Tukey's multiple comparisons test, The Fisher's exact test was performed in R. For IgH repertoire analysis, two-tailed T-test was performed to examine V-gene usage difference.

RNA-seq and microarray datasets, with accession numbers GSE110999, GSE13917, GSE92387, GSE65928, GSE69033, GSE68878, and GSE28887, are available on the NCBI's database Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). The FCRL4⁺ tonsil B cell microarray data (formerly listed as CBX-40 at CIBEX (36) is available at ftp://ftp.ddbj.nig.ac.jp/ddbj_database/gea/experiment/E-GEAD-000/E-GEAD-336/. Supplementary Data File 2 lists gene expression dataset details including sorted populations and platforms.

Previous reports suggested there might be two subpopulations of CD11c⁺ cells in SLE patients based on IgD expression: IgD⁺ activated naïve cells that produce autoantibodies, and IgD⁻ double negative cells (CD27⁻ IgD⁻, DN2) that are poised to become plasma cells (3, 37, 38). To examine how different they are, we sorted CD11c^hi IgD⁺ and IgD⁻ cells and compared their transcriptional relatedness to naïve and conventional memory B cells. Using deposited RNAseq data (Supplementary Data File 2) from SLE patients (Supplementary Table 1) (1), principal component (PC) analysis revealed that CD11c^hi IgD⁺ and IgD⁻ may be separate populations, but they are vastly different from naive (CD11c⁻ CD27⁻ IgD⁺) and memory (CD11c⁻ CD27⁺ IgD⁻) B cells (Figure 1A). PC1 represented 47% of the variance between populations and separated naïve and memory cells from the IgD⁺ and IgD⁻ B cells, and PC2 represented 17% of the variance between populations and separated memory from naïve B cells. When DEGs from CD11c^hi IgD⁺ or IgD⁻ B cells were compared to naïve or memory B cells, the majority of genes were shared between the two subsets (Figure 1B, Supplementary Data File 3). Many DEGs encoding immune-related cell surface and signaling proteins were detected in both CD11c^hi IgD⁺ and IgD⁻ B cells in comparison to naive and/or memory B cells (Table 1), emphasizing both the uniqueness of the total CD11c^hi population and the similarity of the IgD⁺ and IgD⁻ subsets. Although only four samples were analyzed, the overall similarity of the results suggested that the populations were highly representative.

GO BP analysis of the differentially expressed transcripts compared to naïve or memory B cells, demonstrated that CD11c^hi IgD⁺ and IgD⁻ B cells had similarly enriched functional categories, in agreement with there being few differentially expressed transcripts between the two populations isolated based on IgD expression (Figure 1C, Supplementary Data File 4). A second functional clustering program, biologically informed gene clustering (BIG-C) (33), also demonstrated agreement with enriched categories for CD11c^hi IgD⁺ and IgD⁻ B cells compared to naïve and memory B cells, with notable differences of enrichment in autophagy in IgD⁺ B cells and endosome and vesicles in IgD⁻ B cells when compared to naïve B cells, but no significant differences when compared to memory B cells (Supplementary Figure 2A, Supplementary Data File 5). The functional overlap was further corroborated by the similarity in the expression of hundreds of pathways by CD11c^hi IgD⁺ and IgD⁻ B cells using Ingenuity Pathway Analysis (IPA) (Supplementary Figure 2B). Out of 410 pathways, a few cell cycle related pathways had a Z-score difference of more than 1 between IgD⁺ and IgD⁻ cells, with IgD⁻ having lower, negative scores, although not reaching Z score values of −1.8 for true significance. Additionally, mTor and complement signaling related pathways had a Z score difference of more than 1 between IgD⁺ and IgD⁻ cells; these scores were higher in IgD⁻ B cells and were significant.

Dendrogram grouping of transcripts showed intermingling of CD11c^hi IgD⁺ and IgD⁻ B cells and a similar distance of these populations from memory and naive B cells, suggesting patient-specific differences were more important than IgD expression by CD11c^hi cells, and supporting the overlap in DEGs when compared to naïve and memory B cells (Supplementary Figure 2C). Whereas, comparison of CD11c^hi IgD⁺ or IgD⁻ B cells to either naïve or memory B cells detected thousands of DEGs, there was a difference of only 175 DEGs when CD11c^hi IgD⁺ and IgD⁻ cells were compared to each other (Supplementary Data File 6). Pathway analysis using GO BP determined that the CD11c^hi IgD⁺ B cells had increased expression of genes associated with the respiratory burst and cell adhesion, whereas IgD⁻ B cells expressed increased transcripts for intracellular signal transduction, microtubule organization, and positive regulation of transcription (Supplementary Data File 7). Taken together, the CD11c^hi IgD⁺ and IgD⁻ B cells are highly similar and potentially consist of a single population with variable levels of IgH class switch recombination.

Gene set variation analysis (GSVA) was used to examine the co-expression of groups of genes identified as increased or decreased in CD11c^hi B cells compared to both naïve and memory B cells (Supplementary Data Files 3–5). A few germinal center (GC) B cell markers were uniformly increased in both CD11c^hi IgD⁺ and IgD⁻ B cells compared to naïve and memory B cells. Inhibitory protein-tyrosine phosphatases (PTPs), inhibitory dual-specific phosphatases (DUSPs), inhibitory Fc receptor-like (FCRLs), and inhibitory signaling adaptors were also significantly enriched in CD11c^hi IgD⁺ and IgD⁻ cells compared to naïve and memory B cells, as were genes associated with glycolysis. Likewise, gene signatures for lysosome, endocytosis, Ras superfamily of GTPases, cytoskeleton, and Golgi were uniformly enriched in CD11c^hi B cells compared to naïve and memory B cells. However, genes involved in mRNA translation were significantly decreased in both CD11c^hi IgD⁺ and IgD⁻ (Figure 2). Uniform enrichment of groups of inhibitory signaling molecules in CD11c^hi IgD⁺ and IgD⁻ B cells has not previously been described and may suggest an interrupted or inhibited signaling process.

Our transcriptional analysis (Table 1) showed that CD11c^hi IgD⁺ and IgD⁻ B cells from SLE patients have decreased expression of CR2 (CD21), CD27, and CXCR5 consistent with B cell populations from RA and combined variable immunodeficiency (CVID) described in the literature as anergic (7) (CD21⁻ CD27⁻), DN2 from SLE patients (3) (CD27⁻ IgD⁻ CXCR5⁻), atypical from malaria patients (21) (CD21⁻ CD27⁻), and FCRL4⁺ tissue memory B cells from tonsil (36). Additionally, GC B cell markers AICDA, BCL6, DAPP1, and RGS13 were increased in both CD11c^hi IgD⁺ and IgD⁻ cells compared to naïve and memory (Table 1), which suggested this population could be closely related to GC B cells. However, genes for plasma cell markers XBP1 and SLAMF7 were also increased in both IgD⁺ and IgD⁻ subsets compared to naïve and memory B cells, implying CD11c^hi B cells might display an intermediate phenotype with both GC and plasma cell markers.

To determine whether these B cells are a common cell population which arise under different stimuli, we compared gene expression profiles of the above populations with naïve, memory, centroblasts, centrocytes, plasmablasts from tonsil, and plasma cells from bone marrow (39). As can be seen in Figure 3, CD11c^hi B cell transcripts from SLE and RA patients clustered most closely with DN2 (SLE), DN2 (HD), anergic (RA/CVID), and atypical (malaria) B cells compared to either naive (Figure 3A) or memory (Figure 3B) B cells from each data set. Hierarchical clustering also demonstrated that the CD11c^hi cell subset shared more transcripts in common with either naïve or memory B cells, and fewer transcripts in common with GC centroblasts, centrocytes, and plasma cells, even though the CD11c^hi B cells expressed some common GC or plasma cell transcripts when compared to naïve or memory B cells. Of note, FCRL4⁺ tissue memory cells (tonsil) grouped most closely with centrocytes and centroblasts despite over-expressing ITGAX (CD11c). Compared to naïve or memory B cells, the DEGs in centroblasts and centrocytes were associated with DNA repair, chromatin remodeling, and cell cycle genes, but very few of these transcripts were increased in the CD11c^hi B cell subsets when compared to naïve and memory B cells. Similarly, compared to naïve and memory B cells, the DEGs in plasma cells were linked to Ig secretion, unfolded protein response, proteasome, endoplasmic reticulum and Golgi genes, but these gene were not expressed in the CD11c^hi B cell subsets compared to memory and naïve B cells (Supplementary Data File 8). Hierarchical clustering demonstrated that cell populations referred to as atypical (malaria), anergic (RA/CVID) and DN2 (SLE) were most closely related to the CD11c^hi B cell populations and were separated from all other B cell populations. Notably, the CD11c^hi-like populations all shared increased expression of ITGAX (CD11c), FCLR5, and TBX21 (T-bet) along with down regulation of CR2 (CD21) and CD27 compared to naïve (Figure 3C) and memory (Figure 3D). We have previously observed that CD11c^hi cells arise in healthy individuals, albeit at a lower frequency than during SLE (1). Analysis of the DN2 population from healthy donors showed close clustering with the CD11c^hi-like populations from various diseases. This suggested that the population can arise under a variety of circumstances and is not limited to SLE.

To understand how the CD11c^hi IgD⁻-like populations are related, GSVA was used to examine enrichment of gene sets that were detected in CD11c^hi B cells from SLE. Strikingly, B cell populations isolated on the basis of low IgD, CD27, CD21, and CXCR5 had similar GSVA enrichment patterns as CD11c^hi B cells. All five populations had increased GC marker expression and high expression of inhibitory PTPs, DUSPs, FCRLs and signaling adaptors compared to naïve and memory B cells. These individual B cell subsets expressed increased integrin, transporter, lysosome, endocytosis, RAS superfamily, cytoskeleton, and Golgi signatures similar to CD11c^hi B cells from SLE patients. Furthermore, genes representative of IL-21 signaling, previously demonstrated to be increased in CD11c^hi B cells from SLE patients (1), were increased in these B cell subsets as well. However, transcripts for mRNA translation and cell surface/signaling transcripts decreased in CD11c^hi IgD⁺ and IgD⁻ cells from SLE patients were not universally enriched in all the CD11c^hi-like B cells, compared to naïve and memory B cells suggesting that these differences may be specifically related to SLE (Figure 4, Supplementary Data File 9). Notably, all subsets had increased ITGAX, TBX21, SLAMF7, FCRL5 and SYK and decreased CCR7, CD24, CR2, CXCR4, CXCR5, and IL4R compared to naïve and memory B cells.

A population in mice called age-associated B cells (ABCs) is considered by some to be the equivalent of the human DN2 cell population, and the human and murine populations are considered singular in some review articles (40, 41). ABCs are B220⁺CD19⁺CD11b⁺ spleen cells that also may express CD11c, but they are not isolated in a similar manner to any of the human populations. To determine the relationship of mouse splenic ABCs to human peripheral blood CD11c^hi B cells, DE analysis of the T-bet⁺ (Tbx21) CD11b⁺ (Itgam) CD11c⁺ (Itgax) murine ABC subset compared to the mouse follicular (FO) B cell subset from murine spleen of C57BL/6 mice (GSE28887) (27) was carried out. Mouse genes were then converted to their human orthologs, and the DEGs were compared to those of human CD11c^hi B cells differentially expressed to naïve B cells because healthy mice have mostly naïve B cells. The ABC population shared 311 transcripts with human CD11c^hi IgD⁻ B cells, including increased Itgax, Tbx21, Syk, Cd72, Ptpn22, and Sox5, and decreased Cr2, Cxcr5, and Fcer2 (Figure 5A, Supplementary Data File 10), but had no increased expression of GC markers. Seven GO BP categories were in common with human CD11c^hi B cells: cilium morphogenesis, protein dephosphorylation, microtubule-based process, peptidyl-tyrosine autophosphorylation, receptor protein tyrosine kinase signaling, protein autophosphorylation, and protein phosphorylation. However, GO BP categories with the most significant enrichment in mouse ABCs, such as immune system process, innate immune signaling, integrin signaling, and regulation of phagocytosis, were not enriched in human CD11c^hi B cells (Figure 5B, Supplementary Data File 11). Hierarchical clustering analysis using human orthologs of murine genes that were increased in ABCs compared to murine FO B cells showed that mouse ABCs were equally distant from all human B cell subsets (Figure 5C, Supplementary Data File 12). A deeper analysis for potential inhibitors of B cell signaling in ABCs revealed an increase in Siglece, Ptpn22, Cd72, PirB, Pilra, Cd5 and Lair1, suggesting some similarity to the large number of inhibitory signaling genes increased in CD11c^hi-like B cells, but the mouse ABC subset also expressed an increase in multiple cytokines and chemokines, including Il18, Cxcl9, Ccl9, Ccl8, Ccl6, Cxcl10, Il6, Tnfsf13b, and Cxcl13 not detected as increased in any of the human CD11c^hi-like populations. The increased integrin pathways in mouse ABC could reflect their origin as tissue resident cells and may not be an appropriate comparator to human circulating cells. Overall, these analyses suggest that mouse ABCs are distinct from human CD11c^hi-like populations, and caution should be used in interpreting them as equivalent from a functional standpoint at this time.

We next compared the repertoires of expressed V_H genes of five populations of B cells that were sorted from eight SLE patients: CD11c^hi IgD⁺, CD11c^hi IgD⁻, CD11c⁻ CD27⁻ naïve, CD11c⁻ CD27⁺ IgD⁻ memory, and CD27⁺ CD38⁺⁺ plasmablasts/plasma cells (Supplementary Figure 1, Tables 2, 3). Sequencing revealed that CD11c^hi IgD⁺ and IgD⁻ populations from SLE patients were polyclonal and utilized only two genes at a significantly different frequency than the naïve population (Figure 6A), e.g., V3–30 in CD11c^hi IgD⁺ B cells and V3–13 in both CD11c^hi IgD⁺ and IgD⁻ B cells. The V4–34 gene associated with autoimmunity (37) was highly expressed in CD11c^hi B cells, but not at a level different from the naïve population. However, when compared to memory B cells, both CD11c^hi IgD⁺ and IgD⁻ populations manifested greater variance in V-gene utilization with 3 genes over-expressed and 3 genes under-expressed (Figure 6B). When compared to circulating plasma cells, both CD11c^hi subsets showed 6 genes were over-expressed, including V4–34, and 2 V genes were under-represented (Figure 6C). This suggests that the CD11c^hi population contains a diverse repertoire similar to naïve, and they appear to lack significant selection against antigen.

Mutation frequencies were measured to confirm active diversification of V_H genes. Naïve B cells had a frequency of 0.4 × 10⁻² mutations/bp, which likely represents the PCR error rate (Figure 7A). The average mutation frequency in CD11c^hi B cells was substantially elevated for IgD⁺ sequences, 2.8 × 10⁻², and IgD⁻ sequences, 4.7 × 10⁻². The higher mutation frequency in V_H genes from IgD⁻ cells was confirmed in samples from most of the individual patients (Figure 7B), which supports the conclusion that switched cells undergo more mutation events (42). As a comparison, the mutational frequency in memory and plasma cells was even higher at 6.0 × 10⁻² and 7.7 × 10⁻², respectively (Figure 7A). Nonetheless, the data are consistent with the observation that both CD11c^hi IgD⁺ and IgD⁻ B cells have mutational frequencies typical of a GC experience.

To further confirm that CD11c^hi IgD⁺ and IgD⁻ B cell subsets from SLE are a common cell population with various stages of class switching, we analyzed the clonal relationship between the two populations. To create clonally-related trees, DNA sequences with identical V_H CDR3 nucleotide sequences were considered to be related and were grouped together. In each case, the common precursor had undergone somatic hypermutation with additional mutations occurring after the cells branched. The results in Figure 8A show that unique mutations occurred in both IgD⁺ and IgD⁻ branches. Since our analysis suggests that CD11c^hi B cells are GC-emigrants, we wanted to address whether they originate from the same precursor cell as memory and plasma cells. Comparing CD11c^hi to memory B cells, we identified several clonally-related cells, which were independent of IgD status (Figure 8B). Similar clonal relationships were seen with CD11c^hi and plasma cells (Figure 8C). Taken together, this analysis demonstrates that CD11c^hi B cells originate from the same precursor population as other antigen-experienced B cells.