Lupus susceptibility results from the combined effects of numerous genetic loci, but the contribution of these loci to disease pathogenesis has been difficult to study due to the large cellular heterogeneity of the autoimmune immune response. We performed single-cell RNA, BCR, and TCR sequencing of splenocytes from mice with multiple polymorphic lupus susceptibility loci. We not only observed lymphocyte and myeloid expansion, but we also characterized changes in subset frequencies and gene expression, such as decreased CD8 and marginal zone B cells and increased Fcrl5- and Cd5l-expressing macrophages. Clonotypic analyses revealed expansion of B and CD4 clones, and TCR repertoires from lupus-prone mice were distinguishable by algorithmic specificity prediction and unsupervised machine learning classification. Myeloid differential gene expression, metabolism, and altered ligand–receptor interaction were associated with decreased Ag presentation. This dataset provides novel mechanistic insight into the pathophysiology of a spontaneous model of lupus, highlighting potential therapeutic targets for autoantibody-mediated disease.
This article is featured in Top Reads, p.
Systemic lupus erythematosus (SLE) is a chronic multiorgan autoimmune disease characterized by development of autoantibodies and potentially fatal lupus nephritis (1, 2). Current treatments for SLE are limited to broad immunosuppressants with little development of novel targeted therapies in recent decades, underpinned by an incomplete understanding of disease pathogenesis (3–5). Mouse models and human linkage studies have demonstrated that SLE is a complex polygenic disease linked to many individually weak susceptibility alleles that in combination with environmental factors result in a range of disease phenotypes (6, 7). Detailed genetic dissection of these loci has led to the identification of a growing list of susceptibility genes (8–10), but the relationship between susceptibility loci and cellular and molecular mechanisms of disease pathogenesis has not been fully characterized.
Spontaneous mouse models of lupus have provided a valuable tool to both discover and validate susceptibility loci in human disease (5, 9, 11). Congenic mouse studies led to the identification and characterization of the NZM2410-derived Sle lupus susceptibility loci consisting of Sle1, Sle2, and Sle3 (12, 13). Sle1 includes at least seven independent loci, including the Slam/Cd2 haplotype and polymorphic Fcgr2b that in combination lead to loss of tolerance to chromatin (14, 15). Addition of the BXSB-derived y-linked autoimmune accelerating (yaa) locus, which includes Tlr7 to create bicongenic B6.Sle1.yaa (SLE.yaa) mice, leads to the development of highly penetrant fatal glomerulonephritis (16–18). Individual potential therapeutic targets have been identified or validated in SLE.yaa mice, including IL-6 (19, 20), TLR7 (16), CXCR4 (21), IL-21 (22), and BANK1 (23) pathways, but broad characterization of the immune landscape of SLE.yaa mice has not been performed.
Generalized mechanistic studies in lupus-susceptible mice have been difficult due to the complexity of immune cell types and phenotypes that arise to drive disease pathogenesis (11, 24, 25). Therapeutic approaches for SLE have historically focused on the adaptive immune response (4, 26), and the potential role of innate immune cell types such as myeloid cells in disease pathogenesis has only recently been appreciated (19, 27–29). These studies have been unable to capture the heterogeneity of cell states and cell–cell interactions that is made possible by recent advances in single-cell sequencing technologies. Paired single-cell RNA sequencing (scRNA-seq) with BCR and TCR sequencing allows for characterization of innate and adaptive cell states at single-cell resolution. Unbiased approaches to profile the entire immune landscape of lupus-prone mice would not only facilitate discovery of SLE therapeutic targets but also provide biological insight of the autoreactive germinal center.
In this study, we performed scRNA-seq, single-cell BCR sequencing (scBCR-seq), and single-cell TCR sequencing (scTCR-seq) of splenocytes from SLE.yaa or immunized wild-type mice to generate a single-cell dataset of a spontaneous mouse model of lupus. Single-cell sequencing captured vast changes in immune cell heterogeneity and phenotype, including identification of novel myeloid subpopulations associated with lupus. Metabolic modeling revealed opposing changes in metabolic activity between B cells and myeloid cells in lupus-prone mice. BCR tracing identified clonal expansion toward the plasma cell (PC) compartment, and algorithmic analyses of the TCR repertoire were capable of separating CD4 clonotypes based on genotype. Ligand–receptor interaction analysis highlighted signaling relationships between innate and adaptive immune cell types that may contribute to loss of tolerance.
Functional validation using coculture experiments demonstrated decreased ability of autoimmune splenocytes to stimulate Ag-specific CD4 responses. As the first single-cell atlas of lupus-prone mice to our knowledge, our dataset serves as a resource for therapeutic target discovery and biological characterization of autoantibody-mediated disease.
Materials and Methods
C57BL/6J (B6, JAX no. 000664), B6.SJL (CD45.1, JAX no. 002014), and SLE.yaa mice (JAX no. 021569) were obtained from The Jackson Laboratory. OT-II mice were a gift from Gabriel Victora (Rockefeller University) and crossed onto a CD45.1 background. Mice were maintained at Harvard Medical School in an American Association for the Accreditation of Laboratory Animal Care–accredited facility. Fifteen-week-old male mice were used for all experiments unless noted otherwise. All animal experiments were approved by the Institutional Animal Care and Use Committee of Harvard Medical School (IS111).
Immunization was performed by i.p. injection of 100 μg of 4-hydroxy-3-nitrophenyl acetyl hapten conjugated to OVA (NP-OVA, Biosearch Technologies) in 50 μl of HBSS emulsified with 50 μl of Imject alum adjuvant (Thermo Scientific). Three weeks after immunization, mice were boosted with i.p. 100 μg NP-OVA in 100 μl of HBSS (Corning Life Sciences). Mice were sacrificed 10 d after booster immunization.
Spleens were mechanically digested through a 70-μm cell strainer (Corning Life Sciences), incubated in RBC lysis buffer (155 mM NH4Cl, 12 mM NaHCO3, 0.1 mM EDTA) for 3 min at room temperature, and washed with FACS buffer (PBS with 1% heat-inactivated FBS, 1 mM EDTA, and 0.05% sodium azide). Cells were counted and 1 × 106 cells/well were incubated with 50 μl of appropriate Abs and eFluor 780 viability dye (eBioscience) in FACS buffer for 30 min on ice. The following Abs were used: anti-CD45.1 (A20, 1:400), anti-CD45.2 (104, 1:400), anti-CXCR3 (CXCR3-173, 1:200), anti-CD1d (1B1, 1:300), anti-CD24 (M1/69, 1:200), anti-CD19 (1D3/CD19, 1:200), anti-CD21/35 (7E9, 1:300), anti-CD138 (281-2, 1:200), anti-CCR6 (29-2L17, 1:200), anti-CD74 (In1/CD74, 1:200), anti-CD16/32 (2.4G2, 1:300), anti-F4/80 (BM8, 1:300), anti-IA/IE (M5/114.15.2, 1:2000), anti-CD31 (390, 1:300), anti-CD11c (N418, 1:300), anti-CD68 (FA-11, 1:300), anti-CD11b (M1/70, 1:300), anti-CD267, anti-B220 (RA3-6B2, 1:400), anti-CD268, anti-Ly6c (HK1.4, 1:300), anti-Ly6d (49-H4, 1:300), anti-Ly49C/F/I/H (14B11, 1:300), anti-CD122 (TM-β1, 1:300), anti-CD8 (53-6.7, 1:300), anti-CD16.2 (9E9, 1:300), anti-Fcrl5 (FAB6757, 1:300), anti-CD5L (AF2834, 1:300), anti-MARCO (FAB2956A, 1:300), anti–Ki-67 (SolA15, 1:500), anti-CD153 (RM153, 1:200), anti-ICOS (C398.4A, 1:200), anti–Lag-3 (C9B7W, 1:200), anti-CD49b (DX5, 1:300), anti-CD3 (500A2, 1:400), anti-CD44 (IM7, 1:1000), anti-CD69 (H1.2F3, 1:300), anti-CD62L (MEL-14, 1:200), anti-CXCR5 (L138D7, 1:200), anti-PD1 (RMP1-30, 1:200), anti–Sca-1 (D7, 1:400), anti-GL7 (GL7, 1:300), anti-CD4 (GK1.5, 1:400), and anti-Foxp3 (FJK-16s, 1:200). For two-step staining procedures cells were incubated with 50 μl of streptavidin-conjugated fluorophore or Alexa Fluor 488 anti-goat (Thermo Fisher Scientific, A27012, 1:400) in FACS buffer for 20 min on ice. For intracellular staining, cells were fixed with fixation/permeabilization buffer (eBioscience) for 30 min at room temperature, washed with permeabilization buffer (eBioscience), and incubated with 50 μl of intracellular Ab in permeabilization buffer for 30 min at room temperature. Cells were resuspended in FACS buffer and read on a FACSCanto II (BD Biosciences) with 488-, 405-, and 640-nm lasers using FACSDiva (BD Biosciences). Data were analyzed using FlowJo (Tree Star).
Immunofluorescence confocal microscopy
Spleens were perfused with PBS followed by 2% paraformaldehyde (Electron Microscopy Sciences) in PBS. Tissues were fixed in 2% paraformaldehyde for 8 h at 4°C, cryoprotected with 30% sucrose in PBS overnight at 4°C, perfused with 30% OCT (Tissue-Tek) in PBS, and embedded in 100% OCT in standard Cryomolds (Tissue-Tek) in the vapor phase of liquid nitrogen and stored at −80°C. Frozen sections were cut on a cryostat at a thickness of 20 μm and allowed to dry for 60 min at room temperature. Sections were fixed with acetone for 5 min at −20°C, then permeabilized and blocked with 10% normal rat serum (Thermo Fisher Scientific) in immunofluorescence buffer (PBS with 0.2% BSA and 0.3% Triton X-100) for 1 h at room temperature. Slides were stained with primary Ab in immunofluorescence buffer overnight at 4°C and secondary Ab or streptavidin-conjugated fluorophore in immunofluorescence buffer for 4 h at room temperature where appropriate. The following Abs were used: PacBlue anti-CD21/35 (7E9, 1:100), goat anti-CD5L (AF2834, 1:200), biotin anti-CD11c (N418, 1:100), allophycocyanin anti-CD11b (M1/70, 1:100), Alexa Fluor 488 anti-CD74 (In1/CD74, 1:100), FITC anti-CD169, PE anti-CD11b (M1/70, 1:100), Alexa Fluor 647 anti-CD16/32 (2.4G2, 1:100), PacBlue anti-GL7 (GL7, 1:100), Alexa Fluor 488 anti-CD16.2 (9E9, 1:100), and Alexa Fluor 488 anti-goat (Thermo Fisher Scientific, A27012, 1:400). Slides were mounted using Fluoro-Gel (Electron Microscopy Sciences). Images were acquired using a FluoView FV3000 confocal laser scanning microscope (Olympus) with ×10 or ×30 objectives and analyzed using Fiji (ImageJ).
Splenocyte and CD4 coculture
Spleens were harvested into ice-cold MACS buffer (PBS with 0.5% BSA and 2 mM EDTA) and mechanically digested through a 70-μm cell strainer (Corning Life Sciences). Spleens were incubated in RBC lysis buffer (155 mM NH4Cl, 12 mM NaHCO3, 0.1 mM EDTA) for 3 min at room temperature and washed with MACS buffer. OT-II CD45.1 splenocytes were pooled and CD4 cells were isolated using the CD4 T cell isolation kit (Miltenyi Biotec) and labeled with 5 μM CellTrace Violet (Thermo Fisher Scientific) for 20 min at 37°C. Total splenocytes from B6 or SLE.yaa mice were counted and 4 × 105 cells/well were added to round-bottom 96-well plates and incubated with indicated concentrations of NP-OVA (Biosearch Technologies) in complete RPMI (RPMI 1640 with 10% heat-inactivated FCS, 50 μM 2-ME, 2 mM l-glutamine, 100 U/ml penicillin, 100 μg/ml streptomycin, 20 mM HEPES, 1 mM sodium pyruvate, and 100 μM nonessential amino acids) for 1 h at 37°C. Splenocytes were washed with complete RPMI, and 1 × 105 OT-II CD45.1 CD4 cells/well were added to splenocytes or 5 μg/ml anti-CD3 (145-2C11, BioLegend)– and 5 μg/ml anti-CD28 (37.51, BioLegend)–coated wells with 10 ng/ml human IL-2 (PeproTech). After 4 d supernatants were stored at −80°C and cells were subject to flow cytometry as described above.
Droplet-based scRNA-seq, scBCR-seq, and scTCR-seq
The scRNA-seq, scBCR-seq, and scTCR-seq libraries were prepared using the 10x Genomics single-cell immune profiling solution kit. CD45 cells were isolated from 107 splenocytes per sample using CD45 MicroBeads (Miltenyi Biotec) according to the manufacturer’s protocol. Biological replicates were kept separate. Immediately after isolation cells were captured in gel beads-in-emulsions at a targeted recovery of 10,000 cells per sample and multiplet rate of ∼7.6% using a chromium controller (10x Genomics), followed by barcoding, gel beads-in-emulsions cleanup, and cDNA amplification for 13 cycles. For gene expression library preparation, 50 ng of amplified cDNA was fragmented and end repaired, double-sided size selected with SPRIselect beads (Beckman Coulter), adaptor ligated, PCR amplified with indexing primers, and double-sided size selected with SPRIselect beads. For VDJ library construction, BCR and TCR transcripts were enriched from 2 μl of amplified cDNA by PCR, and 50 ng of PCR product was fragmented and end repaired, size selected with SPRIselect beads, adaptor ligated, PCR amplified with indexing primers, and size selected with SPRIselect beads. Sequencing of scRNA, scBCR, and scTCR libraries were performed together on a NovaSeq 6000 (Illumina) to a minimum sequencing depth of 40,000 reads per cell.
Processing and filtering of scRNA-seq data
The Cell Ranger (10x Genomics, version 5.0.0) count pipeline was used to align 5′ gene expression reads to the GRCm38 reference genome (mm10-2020-A). We obtained reads from 30,885 cells with an average of 1,639 genes per cell and 67,601 reads per cell. Individual sample matrices were loaded in Seurat (version 4.0.0) (30) using the Read10x function and filtered for cells with at least 200 genes detected and genes detected in at least three cells using the CreateSeuratObject function, leaving 14,190 cells from B6 mice and 16,326 cells from SLE.yaa mice. Individual samples were merged using the merge function, and S and G2/M cell cycle phase scoring was assigned using CellCycleScoring. To remove batch effects between samples associated with a heat shock gene expression signature, genes annotated with the Gene Ontology biological process (GOBP) term “cellular response to heat” (GO: 0034605) was used to assign a heat shock score using AddModuleScore. Cells with <200 or >3,000 genes detected, <1,000 reads detected, >5% mitochondrial RNA content, >25% rRNA content, an S phase score >0.15, or a G2/M phase score >0.15 were excluded from analysis, with 11,194 cells from B6 mice and 11,091 cells from SLE.yaa mice passing the filters. BCR and TCR variable and constant genes were excluded from scRNA-seq analysis to prevent clustering based on VDJ transcripts.
Unsupervised clustering of scRNA-seq data
Regularized negative binomial regression was performed on cells from B6 or SLE.yaa mice separately using the sctransform normalization method (31) to normalize, scale, and select variable genes, and to regress out mitochondrial RNA content, rRNA content, the number of unique molecular identifiers, and heat shock score. B6 and mixed SLE.yaa datasets were then integrated (32) using SelectIntegrationFeatures, PrepSCTIntegration, FindIntegrationAnchors, and IntegrateData. Following principal component analysis (PCA), clusters were identified using FindClusters to apply shared nearest neighbor–based clustering using the first 25 principal components with k = 20 and resolution = 0.15. The same principal components were used to generate uniform manifold approximation and projection (UMAP) projections.
Diffusion map and pseudotime analysis
Seurat objects were exported to Scanpy (version 1.5.1) using anndata2ri (version 1.0.2). Partition-based graph abstraction was performed using the PAGA (partition-based graph abstraction) function (33) with 15 neighbors and the first 20 principal components. A randomly selected follicular B cell was used as the root cell for diffusion pseudotime computation using the first 10 diffusion components. Diffusion component coordinates and pseudotime values were added back to the Seurat object using CreateDimReducObject and AddMetaData, respectively.
RNA velocity analysis
RNA velocity processing was performed on cellranger counts using velocyto (version 0.17.17) (34), and RNA velocity was computed from loom files using scVelo (version 0.2.3) (35). Seurat objects were exported using anndata2ri (version 1.0.2), and RNA velocity embeddings were added and visualized using scVelo.
Cell cluster annotation
Clusters were annotated based on expression of marker genes for known populations, including Cd19 (B cells), Cd4 (CD4 T cells), Cd8a (CD8 T cells), Cd68 (macrophages), S100a8 (polymorphonuclear cells [PMNs]), Itgax (dendritic cells [DCs]), Klra8 (NK cells), Ifit3 (B cell IFN-stimulated gene [B-ISG]), Siglech (plasmacytoid DCs [pDCs]), Mki67, Xbp1, and Jchain (PCs), Cd63 (basophils), and Vpreb3 and Iglc1 (B-early cells). Cluster determination was confirmed by identifying differentially expressed marker genes for each cluster using FindAllMarkers with the MAST (model-based analysis of single-cell transcriptomics) algorithm (36) and comparing to known cell type–specific marker genes. Cluster names were updated using RenameIdents. B cells were subset using the B, B-early, PC, and B-ISG clusters with the subset function, leaving 7243 cells from B6 mice and 5943 cells from SLE.yaa mice. Normalization, integration, and clustering were performed as above. Clusters were annotated based on expression of marker genes for known populations, including Sell, Cd55, and Ighd (follicular B cells), Mif, Cd83, and Myc (activated B cells), Fcrl5, Cr2, Cd1d1, Sema7a, and Dtx1 (marginal zone [MZ] B cells), Lars2 (transitional stage 2 B [T2-B] cells), Jchain, Xbpi1, Prdm1, and Slamf7 (PCs), Cd24a, Cd93, Vpreb3, Ms4a1, and Niban3 (transitional stage 1 B [T1-B] cells), Ifit3 and Isg15 (B-ISG), Apoe, Spn, and Itgb1 (B1 cells). T cells were subset using the CD4 and CD8 clusters with the subset function, leaving 2648 cells from B6 mice and 1814 cells from SLE.yaa mice. Clusters were annotated based on expression of marker genes for known and novel populations, including Cd4 and Cd40lg (CD4 T cells), Cd8a (CD8 T cells), Ccl5 and Nkg7 (CD8 effector [CD8-eff] cells), Foxp3, Ctla4, and Ikzf2 (regulatory T [TREG] cells), Cd74 (CD74), and Cxcr6 (CD4 effector [CD4-eff] cells). Myeloid cells were subset using the macrophage, PMN, DC, and pDC clusters with the subset function, leaving 1067 cells from B6 mice and 3238 cells from SLE.yaa mice. Clusters were annotated based on expression of marker genes for known and novel populations, including Cd68, Apoc2, and Ace (Mac-1 cells), Ccr2, Lyz, Crip1, and Fn1 (monocytes), S100a8, Slpi, and Csf3r (PMN-1 cells), Cd68, Itgax, Cd72, C1qb, Fcrl5, and Fcgr4 (Mac-2 cells), S100a8, Ngp, and Camp (PMN-2 cells), Cd68, Marco, Mertk, Sirpa, Adgre1, Vcam1, and Cd5l (Mac-3 cells), Cd79a and Cr2 (follicular DCs [FDCs]), Siglech (pDCs), and Cd63 and Cyp11a1 (basophils). CD4 cells were subset using the CD4, TREG, and CD4-eff cell clusters with the subset function, leaving 1544 cells from B6 mice and 1159 cells from SLE.yaa mice. Clusters were annotated based on expression of marker genes for known populations, including Foxp3 and Il2ra (TREG cells), Ccl5, Stat1, and Isg15 (CD4-eff cells), Smg1 (CD4 mem cells), Pdcd1 and Id2 (CD4-exhausted cells), Ccr7 and Il7r (CD4-naive cells), Nr4a1 (CD4-activated cells), Pdcd1, Sostdc1, Icos, Bcl6, and Cxcr5 (TFH cells). CD8 cells were subset using the CD8 and CD8-eff clusters with the subset function, leaving 935 cells from B6 mice and 562 cells from SLE.yaa mice. Clusters were annotated based on expression of marker genes for known populations, including Ccr7 and Sell (CD8-naive cells), Ccl5, Gzmk, and S100a6 (CD8-eff cells), and Slamf6 and Id3 (CD8 stem cells).
Differential gene expression
Differentially expressed genes (DEGs) between SLE.yaa and B6 mice were determined using FindMarkers with the MAST algorithm (36) across all genes. For differential expression analysis within individual clusters or clonotypes, Seurat objects were first subset using the subset function. To compare differential gene expression between different models of autoimmunity, DEGs were computed using scRNA-seq data from follicular CD4 cells isolated from mixed autoimmune chimeras (37), and log fold changes were compared on a gene-by-gene basis to log fold changes observed in SLE.yaa mice.
Analysis of human scRNA-seq data
Reference scRNA-seq data from human renal biopsies were obtained from ImmPort (SDY997) (38) and processed as previously described (37). Briefly, raw count matrices from healthy controls and SLE patients were imported into Seurat and processed as described above. Clusters were annotated based on expression of marker genes for known populations, including CD3E, CD4 (CD4 T cells), CD8A (CD8 T cells), LYZ (macrophages), KLRF1 (NK cells), MS4A1, BCL11A (B cells), and IRF8 (monocytes). Comparison between human and mouse differential expression analysis was performed by determining the homologs of each gene using getLDS from biomaRt (version 2.24.0). From the 15,817 DEGs identified by scRNA-seq of mouse cells and 13,558 DEGs identified by scRNA-seq of human renal biopsies, 10,854 homologs were identified in both datasets.
Gene expression signature scoring
Pathway analysis and gene set enrichment analysis were performed using clusterProfiler (version 3.14.3) (39). DEGs were selected using p-adj <0.01 and ranked according to log2 fold change (log2FC) for enrichment analysis. Ranked gene lists were used to query GOBP (40, 41) and MSigDB (version 7.0.1) (42, 43) signature libraries. For Gene Ontology analysis and annotation, DEGs were selected using p-adj <0.05 and absolute log2FC >0.1 thresholds.
Data processing of scBCR-seq and scTCR-seq libraries
The Cell Ranger (10x Genomics, version 5.0.0) vdj pipeline was used to align BCR and TCR reads to the vdj-GRCm38-alts-ensemble 5.0.0 reference genome (10x Genomics). After BCR alignment we obtained reads from 15,270 cells with an average of 71,917 reads per cell. Only BCRs with full-length and productive H and L chain sequences were included in analysis. After TCR alignment we obtained reads from 5,063 cells with an average of 52,352 reads per cell. Only TCRs with full-length and productive α- and β-chain sequences were included in the analysis.
BCR clonality analysis
BCR repertoire and clonality analysis was performed using the Immcantation suite (version 4.1.0) (44, 45). Clonal clustering and germline reconstruction were performed from Cell Ranger contigs using Change-O (version 1.0.0). Clonal distances were calculated using SCOPer (version 1.1.0), and clonal lineage trees were created using IgPhyML (version 1.1.3). Lineage reconstruction and topology, VDJ gene usage, and repertoire diversity were calculated using Alakazam (version 1.0.2). Mutation profiling, selection pressure scores, and somatic hypermutation modeling were performed using shazam (version 1.0.2). Clonality was matched with gene expression analysis in Seurat by adding clonality information to the metadata using AddMetaData based on cell barcodes. A UMAP clone size topography map was created using clonalOverlay, clonal networks were created using clonalNetwork, and clonal Morisita overlap scores between different clusters were calculated using clonalOverlap in scRepertoire (version 1.1.4) (46).
TCR clonality analysis
TCR clonality analyses were performed as previously described (37). Clonotypes were defined by identical productive CDR3α and CDR3β amino acid sequences, and clone size was calculated by the number of unique cell barcodes belonging to an individual clonotype. Clonality was matched with gene expression analysis in Seurat by adding clonality information to the metadata using AddMetaData based on cell barcodes. Unweighted TCR network analysis between samples and conditions was performed using the qgraph package (version 1.6.9). To evaluate public clonotype environments, non-metric multidimensional scaling was performed using vegan (version 2.5-7) with k = 2. Stress <0.1 was confirmed using a Shepard plot. Public repertoires were compared using immunarch (version 0.6.7).
Unsupervised classification of TCR repertoires
TCR repertoires of each individual chimera were constructed from CDR3αβ and VDJ gene usage from scTCR-seq. TCR featurization was performed using a variational autoencoder (VAE) in DeepTCR (version 1.4.15) (47) as previously described (37). We used a neural network with 256 latent dimensions, k = 5 for the first convolutional layer of the graph, learned latent dimensionality of 64 for amino acids, learned latent dimensionality of 48 for VDJ genes, latent α of 0.001, and three convolutional layers with 32, 64, and 128 neurons, respectively. The VAE was trained using an Adam optimizer with the learning rate = 0.001 until the convergence criterion of a >0.01 decrease in the determined interval was met. For sample-agnostic clustering, a heatmap was created to assess unbiased hierarchical clustering of the sample genotype based on VAE featurization.
Metabolic heterogeneity modeling of single cells was performed using COMPASS (version 0.9.9.6.2) (48). Differential metabolic activity was computed from a COMPASS metareaction consistency scores subset by scRNA-seq–based cluster annotation. Metareactions were filtered for core reactions, defined as reactions with Recon2 confidence of either 0 or 4 and annotated with an enzyme commission number. PCA was performed on COMPASS metareaction space using factoextra (version 1.0.7). Cell cluster annotation and metadata were loaded from Seurat, and the top two principal components were used for two-dimensional mapping. Signature scores were assigned to individual cells using AddModuleScore and gene lists from KEGG (49), GOBP (40, 41) and MSigDB (version 7.0.1) (42, 43). Spearman correlations between signature scores and each individual metareaction or principal component were computed using all individual cells. Correlation coefficients for different signature scores were compared on a metareaction-by-metareaction basis.
Ligand–receptor interaction analysis was performed using CellPhoneDB (version 2.1.7) (50). Mouse gene were converted to human orthologs using getLDS from biomaRt (version 2.24.0). Separate count matrices for each genotype were exported from Seurat, and ligand–receptor interaction scores and p values were computed using CellPhoneDB.
TCR specificity group prediction
The grouping of lymphocyte interaction by a paratope hotspots (GLIPH2) algorithm (51) was used to predict TCR specificity groups as previously described (37). The GLIPH2 mouse TCR dataset was used for reference. A Fisher’s exact test was used to assess the statistical significance of a given motif, and specificity groups were filtered for clusters with significant V gene bias (p < 0.05 by GLIPH2) and significant final score (p < 1 × 10−5 by GLIPH2). Specificity group prediction was matched with gene expression analysis in Seurat using AddMetaData based on clonotype sequences.
CDR3β database construction
To predict Ag specificities, we used our previously published (37) reference database of known CDR3β sequences and their cognate Ags. To test whether TCR sequences identified by scTCR-seq had known Ag specifities, we searched for the presence of each CDR3β sequence in our reference database. To extend Ag predictions to CDR3β sequences not found in the reference database, we performed GLIPH2 analysis on our reference database, identifying annotated Ags associated with CDR3β local motif patterns. Ag prediction was matched with gene expression analysis in Seurat using AddMetaData based on GLIPH2 local motif predictions.
Flow cytometry quantification was compared using an unpaired two-tailed Student t test. COMPASS metabolic scores were compared using Wilcoxon rank sums for each metareaction, and Cohen’s d statistic was used to assess effect size (48). Metareaction correlation with gene signature scores was computed using Spearman’s rank correlation coefficient. Differential gene expression was computed using the MAST algorithm (36). Pathway analysis was performed using DEGs and a Fisher’s exact test. The p values were adjusted for multiple comparisons using a Bonferroni correction.
Bar graphs were created using ggpubr (version 0.4.0.999) or Prism (GraphPad Software), and correlations versus correlation scatter plots were created using ggplot2 (version 3.3.3). Biological theme comparisons, network plots, and gene set enrichment plots were generated using clusterProfiler (version 3.18.1). Heatmaps and hierarchical clustering was performed using pheatmap (version 1.0.12). UMAP and violin plots comparing gene expression across samples and clusters were generated using Seurat (version 4.0.0). Sequence motifs were created using ggseqlogo (version 0.1). Scatter plots, stacked bar graphs, and dot plots were created using ggplot2. Volcano plots were generated using EnhancedVolcano (version 1.8.0), and DEGs were highlighted in red.
All scRNA-seq, scBCR-seq, and scTCR-seq data generated in this study have been deposited in the GEO database and are available under accession number GSE192762 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE192762).
Relevant codes are available through GitHub (https://github.com/egarren/scSLE).
Lupus-prone mice have expanded splenic lymphocyte and myeloid populations
SLE.yaa mice develop spontaneous severe glomerulonephritis with detectable autoantibodies by 4–6 mo of age (16, 18). To profile the entire immune landscape of mice undergoing active loss of tolerance, we performed scRNA-seq, scBCR-seq, and scTCR-seq on CD45 cells isolated from the spleens of two 15-wk-old SLE.yaa mice. Fifteen weeks represents a time point immediately preceding pathologic disease onset in SLE.yaa mice (18), and therefore differences in single-cell profiles are more likely to represent pathogenic drivers of autoimmunity than consequences of established autoimmunity. Two age- and sex-matched B6 mice immunized with NP-OVA were used as controls to represent an active immune response against a foreign Ag, facilitating identification of autoimmune-associated immunophenotypes rather than generalized immune activation. Known marker genes (Fig. 1A) and DEGs (Supplemental Fig. 1A) were used to annotate the 12 clusters identified by scRNA-seq. All major splenic cell types were recovered, including B cells, CD4 cells, macrophages, CD8 cells, PMNs, DCs, PCs, NK cells, and pDCs (Fig. 1B). The relative frequencies of individual clusters varied between immunized B6 and SLE.yaa mice, with a decrease in the relative frequency of B and CD8 cells and an increase in the frequency of macrophages, PMNs, DCs, and PCs in lupus-prone mice (Fig. 1C). Flow cytometric quantification confirmed expansion of these lymphocyte and myeloid populations in the spleens of lupus-prone mice (Fig. 1D). SLE.yaa mice exhibit profound splenomegaly, which accounts for significant increases in absolute total lymphocyte counts despite relative decreases in individual subset frequency.
Recent immunometabolic studies have led to greater appreciation for a regulatory role of metabolism in various disease contexts, including autoimmune disease (52–57). To compare the metabolic states of individual cells in lupus-prone mice to immunized mice, we applied the metabolic modeling algorithm COMPASS (48) to our single-cell dataset. PCA of the COMPASS metareaction space was unable to separate cells by transcriptome-defined clusters (Fig. 2A), but principal component 1 (PC1) displayed a modest ability to separate SLE.yaa from B6 splenocytes (Fig. 2B). By combining COMPASS metareaction principal components with gene expression data, we determined that PC1 correlated with cell glycolysis (Spearman’s rho = 0.477, p < 2.2 × 10−16) and oxidative phosphorylation (Spearman’s rho = 0.514, p < 2.2 × 10−16) gene signature scores (Fig. 2C), suggesting that lymphocytes from lupus-prone mice have increased glycolytic flux. Oxidative phosphorylation is paradoxically increased in T cells from SLE patients (57) and can drive TH17 pathogenicity over TREG cell development (58). To better characterize the consequences of increased oxidative phosphorylation signature scores in our dataset, we compared the correlations between individual metareactions with other signatures scores, identifying negative relationships with BCR signaling (Spearman’s rho = −0.762, p < 2.2 × 10−16) and the complement cascade (Spearman’s rho = −0.710, p < 2.2 × 10−16), and a positive relationship with Ag presentation (Spearman’s rho = 0.751, p < 2.2 × 10−16) (Supplemental Fig. 1B).
COMPASS metareactions provide greater resolution of the metabolic network in individual cells than general pathway analysis (48). To identify metabolic pathways and reactions associated with lupus in individual cell types, we compared differential metabolic activity on a metareaction-by-metareaction basis in each individual cluster. B and CD4 cells from lupus-prone mice displayed increased activity in most metareactions, whereas CD8 cells had a similar number of upregulated and downregulated metareactions, and macrophages had decreased activity in most metareactions (Fig. 2D). The increase in metabolic activity in B and CD4 cells was so profound that nearly all metareactions assigned to individual metabolic pathways had increased activity, as was the case for decreased activity in macrophages (Fig. 2E). In contrast, individual metabolic pathways in CD8 cells included both upregulated and downregulated metareactions, suggesting more nuanced changes in metabolism than a broad pathway level change. For example, the transaldolase step of the pentose phosphate pathway was upregulated in CD8 cells from lupus-prone mice (Cohen’s d = 0.228, p-adj = 3.23 × 10−5), whereas G6PD initiation of the pentose phosphate pathway was downregulated (Cohen’s d = −0.164, p-adj = 0.009). The pentose phosphate pathway, specifically the rate-limiting G6PD step, influences CD8 activation and polarization (59, 60). These findings suggest that although B and CD4 cells have increased general metabolic activity in lupus-prone mice, macrophages from lupus-prone mice are less metabolically active.
Human SLE is a heterogenous disease linked to many loci beyond those homologous to Sle1 and yaa (7, 9, 61). To determine whether the differential expression observed in lupus-prone mouse splenocytes is also present in human disease, we compared our dataset with scRNA-seq performed on renal biopsies from SLE patients (38, 62). Unbiased clustering of CD45 cells from SLE patients or healthy controls identified seven distinct clusters, all of which were also present in our mouse dataset (Supplemental Fig. 1C). Few B cells were recovered from healthy controls (Supplemental Fig. 1D), making differential gene expression analysis difficult for this cluster. We were able to perform differential gene expression analysis on CD4 cells, CD8 cells, and macrophages between SLE and healthy control biopsies (Supplemental Fig. 1E), identifying decreased IL7R in SLE CD4s, decreased SOD2 in SLE CD8s, and decreased CXCR4 in SLE macrophages, among other DEGs. To identify DEGs present in both datasets, we compared fold changes of all genes in CD4 cell, CD8 cell, and macrophage clusters (Supplemental Fig. 1F). CD4 cells in both mouse and human disease upregulated Tigit, Maf, Ikf2, Itgb1, and Ahnak, although the coinhibitory molecules Lag3 and Pdcd1 were only upregulated in CD4 cells from lupus-prone mice. CD8 cells in both mouse and human disease upregulated Nkg7, Itgb1, Eomes, and Cx3cr1, and macrophages in both mouse and human disease downregulated Fcgr2b. These results demonstrate that despite differences in both tissue and species, our mouse dataset captures some gene expression changes observed in human disease.
B cells are clonally expanded with altered subsets in lupus-prone mice
B cells are activated and produce autoreactive IgG in SLE.yaa mice (16); however, deeper characterization of B cell subsets and clonal evolution in lupus-prone mice has not been performed. To this end, we subset B cell clusters from our dataset, and after additional clustering identified eight clusters of B cells present in both B6 and SLE.yaa mice (Fig. 3A). As observed by flow cytometry (Fig. 1D), PCs were markedly increased in SLE.yaa mice relative to B6 mice, accompanied by a compensatory decrease in the relative frequency of follicular and MZ B cells (Fig. 3B). This transition might represent the follicular entry, CD4 engagement, and eventual Ab secretion of MZ B cells observed in lupus-prone mice (63–65). As we were able to capture distinct states of B cell development among splenic B cells, we next performed pseudotime (66) (Fig. 3C) and RNA velocity (34) (Fig. 3D) analyses to compare their potential developmental trajectories. T2-B cells differentiated into either follicular or MZ B cells (Fig. 3D), a fate decision dependent on both BCR signaling strength and BAFF availability (67). Both PC and B-ISG clusters represented later pseudotime states (Fig. 3C), with trajectories initiating from both follicular and MZ B cell clusters by RNA velocity and PAGA graph (Fig. 3E) analyses. PCs originate from both germinal center and extrafollicular compartments in SLE, although their origin likely influences their autoreactivity and pathogenicity (68–70). These observations confirm the emergence of an Ab-secreting cell population in SLE.yaa mice and provide greater resolution to the accompanying changes in B cell states in lupus-prone mice.
Beyond changes in relative abundances of B cell subsets, autoimmunity is likely driven by changes in gene expression and phenotypes of individual B cell subsets. We performed differential gene expression analysis to compare phenotypes of follicular, activated, MZ, or T2 B cells from SLE.yaa and B6 mice (Fig. 3F). Follicular B cells from lupus-prone mice increased expression of Ly6a, Serpina3g, Slamf6, Psmb9, Cxcr3, Stat1, Ccr6, Tbx21, Itgax, and Pltp and decreased expression of Dnaja1, Cd24a, Cr2, Tnfrsf13b, and Tnfrsf13c; activated B cells increased expression of Ly6a, Gimap4, Psmb9, Stat1, Cxcr3, Ccr6, and Serpina3g and decreased expression of Cr2, Cd24a, Tnfrsf13c, and Fcmr; MZ B cells increased expression of Blvrb, Gimap4, Iigp1, Jun, Socs1, Stat1, Vim, Cd83, and Ccr6 and decreased expression of Cr2, Dph5, Myof, S1pr3, and Marcks; and T2-B cells increased expression of Jchain, Ly6a, Tap1, Ly6c2, and Cxcr3 and decreased expression of Cd24a, Dnaja1, Cd74, and Tnfrsf13c, among other DEGs. Cd83, Gimpa4, Iigp1, Jun, Ly6a, Ly6c2, Pltp, Psmb9, Serpina3g, Socs1, Stat1, and Tap1 are ISGs (71), likely reflecting the excessive IFN response in SLE.yaa mice (72). Tmsb4x expression was increased in most B cell clusters, although this gene is expressed on the yaa duplication (17) and therefore likely represents an artifact of our model. Slamf6 is a transmembrane protein that signals via SAP (signaling lymphocyte activation molecule–associated protein) to promote CD74/MIF (macrophage migration inhibitory factor) prosurvival signaling between B and CD4 cells (73) and therefore represents an appealing potential pathway for further investigation. Increased expression of Tbx21 and Itgax and decreased expression of Cr2 are potentially reflective of the emergence of a subcluster of CD21−CD11c+T-bet+ age-associated B cells, which have been described in both human and murine SLE (74). Collectively, these DEGs represent candidate drivers or consequences of SLE in specific B cell subsets not previously appreciated at single-cell resolution.
To investigate the biological relevance of DEGs in B cells from lupus-prone mice, we performed Gene Ontology and pathway analyses. Biological theme comparison associated DEGs from follicular B cells in lupus-prone mice with endoplasmic reticulum protein processing, Ag processing, Th17 differentiation, and BCR signaling pathways (Fig. 3G). Network visualization demonstrated that these associations were driven by decreased Cr2 and Dnaja1 and increased Stat1 and Ifi30 expression (Fig. 3H), likely due to both increased IFN signaling and a transition to Ab-secreting phenotypes. Gene Ontology analysis of DEGs in activated B cells from lupus-prone mice revealed similar pathway activation, namely protein processing, Th17 differentiation, BCR signaling, and Ag processing (Fig. 3G). B cell Ag presentation is a necessary driver of autoimmunity in lupus-prone mice and can facilitate B cell escape of tolerance (75). These findings suggest that both follicular and activated B cells experience alterations in Ag presentation and protein-processing pathways in SLE.
To confirm the B cell state and gene expression differences observed in our dataset, we performed flow cytometry on spleens of 15-wk-old SLE.yaa mice and age-matched immunized B6 controls. We confirmed that SLE.yaa mice have an ∼4-fold decrease in MZ B cells (Fig. 3I). Differential gene expression was validated at the protein level in total B cells, germinal center cells, and PCs (Fig. 3J), including decreased expression of CD267 (also known as transmembrane activator and CAML interactor [TACI], protein name for Tnfrsf13b), decreased expression of CD268 (also known as BAFF receptor, protein name for Tnfrsf13c), increased expression of stem cell ag-1 (Sca-1, protein name for Ly6a), decreased expression of CD74, decreased expression of CD21/35 (protein name for Cr2), decreased CD24, increased CCR6, and increased CXCR3. CD267 and CD268 recognize BAFF and a proliferation-inducing ligand (APRIL), which are drivers and validated therapeutic targets for SLE (76, 77), and therefore their downregulation likely represents negative feedback from excessive BAFF and APRIL signaling. CD74 deficiency ameliorates disease severity in lupus-prone mice and is an experimental target in human disease (78, 79). CD21/35 deficiency increases susceptibility to SLE and is correlated with disease severity (80, 81). Polymorphic CD24 is associated with SLE, although its pathogenic mechanism is unknown (82). CCR6 expression correlates with SLE severity and through interactions with its cognate ligand CCL20 regulates germinal center kinetics (83, 84). CXCR3 is upregulated on memory B cells in response to IFN signaling (85), although its relevance to SLE has not been studied. These observations provide validation of our transcriptional dataset at the protein level.
We next paired our scRNA-seq analysis with scBCR-seq on B cells from lupus-prone mice. Clonal topography analysis indicated a shift in clonal abundance from MZ, T1-B, and T2-B cells in B6 mice toward the PC compartment in SLE.yaa mice, accompanied by an increase in mutational frequency in SLE.yaa PCs (Fig. 4A). Decreased clonality within the MZ compartment supports the hypothesis of MZ B cell entry into the PC compartment, particularly given the increased cross-reactivity of MZ BCRs and decreased activation threshold relative to follicular B cells (86, 87). Clonal overlap analysis revealed the greatest Morisita overlap index between the follicular and activated B cell clusters, and the least overlap between the T1-B and B1 or T1-B and PC clusters (Fig. 4B), supporting the proposed intercluster trajectories modeled by RNA velocity (Fig. 3D). Clonal network analysis was consistent with overlap score calculations (Fig. 4C), identifying the greatest amount of MZ B cell clonotype sharing with PCs. As expected, SLE.yaa B cells had markedly decreased diversity compared with B6 B cells (Fig. 4D) and increased somatic hypermutation mutability (Fig. 4E), likely representing clonal expansion, affinity maturation, and increased plasma cells. These clonal analyses represent a markedly more robust humoral response in lupus-prone mice compared with immunized mice.
T cells are clonally expanded and express markers of exhaustion in lupus-prone mice
The earliest studies of SLE.yaa mice described changes in CD4 molecular signatures toward a TFH cell phenotype with altered cytokine production (18). In this study, we use our dataset to extend these findings toward increased pathways, genes, and T cell states with both single-cell transcriptomic and clonal resolution. We subset T cell clusters from our dataset, and after additional clustering identified six clusters of T cells present in both B6 and SLE.yaa mice (Fig. 5A). Interestingly, we identified a novel T cell cluster that was defined by expression of Cd3e, Cd74, and Cd79a (Supplemental Fig. 2A) whose frequency was decreased in SLE.yaa mice (Fig. 5B). CD74 is involved in Ag presentation by MHC class II with additional roles in prosurvival signaling (88, 89), and CD79a has been identified in T cell lymphoma (90, 91), although the functional relevance of this cluster remains unclear. Among the other clusters, CD8 cells were relatively decreased and TREG cells were relatively increased in SLE.yaa mice (Fig. 5B), consistent with our flow cytometry immune profiling (Fig. 1D). We subset and clustered CD4 cells, identifying eight different cell states including TREG and TFH cells (Supplemental Fig. 2B). Naive CD4 cells were decreased whereas TFH, exhausted, and activated CD4 cells were increased among CD4 cells from lupus-prone mice (Supplemental Fig. 2C). Similarly, CD8 cells were further subset into naive, effector, and stem-like CD8 cells (Supplemental Fig. 2E), with a relative increase in effector and compensatory decrease in naive CD8 cells among CD8 cells from lupus-prone mice (Supplemental Fig. 2F). There were too few cells in these CD4 and CD8 subsets to perform differential expression analysis. These findings describe changes in T cell subset frequencies, most notably decreased CD8 cells and increased exhausted CD4 cells, that would likely not be appreciable by bulk sequencing.
To study gene expression changes within individual T cell subsets, we performed differential gene expression analysis between SLE.yaa and B6 mice (Fig. 5C). CD4 cells from lupus-prone mice increased expression of Ly6a, Itgb1, Lag3, Srgn, Maf, Eea1, Pdcd1, Cxcr3, Tnfsf8, and Tigit; CD8 cells increased expression of Ly6a, Stat1, Ptms, AW112010, Cxcr3, Pdcd1, and Lag3; effector-like CD8 cells increased expression of Nkg7, Gzmk, Id2, and Ly6a; and TREG cells increased expression of S100a11, Lag3, Tigit, and Srgn, among other DEGs. AW112010, Id2, Itgb1, Lag3, Ly6a, Maf, Nkg7, Pdcd1, Ptms, Srgn, and Stat1 are ISGs (70, 71, 92), and Pdcd1, Tigit, and Lag3 are coinhibitory molecules that are upregulated in exhausted T cells (93). To further study the biological relevance of these DEGs, we performed gene set enrichment analysis, which identified significant associations between DEGs in CD4 and CD8-eff cells with Ag response and respective cell states in CD8 and TREG cells (Fig. 5D). To determine whether the transcriptional changes observed in CD4 cells from lupus-prone mice are recapitulated in another model of autoimmune disease (94), we compared our dataset with scRNA-seq and scTCR-seq data from follicular T cells isolated from mixed 564Igi autoimmune chimeras (37). CD4 cells from SLE.yaa mice and follicular T cells from mixed autoimmune chimeras both upregulated Ccl5, Ly6a, Lag3, and Tnfsf8, although Ahnak, Itgb1, and Ikzf2 were only upregulated in CD4s from SLE.yaa mice (Supplemental Fig. 2H).
We validated these cell state and gene expression changes by flow cytometry, including increased CD4+CXCR5+PD1+ follicular T cells, increased CD44+PD1+ activated CD4 and CD8s, and increased CD44+CD62L− effector memory CD4 and CD8 cells with compensatory decreases in CD44−CD62L+ naive CD4 and CD8 cells in 15-wk-old SLE.yaa mice (Fig. 5E). Differential gene expression was validated at the protein level in conventional T (TCON; CD4+Foxp3−), TREG (CD4+Foxp3+), TFH (CD4+CXCR5+PD1+Foxp3−ICOS+), T follicular regulatory (TFR; CD4+CXCR5+PD1+Foxp3+ICOS+), extrafollicular (CD4+CD62L−PSGL1−), and CD8 cells (Fig. 5F), including increased expression of ICOS, CD153 (protein name for Tnfsf8), CXCR3, Sca-1, PD1 (protein name for Pdcd1), and Lag-3. CD153+ follicular T cells are sufficient to induce spontaneous germinal centers via osteopontin secretion (95), whereas CXCR3 expression on CD4 cells regulates migration toward the IFN-inducible ligands CXCL9, CXCL10, and CXCL11 (96). Flow cytometry confirmed a relative increase in TREG cells in SLE.yaa mice, and Foxp3 expression was not increased in extrafollicular CD4 cells or follicular T cells (Fig. 5F). We also identified decreased Ki-67 expression in TFH cells, but increased Ki-67 expression in TCON and TREG cells from SLE.yaa mice (Fig. 5F). These gene and protein expression changes suggest that both CD4 and CD8 cells in lupus-prone mice experience increased activation and accompanying exhaustion, with notable proliferation differences between TFH and TREG cells.
In addition to T cell phenotypic differences, the TCR repertoires of SLE patients are distinct (97). To examine the effects of lupus susceptibility loci on T cell clonality, we also performed scTCR-seq on T cells from lupus-prone mice. After filtering for productive and paired TCRα- and TCRβ-chains, we retained 1810 cells representing 1501 unique clonotypes from SLE.yaa mice and 2191 cells representing 2163 unique clonotypes from B6 controls. When compared with scTCR-seq data from follicular T cells isolated from mixed 564Igi autoimmune chimeras, we did not identify any shared TCR clonotypes, except for a single clone that was identified in both a SLE.yaa mouse and a control non-autoimmune chimera (Supplemental Fig. 2I). SLE.yaa mice had a greater number of expanded clones, with clonal expansion observed in CD4, TREG, and effector CD8 cells (Fig. 6B–D, Supplemental Fig. 3A), as observed in humans (98, 99). To compare TCR repertoires from SLE.yaa and B6 mice, we performed both repertoire-wide and clonotype-level analyses. Unbiased hierarchical clustering of variable gene usage separated SLE.yaa from B6 TCRs (Supplemental Fig. 3B), although CDR3 consensus sequences were not significantly different between SLE.yaa and B6 mice (Supplemental Fig. 3C). Machine learning classification of TCR repertoires using variable autoencoder featurization (47) was similarly able to distinguish SLE.yaa and B6 repertoires (Supplemental Fig. 3D), suggesting that SLE.yaa and B6 mice have fundamentally different TCR repertoires. Individual samples had limited T cell clonal overlap (Fig. 6E), but public clonotype size was capable of separating SLE.yaa and B6 mice (Supplemental Fig. 3E). Of the 3663 unique clonotypes identified, only two were shared among SLE.yaa mice and none was shared among B6 mice. One of these clonotypes (CAASDYGSSGNKLIF + CASSFSSQNTLYF) was shared among both SLE.yaa and B6 mice. These results suggest that SLE.yaa T cells are clonally expanded and distinguishable at both the repertoire and individual clonotype level.
To study whether T cell clonal expansion is associated with differential gene expression in lupus-prone mice, we paired scRNA-seq with scTCR-seq data. We compared DEGs in expanded compared with nonexpanded T cells between SLE.yaa and B6 mice (Fig. 6F). SLE.yaa expanded T cells expressed increased Eea1, Pdcd1, Maf, and Ikzf2 and decreased Trat1 relative to B6 expanded T cells. These findings suggest that the increased markers of exhaustion we observed in SLE.yaa T cells are not trivially due to increased activation and clonal expansion, but rather to intrinsic differences in clonal phenotypes. Further correlational comparisons using fold changes revealed that Rora expression was increased in expanded B6 T cells but decreased in SLE.yaa expanded T cells (Fig. 6G). Rora is a nuclear receptor involved in CD4 polarization (100, 101), suggesting that SLE.yaa T cells have altered polarization while undergoing clonal expansion. These analyses complement our DEG analysis by identifying additional candidate genes that are negatively associated with T cell clonal expansion in SLE.yaa mice.
Due to the cross-reactive nature of the TCR, clonotype level analysis might not be sensitive enough to capture Ag-specific differences in repertoire (102, 103). We therefore applied the GLIPH2 algorithm (51) to our scTCR-seq dataset to predict Ag specificity groups of SLE.yaa T cells. After filtering, GLIPH2 yielded 199 specificity groups representing 531 unique clonotypes. Integrating these specificity groups with our scRNA-seq dataset allowed us to predict the specificity of 187 cells (10.3%) from SLE.yaa mice and 203 cells (7.7%) from B6 mice. Nonmetric multidimensional scaling analysis (Supplemental Fig. 3F) and unweighted network analysis (Supplemental Fig. 3G) separated SLE.yaa and B6 mice based on specificity group sharing between mice, suggesting that SLE.yaa and B6 T cells have distinct specificities. Shared specificity group size between SLE.yaa and B6 mice also showed limited specificity group overlap between genotypes, although the largest specificity groups were shared between SLE.yaa and B6 mice (Fig. 6H). Unweighted network analysis demonstrated that specificity groups with clonotypes from both SLE.yaa and B6 mice were polyclonal (Supplemental Fig. 3H), confirming that GLIPH2 can capture predicted Ag overlap not appreciable from CDR3 sequences alone. These results suggest that although Ag specificity predication can separate the TCR repertoires of SLE.yaa and B6 mice, there is a subset of expanded clonotypes that share predicted Ag specificity between SLE.yaa and B6 mice.
To predict the Ag specificity of GLIPH2 groups, we annotated these groups using GLIPH2 analysis of a database of CDR3β sequences with known Ag specificities (37). Database annotation could predict Ag specificities of 41 of 199 (20.6%) GLIPH2 specificity groups. Successful specificity group annotation was limited to specificity groups uniquely present in either SLE.yaa or B6 mice, with few of the common specificity groups represented in the known Ag database (Fig. 6H). Notably, an autoimmune-related specificity group was identified among both SLE.yaa and B6 enriched specificity groups (Fig. 6H). Superimposition of UMAP visualization with disease-related Ag predictions revealed similarly broad distributions of specificity groups among SLE.yaa and B6 T cells (Fig. 6I). These data suggest that both SLE.yaa and B6 T cells have distinct specificities for viral- and autoimmune-related Ags, although we were unable to predict the specificities of SLE.yaa and B6 T cells predicted to react to similar Ags.
Myeloid subset distribution is altered in lupus-prone mice
The earliest studies of SLE.yaa mice described an expanded myeloid compartment (16) likely due to IL-6–mediated myelopoiesis (19), but the heterogeneity of this compartment has not been profiled. We used our dataset to subcluster splenic myeloid cells from SLE.yaa mice, identifying nine clusters of myeloid cells (Fig. 7A). Macrophages, determined by Cd68 and Adgre1 expression, were represented by three novel subclusters, which we named Mac-1, Mac-2, and Mac-3. Mac-1 cells were defined by Apoc2 and Ace, Mac-2 cells were defined by Cd72, C1qb, Fcrl5, and Fcgr4, and Mac-3 were defined by Vcam1 and Cd5l (Fig. 7B). Fcrl5 is an IgG receptor found on dysfunctional atypical memory B cells (104), but it has not been previously described on myeloid subsets, whereas CD5L is a glycoprotein secreted by macrophages that can regulate lipid metabolism and inflammation (105–107). We also identified two subclusters of PMNs: PMN-1 was defined by expression of Slpi, Csf3r, Il1b, and Cxcr2, whereas PMN-2 was defined by Ngp, Camp, and Ltf (Fig. 7B). Myeloid clusters were prominently altered in SLE.yaa mice, most notably by the emergence of Mac-2 and Mac-3 cells, which were nearly absent from B6 mice (Fig. 7C). The relative frequency of Mac-1 cells was also increased in SLE.yaa mice, whereas monocytes and PMN-1 were decreased. In addition to describing shifting myeloid subset distribution, these subset level changes identify novel myeloid clusters that arise in lupus-prone mice.
To study gene expression changes within individual myeloid subsets, we performed differential gene expression analysis between SLE.yaa and B6 mice (Fig. 7D). Mac-1 cells from lupus-prone mice increased expression of Hp, Txn1, Fabp4, and Pecam1 and decreased expression of Cd74, H2.Eb1, and H2.Aa; monocytes increased expression of Gpx1, Hp, Apoc2, and Actg1 and decreased expression of H2.DMa, Fcgr3, Irf1, and Slamf7; PMN-1 cells increased expression of Ifi27l2a, Marcksl1, C3, Lcn2, and Slfn4 and decreased expression of Fox, Pbx1, D8Ertd738e, and Tgbi; and Mac-2 cells increased expression of AW112010, S100a11, and Srgn and decreased expression of Rps24 and Rps7 among other DEGs. Too few cells were present in Mac-3, FDC, pDC, and basophil clusters to perform differential gene expression analysis. Apoc2, AW1102010, C3, Fabp4, Hp, Ifi27l2a, Lcn2, Macksl1, Pecam1, Slfn4, and Srgn are ISGs (71). Biological theme comparison associated DEGs from Mac-1 cells in lupus-prone mice with phagocytosis, immune effector process, and Ag processing and presentation, among other pathways (Fig. 7E). Monocyte DEGs were associated with bacterial responses, cell–cell adhesion, phagocytosis, and inflammatory response to antigenic stimulus, PMN-1 DEGs were associated with viral responses, and Mac-2 DEGs were associated with ribosome biogenesis and assembly (Fig. 7E). These findings suggest that myeloid cells from lupus-prone mice are broadly inflammatory and have altered Ag presentation pathway activity, with notable subset level differences in IFN response and ribosome processing.
We next used flow cytometry to validate these myeloid cell clusters and gene expression changes. To broadly identify splenic myeloid cells, we used dual staining for CD11c (protein name for Itgax) and CD11b (protein name for Itgam). We identified two distinct populations that were double positive for CD11c and CD11b (Fig. 8A), which we named DP1 (CD11chiCD11b+) and DP2 (CD11cmedCD11b+). As previously described (16, 19) and confirmed by scRNA-seq (Fig. 7C), myeloid cells were expanded in SLE.yaa mice, with increased relative frequency of DP1, DP2, and CD11c populations (Fig. 8B). Within the DP1 population, we also identified the emergence of a Fcrl5+ and CD16.2+ (protein name for Fcgr4) population in SLE.yaa mice likely representing our Mac-2 cluster, and a CD5L+ population, likely representing our Mac-3 cluster (Fig. 8A). Differential gene expression was validated in macrophages (F4/80+CD68+) and DP1, DP2, and CD11c cells (Fig. 8B), including increased expression of CD31 (protein name for Pecam1), decreased expression of CD16/32 (protein name for Fcgr3), and decreased expression of CD74 in SLE.yaa mice. As observed by scRNA-seq, not all protein level changes were significantly altered in all myeloid subsets, suggesting subset level differences in myeloid phenotypes. Longitudinal flow cytometric profiling revealed that these myeloid subpopulations emerged between 12 and 15 wk of age in SLE.yaa mice and persisted up to 21 wk (Supplemental Fig. 4A), suggesting that Mac-2 and Mac-3 represent stable states that emerge prior to development of pathologic autoimmunity.
To visualize the spatial distribution of myeloid subsets observed by scRNA-seq, we performed immunofluorescence of SLE.yaa spleens for markers defined by scRNA-seq. Follicular architecture was disrupted in SLE.yaa mice, with loss of germinal center organization and CD169+ MZ macrophages (Supplemental Fig. 4B). CD11c, CD11b, and double-positive myeloid cells were no longer excluded to interfollicular regions, but they were instead dispersed throughout the spleen, including within white pulp (Supplemental Fig. 4B). Increased expression of CD5L and CD16.2 and decreased expression of CD74 and CD16/32 among myeloid cells was validated by confocal microscopy (Fig. 8C). Myeloid subsets expressing CD5L or CD16.2 were widely distributed in SLE.yaa spleens, located both near and far from germinal centers (Supplemental Fig. 4B). CD74 was also expressed by B cells, which also downregulated its expression in SLE.yaa spleens (Fig. 8C), confirming B cell scRNA-seq observations (Fig. 3F). FDCs were not efficiently captured by scRNA-seq (Fig. 7A) but were easily visualized by immunofluorescence using the complement receptor CD21/35 (Supplemental Fig. 4B). In contrast to our myeloid observations, we identified robust expression of CD5L on FDCs in immunized B6 spleens, which was lost in SLE.yaa spleens (Fig. 8C). CD5L expression has previously been described in macrophages and lymphocytes (108), but, to our knowledge, this is the first observation on FDCs in situ. These histologic findings not only validate our scRNA-seq observations, but they also lend spatial insight toward follicular disorganization in lupus-prone mice.
Myeloid cells from lupus-prone mice are less capable of presenting Ag
Given the profound myeloid subset alterations and potentially disrupted Ag processing pathways in lupus-prone mice, we next sought to identify the functional consequences of the myeloid gene expression changes we observed by scRNA-seq. Confirming the downregulation of Ag presentation–related genes (H2.Aa, H2.Eb1, H2.Ab1, H2.DMa, Cd74) in Mac-1, Mac-2, and Mac-3 cells from SLE.yaa mice, we observed decreased I-A/I-E expression in macrophages, DP1 cells, and CD11c cells by flow cytometry (Fig. 9A). MHC class II expression was most decreased in DP1 cells (Fig. 9A), which include subsets representative of the Mac-2 and Mac-3 clusters identified by scRNA-seq (Fig. 8A). These data suggest that the novel Fcrl5- and Cd5l-expressing macrophage populations that emerge in SLE.yaa mice present less Ag via MHC class II.
To test whether lower I-A/I-E expression leads to decreased ability to activate Ag-specific CD4 responses, we performed an in vitro coculture experiment with Ag-loaded splenocytes and CD4 cells. Splenocytes from SLE.yaa mice loaded with OVA were less capable of stimulating OT-II CD4 T cell proliferation in a concentration-dependent manner (Fig. 9B). In the absence of Ag, coculture with SLE.yaa splenocytes led to a small but significant increase in CD4 proliferation, albeit to a much smaller degree than Ag or polyclonal TCR stimulation with anti-CD3 and anti-CD28 (Fig. 9B). Similar decreases in CD4 responsiveness to SLE.yaa splenocyte Ag presentation were observed based on CD44, CD69, ICOS, and PD1 expression (Fig. 9B). These results suggest that although SLE.yaa splenocytes might provide Ag-independent stimulation to CD4 cells, they have reduced ability to process and present Ag to stimulate cognate CD4 responses.
To identify cellular interactions beyond Ag presentation that might contribute to loss of tolerance, we applied CellPhoneDB (50) to our dataset. Interaction scores were calculated for known receptor–ligand pairs between all clusters identified by scRNA-seq. Interaction scores were largely similar between SLE.yaa and B6 mice, with the largest differences observed in the CEACAM1/CD209 and CSF1/SIRPA pathways (Fig. 10A). CEACAM1 is an adhesion molecule with isoform- and posttranslational modification–dependent ability to modulate inflammation (109), whereas CD209, also known as DC-SIGN, is a C-type lectin that mediates DC rolling and phagocyte activation (110). DC-SIGN binding to CEACAM1 through Lewis(x) moieties promotes neutrophil–DC interactions that modulate T cell responses (111). CSF1 is a growth factor that promotes macrophage differentiation (112), whereas SIRPA is an inhibitory receptor that prevents phagocytosis (113). The difference in interaction scores in these pathways emphasizes the prominence of altered myeloid signaling in SLE.yaa mice.
To gain insight into broad differences in cellular interactions, the total of all interaction scores for each cluster–cluster combination was compared between genotypes. SLE.yaa mice had increased reaction scores between most cell types, with the greatest increases observed between DCs and CD4 cells and between macrophages and B cells (Fig. 10B). To compare these altered interactions at the individual receptor–ligand level, we visualized all significant curated interactions (Supplemental Fig. 5). These analyses highlighted specific interactions of biological relevance altered in SLE.yaa mice, such as persistence of CD40/CD40LG interactions between CD4 and B cells but not between CD4 cells and DCs in SLE.yaa mice (Fig. 10C). CECAM1/CD209 signaling between DCs and B cells also decreased in SLE.yaa mice, whereas CSF1/SIRPA signaling between DCs and B-ISG and CD4 cells increased (Fig. 10C). These findings highlight altered cellular interactions, namely between innate and adaptive immune cells, that potentially contribute to the myeloid dysfunction observed in SLE.yaa mice.
In this study, we present a transcriptomic and clonal dataset of lupus-prone mice at single-cell resolution. Paired scRNA-seq, scBCR-seq, and scTCR-seq of splenocytes from SLE.yaa mice identified cluster, transcriptomic, and clonal differences among myeloid and lymphocyte populations in autoimmune disease. Whereas both preclinical and therapeutic studies of SLE have focused on adaptive immune cells (4, 26), our findings highlight a disrupted and dysfunctional myeloid compartment as a potential therapeutic target in SLE. We not only confirm previously described myeloid expansion in SLE.yaa mice (16, 18, 19, 114), but we also identify two novel myeloid subpopulations defined by Cd5l and Fcrl5 expression that emerge in SLE.yaa mice. CD5L signaling can promote an anti-inflammatory cytokine response to TLR activation and M2 macrophage polarization (105, 106), suppress DAMP-mediated inflammation (115), and restrain pathogenic Th17 polarization (107), although a pathogenic role in autoantibody disease has not been described. Our observation of paradoxical loss of CD5L expression on FDCs in lupus-prone mice suggests cell type–dependent regulation in autoimmune disease.
Our newly identified myeloid clusters also exhibited decreased Ag presentation–related genes, broadly suppressed metabolism, and decreased inflammatory receptor–ligand interactions, supporting the hypothesis that these myeloid subsets are less inflammatory despite being expanded in SLE.yaa mice. This is consistent with our observation of reduced Ag-specific CD4 activation by SLE.yaa splenocytes, which supports previous reports of myeloid dysfunction in SLE patients (116–120). Although prior studies attribute this to serologic factors present in SLE (117, 121, 122), our coculture experiments are consistent with myeloid cell–intrinsic dysfunction, although they are limited by use of a heterogenous APC population. Alternatively, reduced Ag-specific CD4 activation in our coculture assay might be due to differences in APC coinhibitory receptor expression or frequency in SLE.yaa mice. We hypothesize that emergence and follicular entry of a suppressive myeloid subset in SLE.yaa mice are a homeostatic response to excessive chronic inflammation and accumulation of cellular debris, explaining the decreased adaptive immune responses to foreign pathogens or immunization observed in SLE patients (123–127). Although detailed functional characterization and human validation of these myeloid subsets are necessary, these results introduce a novel therapeutic approach for SLE, with particular emphasis on myeloid-directed therapy.
The cell type agnostic nature of our dataset allowed us to assess a range of potential mechanisms driving SLE.yaa pathology beyond the myeloid compartment. In contrast to previous reports of MZ B cell expansion in B6.Sle1.Sle2.Sle3 mice (65), we observed MZ B cell contraction and upregulation of CCR6 in SLE.yaa mice. RNA velocity and clonotypic analyses suggested that this is due to MZ B cell conversion to PCs, supporting previous reports of activation of extrafollicular autoreactive MZ B cells toward Ab-secreting cells (63–65). Given their lowered activation threshold, polyreactive BCRs, and the ability to rapidly differentiate to PCs without T cell help (69, 86, 87), we hypothesize that MZ B cell differentiation to PCs represents a major source of autoantibodies in SLE.yaa mice. We also observed both previously described and novel phenotypic changes in other B cell compartments, including broadly increased glycolysis and oxidative phosphorylation, increased expression of Slamf6, Cxcr3, Tbx21, Itgax, and Ccr6, and decreased expression of Cr2, Cd74, Cd24, Tnfrsf13b, and Tnfrsf13c at the gene and protein levels. Recapitulation of previous findings, including therapeutic targets under active preclinical investigation such as CD74 (78, 79) or recently U.S. Food and Drug Administration–approved approaches such as BAFF targeting with belimumab (76, 77), supports the therapeutic potential of our newly identified potential targets for SLE.
CD8 abnormalities are observed in SLE patients and postulated to contribute to autoimmunity (128), including cytolytic deficiency of peripheral CD8 cells (129–131), expansion of effector memory tissue-infiltrating CD8 cells (132–134), and impaired suppressive activity of regulatory CD8 cells (135–139). Consistent with these findings, we observed a relative increase and clonal expansion of effector but not other CD8 subsets in lupus-prone mice. However, in contrast to circulating effector CD8s (129), these splenic effector CD8s appeared to have increased cytotoxicity due to increased Eomes and Gzmk expression. This supports more recent scRNA-seq observations from SLE patient kidney biopsies in which tissue-resident CD8 cells increased GZMB and GZMK expression (38), highlighting the modulatory ability of tissue residence and potential for in situ therapeutic selectivity. Detailed single-cell metabolic modeling revealed that in contrast to CD4 cells, which broadly increased metabolic activity, CD8 cells exhibited nuanced tuning of individual metabolic reactions. Oxidative stress, impaired mitochondrial respiration, and mTORC1 signaling contribute to CD8 exhaustion and dysfunction in SLE (140–145), supporting our observation of elevated NAD metabolism, reactive oxygen species detoxification, and fatty acid oxidation in CD8 cells. Given the CD8-specific changes in gene expression and metabolism in our dataset, we hypothesize that subtle metabolic rewiring influenced by tissue residence drives CD8 dysfunction and conversion to effector phenotypes, representing a unique therapeutic approach for SLE.
CD4 cells were also transcriptionally and clonally distinct in SLE.yaa mice, consistent with prior observations in other mouse models of lupus and human disease (37, 146–148). Metabolic profiling identified increased CD4 glycolysis and oxidative phosphorylation, a previously validated therapeutic target for murine lupus (149). Although TREG cells are dysfunctional in NZB × NZW F1 and MRL/lpr mouse models of lupus (150–152), we observed paradoxically increased TREG cell and decreased TFH cell proliferation in SLE.yaa mice, which more closely resembles the heterogenous TREG cell phenotypes of human SLE (Refs. 153–156 and T. Sasaki, S. Bracero, J. Keegan, L. Chen, Y. Cao, E. Stevens, Y. Qu, G. Wang, J. Nguyen, S. E. Alves, et al., manuscript posted on bioRxiv, DOI: 10.1101/2021.11.08.467791). CD4 cells also upregulated CXCR3 and CD153, which regulate CD4 trafficking (96) and germinal center initiation (95), respectfully, representing potential therapeutic targets for SLE. Despite our ability to distinguish SLE.yaa from wild-type TCR repertoires in silico, we also observed convergence of algorithmically predicted TCR specificity between autoimmune and non-autoimmune mice, consistent with our previous profiling of follicular T cells in a bone marrow chimera model of autoantibody disease (37). We hypothesize that autoimmune and non-autoimmune TCR specificity convergence reflects TCR cross-reactivity and is responsible for the ability of peripheral CD4 cells to provide help to autoreactive B cells despite central tolerance mechanisms (157–160). Reflecting our prior hypothesis of MZ B cell follicular entry and PC conversion, we predict that polyreactive MZ B cells with a lower activation threshold help polarize cross-reactive CD4 cells toward self. This model is consistent with molecular mimicry as a potential etiology of SLE, in which TCRs reactive to both foreign and self Ags may trigger autoreactive B cell responses following infection (161, 162).
This dataset, however, is subject to several limitations. In this study we profile splenocytes from SLE.yaa mice, which might not necessarily capture the immune cell phenotypes of human disease. Immune cell profiling of the blood or inflamed tissues such as the kidney, alternative mouse models of lupus (particularly those in which female mice may be used), different stages and severity of disease, and human biopsies are likely to yield both complementary and unique observations. Indeed, changes in splenic immune cell populations or phenotypes might reflect leukocyte migration or emigration from the spleen to other secondary lymphoid organs or peripheral tissues (163, 164). Nevertheless, even when comparing our dataset to scRNA-seq performed on human renal biopsies (38), we observed a set of concordant gene expression changes. By enriching for CD45+ cells prior to sequencing, we were unable to capture CD45− splenocyte populations such as FDCs and the findings might also be biased by modulation of CD45 expression among certain cell types. Our conclusions in the present study might represent pathogenic drivers of autoimmunity or homeostatic feedback from already established autoimmunity. Furthermore, how the polymorphisms in the Sle1 and yaa loci, namely Slam/Cd2, Fcgr2b, and Tlr7, directly contribute to the phenotypes described in the present study also remains unknown. Future experiments are necessary to test these mechanistic hypotheses and potential therapeutic targets.
In this study, we present extensive transcriptional and clonotypic changes in lupus-prone mice, which may serve as a resource for ongoing profiling of SLE (165, 166). Characterization of novel myeloid subsets, CD8 metabolism, and MZ B cell loss highlights disease mechanisms for potential therapeutic targeting. Profiling of the TCR repertoire led to a model for B cell–T cell collaboration in loss of germinal center tolerance. Beyond its clear role in SLE pathogenesis, autoreactive germinal center biology and the insight gained in this study is broadly applicable to both autoimmune disease and our understanding of mechanisms of peripheral tolerance. Collectively, these findings emphasize the complex heterogeneity of the SLE immune landscape while providing potential novel therapeutic targets, particularly myeloid cell states, for further interrogation.
We thank H. Leung of the Optical Microscopy Core and J. Moore of the Flow and Imaging Cytometry Resource at the PCMM for technical assistance.
E.H.A.-G. was supported by National Institutes of Health/National Institute of General Medical Sciences Grant T32GM007753 and by National Institutes of Health/National Institute of Allergy and Infectious Diseases Grants T32AI007529 and F30AI160909. M.C.C. is supported by National Institutes of Health/National Institute of Arthritis and Musculoskeletal and Skin Diseases Grant R01 AR074105 and National Institutes of Health/National Institute of Allergy and Infectious Diseases Grant R01 AI130307.
Michael C. Carroll is a Distinguished Fellow of AAI.
Conceptualization, E.H.A.-G.; methodology, E.H.A.-G.; software, E.H.A.-G.; validation, E.H.A.-G.; formal analysis, E.H.A.-G.; investigation, E.H.A.-G.; data curation, E.H.A.-G.; writing – original draft, E.H.A.-G.; visualization, E.H.A.-G.; supervision, M.C.C.; funding acquisition, M.C.C.
The sequencing data presented in this study have been submitted to the Gene Expression Omnibus under accession number GSE192762.
The online version of this article contains supplemental material.
Abbreviations used in this article:
B cell ISG
differentially expressed gene
follicular dendritic cell
Gene Ontology biological process
log2 fold change
model-based analysis of single-cell transcriptomics
4-hydroxy-3-nitrophenyl acetyl hapten conjugated to OVA
adjusted p value
principal component 1
principal component analysis
single-cell BCR sequencing
single-cell RNA sequencing
single-cell TCR sequencing
systemic lupus erythematosus
transitional stage 1 B
transitional stage 2 B
T follicular regulatory
uniform manifold approximation and projection
y-linked autoimmune accelerating
The authors have no financial conflicts of interest.