IL-23 promotes autoimmune disease, including Th17 CD4 T cell development and autoantibody production. In this study, we show that a deficiency of the p19 component of IL-23 in the autoimmune BXD2 (BXD2-p19−/−) mouse leads to a shift of the follicular T helper cell program from follicular T helper (Tfh)–IL-17 to Tfh–IFN-γ. Although the germinal center (GC) size and the number of GC B cells remained the same, BXD2-p19−/− mice exhibited a lower class-switch recombination (CSR) in the GC B cells, leading to lower serum levels of IgG2b. Single-cell transcriptomics analysis of GC B cells revealed that whereas Ifngr1, Il21r, and Il4r genes exhibited a synchronized expression pattern with Cxcr5 and plasma cell program genes, Il17ra exhibited a synchronized expression pattern with Cxcr4 and GC program genes. Downregulation of Ighg2b in BXD2-p19−/− GC B cells was associated with decreased expression of CSR-related novel base excision repair genes that were otherwise predominantly expressed by Il17ra+ GC B cells in BXD2 mice. Together, these results suggest that although IL-23 is dispensable for GC formation, it is essential to promote a population of Tfh–IL-17 cells. IL-23 acts indirectly on Il17ra+ GC B cells to facilitate CSR-related base excision repair genes during the dark zone phase of GC B cell development.
Follicular T helper (Tfh) cells have been reported to express well-characterized transcription factors and different cytokines that define this Th cell subset as well as molecules that promote T–B cell interactions (1, 2). Upregulation of the Bcl6 transcription factor and the cytokine IL-21 are key initial events that enable the establishment of the Tfh cell program and upregulation of adhesion molecules, leading to effective, long-term, and productive interactions between T and B cells in the germinal center (GC) (3–6). The spectrum of CD4+CXCR5+PD-1+ICOS+CD40L+ Tfh cells has increased over the past several years to include multiple subpopulations (5, 7, 8). Although classical Tfh cells produce IL-21, sublineages and pathogenic Tfh cells, including IFN-γ–producing Tfh–IFN-γ, IL-4–producing Tfh–IL-4, and IL-17–producing Tfh–IL-17, have been identified in both humans and mice (9–11). In BXD2 autoimmune mice, IL-17 was shown to contribute to spontaneous GC formation and autoantibody production by enhancing and prolonging the interaction between Tfh cells and B cells through the desensitization of chemokine receptor responses through regulator of G-protein signaling (RGS) 13 and RGS16 (5, 12–15). Such a response promotes the migration arrest of CXCR4+- or CXCR5+-responding cells and thereby is critical to enable prolonged interactions of GC B cells with other cells for GC B cell affinity maturation. Similarly, Tfh–IFN-γ cells can promote the development of autoreactive B cells outside the GC through the T-bet transcription program (16). Despite these previous insights, it is not clear if these different Tfh cytokines operate concomitantly or sequentially in different niches or stages of GC B cell development to facilitate pathogenic autoantibody formation. It is also unknown if the functions of one lineage can be replaced by the others.
It is also perplexing that activation-induced cytidine deaminase (AID) expression and activity primarily occurs in the highly proliferative CXCR4+ dark zone (DZ) GC B cells, leading to extensive somatic hypermutation (SHM) (17, 18). In contrast, Tfh cells interact primarily with CXCR5+ GC B cells in the light zone (LZ) for fine-tuning of selection and plasma B cell development (19–21). We previously have shown that the spontaneous autoreactive GCs in BXD2 mice lead to the development of autoantibodies through the upregulation of AID, which promotes SHM and class-switch recombination (CSR) (22, 23).
In addition to AID, both SHM and CSR require Bcl6-directed programs of both B cell survival (24) and DNA damage and repair genes (25–27). For CSR, nonhomologous end joining (NHEJ) first requires transcription initiating upstream of the IgH constant region undergoing CSR, followed by AID-mediated dsDNA breaks and then double-stranded break repair, resulting in a synapse and joining of the two switch regions (28). Known classical NHEJ players include Ku70 (XRCC6), Ku80 (XRCC5), XRCC4, DNA ligase IV, and DNA-PKcs (29), whereas the nonclassical (or alternative) NHEJ pathway uses XRCC1 and DNA ligase III for DNA end joining (30, 31). Other DNA damage/repair gene families, including histone variants H2a and H3a, which encode core histone proteins with substantial variations in their amino and C-terminal regions, DNA polymerase δ (POLD) (32, 33), and nucleotide excision repair-like (NEIL) (34, 35), have been less well studied in CSR. Although AID is considered the master regulator of DNA breaks and NHEJ that enable CSR and we previously showed that a deficiency of IL-17RA in the lupus-prone BXD2 mice affected the expression of AID in vivo and diminished the development of IgG autoantibodies (12), how IL-23 or IL-17 regulate AID-induced NHEJ to promote CSR in immunized or spontaneously autoreactive GC B cells is largely unknown.
IL-23 is an IL-12–related cytokine composed of a p40 subunit that it shares with IL-12 and a unique p19 subunit (36). In this study, we have generated BXD2-p19−/− mice to determine if IL-23 acts through modulating Tfh cell lineage development to regulate spontaneous GC formation and class-switched autoantibody production in BXD2 mice (5, 12, 14, 15, 22, 23, 37, 38). We found that although BXD2-p19−/− mice did not have an impaired spontaneous GC response, there was a significantly lower IgG2b response that could not be compensated for by other Tfh cell subsets. By performing single-cell transcriptomics to elucidate mechanisms, we showed that IL-23 promotes AID-mediated CSR by upregulating histone variant genes, including H1fx, H2afv, H2afx, H2afz, and H3f3a as well as base excision repair genes Pold4 and Neil1. Expression of Il17ra, but not Ifngr1, Il4ra, or Il21r, was synchronized with the expression of Aicda and other GC program genes. Our results suggest a novel concept that various subsets of Tfh cells act at different stages of GC B cell development and that IL-23–generated Tfh–IL-17 cells are the important GC Tfh cell program to ensure efficient AID-mediated CSR in spontaneously autoreactive GCs in BXD2 mice.
Materials and Methods
C57BL/6 (B6) and BXD2 recombinant inbred mice were purchased from The Jackson Laboratory. B6–IL-23 p19−/− mice were kindly provided by Merck Research Laboratories (Palo Alto, CA) (39, 40). B6-Il17ra−/− mice were obtained from Amgen (Thousand Oaks, CA). B6–IL-23 p19−/− and B6-Il17ra−/− mice were backcrossed with BXD2 mice using a marker-associated speed congenic approach for seven generations to generate BXD2-p19−/− and B6-Il17ra−/− mice (5).
Administration of adenovirus LacZ and adenovirus IL-23
Adenovirus (Ad) IL-23 and AdLacZ (2 × 109 PFU per mouse), generous gifts from Dr. J. Kolls (University of Pittsburgh) (41), were administered i.v. Spleen tissue and cells were analyzed 5 d later.
Flow cytometry analysis and sorting
Single-cell suspensions of spleen were harvested and subjected to standard cell-surface staining or intracellular staining following the manufacturer’s instructions. Anti-mouse Abs included BioLegend Pacific Blue anti-CD19 (clone 6D5), PE anti–IL-23R (clone 12B2B64), allophycocyanin anti–IFN-γ (clone XMG1.2), PE-Cy7 anti-IgM (clone RMM-1), BD Biosciences PE anti-CD95 (clone Jo2), PE-Cy7 anti-CXCR5 (clone 2G8), Brilliant Violet 510 anti-IgD (clone 11-26c.2a), Thermo Fisher Scientific eFluor 660 anti-GL7 (clone GL7), FITC anti–PD-1 (clone J43), and Alexa Fluor 647 anti–IL-17A (clone eBio17B7). Peanut agglutinin (PNA) was conjugated with biotin (Vector Laboratories, Burlingame, CA) and detected by Alexa Fluor 488–conjugated streptavidin (Invitrogen). Dead cells were excluded from analysis with allophycocyanin–eFluor 780 Organic Viability Dye (eBioscience).
For cytokine-producing T cell analysis, cells were stimulated for 5 h with PMA (50 ng/ml; Sigma-Aldrich) and ionomycin (750 ng/ml; Sigma-Aldrich) in the presence of GolgiPlug (BD Biosciences). Cells were stained for surface markers and then fixed and permeabilized with Cytofix/Cytoperm solution (BD Biosciences) before intracellular staining (5).
Data were acquired using standard flow cytometry (LSR II; BD Biosciences) and analyzed with FlowJo_v10 software. FACS sorting was performed using a FACSAria cell sorter (BD Biosciences). Both instruments were located in the University of Alabama at Birmingham (UAB) Comprehensive Flow Cytometry Core in the Shelby Biomedical Research Building.
Histology and confocal imaging analysis
Frozen sections of mouse spleens were subjected to histology and confocal imaging. Immunofluorescent staining was performed as previously described (14). Images were taken using a Nikon A1 at the UAB High Resolution Imaging Facility. Imaging data analysis was carried out using the ImageJ software, version 1.47, developed by the U.S. National Institutes of Health (42).
Real-time quantitative RT-PCR
The expression of genes in sorted GC B cells and Tfh cells was determined using quantitative RT-PCR (qRT-PCR) following the previously described method (12, 43). Spleen cells were purified using flow cytometry, gated to either live CD4+PD-1+CXCR5+ Tfh cells or live CD19+GL7+Fas+ GC B cells. Primers used are listed in Supplemental Table I.
Single-cell library generation and sequencing
Spleen cells were FACS sorted into live CD19+GL7+Fas+ GC B cells. Single cells were captured into Gel Bead-In EMulsions using a 10× chromium controller, from which single-cell 5′-biased transcriptome libraries and V(D)J libraries were prepared using a 5′ Library and Gel Bead Kit (PN-1000014) and Single Cell V(D)J Enrichment Kit for Mouse B Cells (PN-1000016, 10× Genomics), according to the 10× Genomics manual. Final 5′-biased transcriptomic libraries were sequenced using an Illumina NextSeq 500 with a targeting minimum of 50,000 read pairs/cell and 26 × 98–bp read lengths for the sequencing cycles. Final enriched V(D)J libraries were sequenced using an Illumina MiSeq with a minimum of 5000 read pairs/cell and 150 × 150–bp read lengths.
Processing of single-cell sequencing raw data
Raw sequencing data were processed with the Cell Ranger pipeline software (v.3.0.2; 10× Genomics). Raw base call files generated by NextSeq 500 or MiSeq (Illumina) were converted to FASTQ files using cellranger mkfastq with default parameters. For single-cell 5′ transcriptome data, cellranger count was run with –transcriptome = refdata-cellranger-mm10-3.0.0 for each sample, in which reads were aligned against the mouse genome (mm10) using the spliced transcripts alignment to a reference aligner (44), and the unique molecular identifiers were counted for each gene. The gene expression matrix in the outputs of cellranger count was used for downstream analysis. Single-cell V(D)J BCR data were processed by the cellranger vdj pipeline to perform quality control, assembling, quantification, and annotation of paired V(D)J transcript sequences. Briefly, after trimming and filtering, a de novo assembly algorithm from cellranger vdj pipeline was applied to the reads of each cell barcode. A set of assembled contiguous sequences of overlapping DNA fragments (contigs), which represent the best estimation of transcript sequences, were generated via the Smith–Waterman algorithm (45) in cellranger vdj pipeline. The assembled contigs in each cell were aligned against all of the sequences from mouse reference genome assembly mm10 in this pipeline. The filtered contig annotation files were used for downstream analysis.
Classification of B cells using single-cell transcriptome data and single-cell BCR data
Seurat (46) (v.3.0.0), implemented using the R package, was applied to exclude low-quality cells in different single-cell experiments. Cells that expressed fewer than 200 genes were filtered out. Genes that were not detected in at least three single cells were excluded. Based on these criteria, we retained the following numbers of genes/cells in each sample: 14,255 genes/2437 cells in BXD2-p19−/− no. 1, 13,366 genes/1084 cells in BXD2-p19−/− no. 2, 14,233 genes/2586 cells in BXD2-p19−/− no. 3, 14,479 genes/3158 cells in BXD2-wild type (WT) no. 1, and 14,195 genes/ 1936 cells in BXD2-WT no. 2. The filtered contig annotation files generated from single-cell V(D)J BCR data were merged with filtered single-cell transcriptome databased on barcode information in R. Single cells with overlapped barcodes between BCR data and transcriptome data were identified as B cells. The transcriptome datasets for B cells were loaded into the R package Seurat. The Seurat data matrix object for each dataset was created using the “CreateSeuratObject” function. The Seurat objects for B cells were filtered using the function “subset.” The individual B cells were selected for positive expression of the coding gene for Ms4a1 (CD20) and CD19, key B cell markers. Only cells that expressed the VDJ regions of BCR as well as Ms4a1 (the coding gene for Cd20) and Cd19 were used for further analysis. The numbers of single cells used for the transcriptomics analysis for each mouse are 566, 299, 984, 1136, and 567 for BXD2-p19−/− no. 1, BXD2-p19−/− no. 2, BXD2-p19−/− no. 3, BXD2-WT no. 1, and BXD2-WT no. 2, respectively.
Identification of differentially expressed genes
The procedures for normalization, variable gene selection, scaling the data, and identifying differentially expressed genes were performed by following the Seurat documentation (v.3.0.0) (47, 48). In brief, the data were normalized using the function “NormalizeData,” which uses a global-scaling normalization method “LogNormalize” to normalize the gene expression measurements for each cell to the total gene expression. The scale factor was 10,000. The highly variable genes were identified by the “FindVariableFeatures” function, and the method of selection was “vst.” Next, the data were scaled using the function “ScaleData” in Seurat. The guided clustering was performed based on the information of samples and cells. The B cell cluster was identified as described above. The samples BXD2-p19−/− no. 1, BXD2-p19−/− no. 2, and BXD2-p19−/− no. 3 were identified as “ko” cluster, and the samples BXD2-WT no. 1 and BXD2-WT no. 2 were identified as “wt” cluster. Differential gene expression analyses were performed using the Seurat function “FindMarkers.” To identify differentially expressed genes, the Wilcoxon rank-sum test was performed with the default threshold of 0.25 for log2 fold change and a filter for the minimum percentage of cells in a cluster of >25%. Differentially expressed genes were isolated by comparing significantly upregulated genes and downregulated genes defined as having an adjusted p value of <0.05. The top differentially expressed genes were visualized by heatmap using the ComplexHeatmap R package.
Single-cell trajectory analysis
R package Monocle (2.6.4) (47, 48) was used to investigate the developmental trajectories of single cells by a method similar to the procedure described in Qiu et al. (47, 48). The Seurat objects from the last step were loaded into Monocle. To reduce the bias caused by technical variation, sequencing depth, and capture efficiency, data normalization was performed using the functions “estimateSizeFactors” and “estimateDispersions.” The low-quality cells (i.e., doublets, triplets or libraries that were accidentally formed from multiple cells) were filtered out. The minimum expression threshold was set at 0.1. Genes were filtered out based on the average expression level and unusual variations across cells. For single-cell trajectories in B cells, a set of ordered genes that define B cell development were collected to order cells for supervised trajectories. These genes included Ighm, Ighd, Ighg2b, Ighg2c, Ighg1, and other key genes involved in GC development pathways. The expression profiles were reduced to two dimensions using the DDRTree algorithm in the function “reduceDimension.” Cells were ordered into a trajectory using the “orderCells” function. The function “plot_pseudotime_heatmap” was used to visualize modules of genes that covaried across pseudotime. The “plot_genes_in_pseudotime” function was used for individual graphs to depict the kinetic trend of the expression levels of genes with significant changes through pseudotime.
Kyoto Encyclopedia of Genes and Genomes pathway analysis
The clusterProfiler (49) package (v3.10.1) was applied to identify enriched pathways based on top differentially expressed genes (adjusted p value <0.05). The ranked gene list was prepared following the the clusterProfilter documentation. The number of selected genes associated with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was assessed via the function “enrichKEGG.” The p value cutoff was 0.05. The enrichment result was visualized by bar and dot chart using “barplot” and “dotplot” functions, respectively. The “cnetplot” function was used to visualize the enrichment map and extract the potential biological complexities associated between genes and diseases.
SeqGeq single-cell RNA sequencing data analysis
Transcriptome gating analysis was carried out using the SeqGeq platform (v1.6.0; BD Life Sciences – Informatics, Ashland, OR). Following the use of gene expressed versus library size quality control, B cells were gated as BCR gene-positive cells. SeqGeq analysis was used to gate on cell subsets based on the expression of gene-encoding Tfh cytokine receptors and then enabled analysis in quantitation of expression of significantly altered gene pairs including the histone variant genes, DNA excision repair genes, and Rgs13 and Aicda. The axis labels show that this indicates relative expression of the indicated genes on a hyperbolic arcsine (Arcsinh) transformation scale (50, 51). The arcsinh values are calculated by applying the arcsinh equation divided by the scale argument to the measured intensity value data that are displayed in a linear-like fashion. The statistical analysis was calculated as a Chi-Square significance of each percentage relative to the percentage for the indicated gene pair for GC B cells expressing the Il17ra+ gene in WT BXD2 compared with the same gene pair in other cytokine receptor–positive B cells in both WT and BXD2-p19−/− mice.
Results are shown as the mean ± SD or mean ± SEM. Pearson’s normality test was used to determine the normal distribution of each dataset. A two-tailed, unpaired Student t test was used when two normally distributed groups of datasets were compared for statistical differences. The Mann–Whitney nonparametric test was used when the two datasets for comparison were not normally distributed. An ANOVA test was used when more than two groups were compared for statistical differences. Multiple comparisons were performed using Tukey multiple comparisons test. Fisher exact test was used for contingency tests. Analysis of correlation between variables was performed using linear regression in GraphPad Prism software. The p values of <0.05 were considered significant.
The single-cell transcriptome data and single-cell V(D)J BCR sequencing data used in this study have been deposited in the Gene Expression Omnibus database under the accession number GSE145922 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE145922). Raw files used in the figures that support the findings of this study are available from the corresponding authors upon reasonable request.
IL-23 regulates the Tfh cell switch from Tfh–IFN-γ to Tfh–IL-17 but does not change Tfh cell frequency
We have generated p19-deficient BXD2 mice (BXD2-p19−/−) to determine if the lack of IL-23 can alter Tfh cell lineage development and if such an alternation can influence spontaneous GC and autoantibody production. The percentage and the absolute number of classical Tfh cells (CD4+PD-1+CXCR5+) were similar in BXD2 and BXD2-p19−/− mice (Fig. 1A, 1B). Importantly, however, the IL-23 deficiency significantly suppressed the expression of IL-17 while enhancing the expression of IFN-γ among Tfh cells in the BXD2-p19−/− mice (Fig. 1C, 1D). In addition, we identified that Tfh cell phenotype between BXD2-p19−/− and BXD2-Il17ra−/− mice are not identical in that there was no decrease in the percentages of IL-17 or increase in the percentage of IFN-γ in IL-17RA–deficient BXD2 mouse Tfh cells compared with WT BXD2 mice (Fig. 1C, 1D). These results suggest that although IL-23 promotes Tfh–IL-17 development, the effects of IL-23 and IL-17 on autoreactive GC development may be different.
We next determined if excessive production of IL-23 in vivo can promote Tfh–IL-17 development. This was tested using an IL-23–expressing Ad system. We previously showed that after AdIL-23 administration to B6 mice, the level of serum IL-23 increased over baseline by day 3, achieved its maximal level by day 5, and then returned to baseline by day 10 (52). Consistent with the pro–IFN-γ effects of BXD2-p19−/− on Tfh cells described above, the overexpression of IL-23 in B6 mice promoted the expression of Tfh–IL-17 cytokine genes (Il17a and Il17f) (Fig. 1E) but hindered the expression of the Tfh–IFN-γ transcriptional factor (Tbet) and cytokine (Ifng) genes in Tfh cells after 5 d of AdIL-23 administration (Fig. 1E). AdIL-23 also suppressed the expression of canonical Tfh cells (Bcl6 and Il21) as well as Tfh–IL-4 transcription factor and cytokine (Gata3 and Il4) genes (Fig. 1E). These results suggest a pro-Tfh–IL-17 effect of IL-23 and that such effects can lead to compensatory alteration in other lineages of Tfh cells in the absence of IL-23.
IL-23 is important for B cell CSR
We have shown that IL-23 positively drives its own signal through the induction of the IL-23R (52). The levels of Il23r in GC B cells in both AdLacZ- and AdIL-23–administered B6 mice were very low. In contrast, there was an almost 100-fold higher level of Il23r in Tfh cells compared with GC B cells (Fig. 2A). Further, AdIL-23 administration significantly induced the expression of Il23r in Tfh cells of both B6 and BXD2-p19−/− mice (Fig. 2A, 2B), suggesting that Tfh cells, but not GC B cells, are the major target of IL-23.
Having established a positive effect of IL-23 in promoting Tfh–IL-17 and diminishing other lineages of Tfh cells, we next determined if the lack of IL-23 would influence the spontaneous GC response and autoantibody formation in BXD2 mice. Comparison of the serum autoantibody levels in B6, BXD2, BXD2-p19−/−, and BXD2-Il17ra−/− mice revealed that there were lower levels of IgG anti-DNA and rheumatoid factor (RF) in BXD2-p19−/− mice compared with BXD2 mice, although the anti-DNA and RF IgM autoantibody levels were significantly increased (Fig. 2C). Interestingly, these results cannot be fully explained by the lower signaling through IL-17RA because the IgG and IgM anti-DNA and RF IgM levels in serum of BXD2-Il17ra−/− mice were decreased, as we previously reported (Fig. 2C) (5, 12). Consistent with the serum autoantibody levels, there was an increase in IgM-producing plasmablasts (PBs) and a decrease in IgG-producing PBs in the BXD2-p19−/− mice, whereas both IgM- and IgG-producing PBs were diminished in BXD2-Il17ra−/− mice compared with BXD2 mice (Fig. 2D–F). The results thus far suggested that IL-23 promoted CSR in BXD2 mice.
IL-23 is not required for GC formation
Unlike BXD2-Il17ra−/− mice, which had smaller GCs and a lower frequency of Fas+GL7+ GC B cells, as we previously reported (5, 12), the size of spontaneous GCs and percentages of Fas+GL7+ GC B cells were slightly but not significantly increased in BXD2 mice lacking IL-23–p19 (Figs. 2D, 3A, 3B). The results indicated that IL-23 is dispensable for GC formation. Consistent with this, the administration of AdIL-23 to BXD2-p19−/− mice did not change the frequency and numbers of GC B cells (Fig. 3C, 3D). There was detection of IgD+PNA+GL-7+ GC B cells in AdIL-23–administered BXD2-p19−/− mice, suggesting that founder GC B cells were recruited to the GC (53). However, there was an increased percentage of PNA+GL-7+ IgM−IgD− switched GC B cells in BXD2-p19−/− mice administered AdIL-23 compared with the AdLacZ group of mice at day 5 (Fig. 3E). Taken together, these results suggested that IL-23 is dispensable for GC formation, but AdIL-23 can induce B cell CSR.
Decreased expression of IgG2b in the spleens of BXD2-p19−/− mice
To further delineate the effects of IL-23 and therefore Tfh–IL-17 in regulating GC B cell CSR, single-cell RNA sequencing (scRNA-seq) was carried out using sorted GC B cells obtained from two WT BXD2 mice (1703 cells) and three BXD2-p19−/− mice (1849 cells). Consistent with the circulating autoantibody results, a Seurat-based differentially expressed gene analysis revealed that Ighm and Cd38 were two highly overexpressed genes, whereas Ighg2b was a highly underexpressed gene in BXD2-p19−/− mice compared with WT BXD2 mice (Fig. 4A). Elevated IgM and reduced levels of IgG2b in the serum BXD2-p19−/− mice compared with WT BXD2 mice were verified using a conventional ELISA analysis (Fig. 4B). There was also a significantly lower intensity of IgG2b plasma B cells in the periarteriolar lymphatic sheaths (PALS) in the spleen of BXD2-p19−/− mice, confirming that the production of IgG2b in BXD2-p19−/− mice was significantly compromised (Fig. 4C, 4D).
IL-23 promotes the induction of DNA excision repair genes in BXD2 Ighg2b+ GC B cells
Because lower Ighg2b and IgG2b is a signature of BXD2-p19−/− mice, we took advantage of scRNA-seq transcriptomics analysis and determined the top 100 downregulated genes in Ighg2b+ GC B cells of BXD2-p19−/− versus WT mice (Fig. 5A). Pathway analysis confirmed that cell cycle genes (Ccna2 and Ccnb2), systemic lupus erythematosus autoantibody-targeted histone variant genes (H1fx, H2afv, H2afx, H2afz, H3f3a, and H3f3b), and DNA base excision repair genes (Pold4 and Neil1) are among the top downregulated pathways in Ighg2b+ GC B cells derived from BXD2-p19−/− mice (Fig. 5B, 5C). Violin plot analysis further demonstrated that histone variant genes (H1fx, H2afv, H2afx, H2afz, and H3f3a), DNA excision repair genes (Pold4 and Neil1), calmodulin-related and cell cycle genes (Calm1 and Calm2), and Rgs13 were all significantly downregulated in Ighg2b+ GC B cells in BXD2-p19−/− mice versus WT BXD2 mice (Fig. 5D).
Il17ra synchronized with the expression of GC program genes in spontaneous autoimmune GCs of BXD2 mice
To further analyze how different Tfh cytokine responses can be associated with the induction of Ighg2b in GC B cells, we have applied a pseudotime trajectory analysis of single-cell transcriptomes of GC B cells derived from WT and BXD2-p19−/− mice (Fig. 6A). The GC B cell developmental stages (pseudotime 0–3.5) were ordered as early (0–1), mid (1–3), and late phase (>3) by pseudotime analysis based on the developmental stages, which were based on the downregulation of Ighm and Ighd and the upregulation of Ighg BCR genes (Fig. 6A, 6B). This organization enabled us to determine GC B cell gene expression from the earliest (Ighm+Ighd+) to the latest plasma cell phase (Ighg+Prdm1+) of development. The GC B cell location was ordered based on the expression of the DZ chemokine receptor gene Cxcr4 and the LZ chemokine receptor gene Cxcr5 (Fig. 6A, 6C).
The pseudotime expression of Tfh cytokine receptor genes, GC B cell program signature genes, and plasma B cell program signature genes was then analyzed. Il4ra expression was synchronized with the expression of Ighm and Ighd in the early phase of GC B cell development in BXD2 mice (Fig. 6A, 6C), a result consistent with transcriptomic analysis of human naive IgD+ B cells as reported by others (54–56). In WT and BXD2-p19−/− mice, Il17ra gene expression was synchronized with the downregulation of Ighm, Ighd, and Il4ra. This was also a stage during which there was initially a slow but progressive upregulation of Ighg1, Ighg2b, and Ighg2c (Fig. 6A, 6B) and with the expression of the DZ chemokine receptor gene Cxcr4 (Fig. 6A, 6C). Consistent with our previous findings (12–15), Il17ra was also synchronized with GC program signature genes including Aicda, Bcl6, and Rgs13 (Fig. 6A, 6C, 6D). Ifngr1 and Il21r exhibited steady but low levels from the Ighm+ and Ighd+ early preswitch phase to the midphase of GC B cell development. A late-phase increase of Ifngr1 and Il21r was correlated with a further elevation of the Ighg genes, the plasma cell program genes Prdm1 and Id2, and the LZ chemokine receptor gene Cxcr5 (Fig. 6A, 6C, 6E).
Although cytokine receptor genes (Il17ra, Ifngr1, and Il21r), chemokine receptor genes, GC signature genes, and plasma cell signature genes showed a similar pseudotime expression kinetics in WT BXD2 and BXD2-p19−/− GC B cells (Fig. 6A, 6C), the lower expression of Ighg2b in BXD2-p19−/− GC B cells became apparent at the mid–GC phase and continued through the late plasma phase of development (Fig. 6A, arrow, Fig. 6B). Also, in BXD2-p19−/− GC B cells, Il4ra expression kinetics were shifted to the late phase of GC B cell development and were synchronized with the expression of plasma cell program genes (Fig. 6A, 6E).
Il17ra synchronized with the expression of histone variant genes and DNA excision repair genes at the GC program phase
The above results suggest that Il17ra is the only Tfh cytokine receptor gene that synchronized with the GC B cells during the Aicda+ DZ phase of GC B cell development. We therefore determined if the major downregulated pathway genes (histone variants, cell cycle genes, and excision repair genes) in Ighg2b+ GC B cells from BXD2-p19−/− mice also exhibited the same expression kinetics as Il17ra genes. Indeed, the Monocle-based pseudotime analysis showed that all major histone variant genes including H1fx, H2afv, H2afx, H2afz, and H3f3a, calmodulin-related and cell cycle genes including Calm1, Calm2, Cdca3, and Ccna2 genes, and DNA excision repair genes Pold4 and Neil1 were expressed at the Aicda+/Cxcr4+ GC DZ phase of development in GC B cells derived from both WT and BXD2-p19−/− mice (Fig. 7A).
To determine if Il17ra+ GC B cells indeed had higher transcript levels of these GC B cell genes that may be involved in the DZ cell cycle and AID-mediated excision repair, a SeqGeq gating analysis was used to identify GC B cells that expressed only Il17ra, Ifngr1, Il21r, or Il4r (Fig. 7B, 7C, upper panels). There was a higher percentage of Rgs13+Aicda+, H1fx+H2afx+, Calm1+Calm2+, and Pold4+Neil1+ cells among the Il17ra+ GC B cells compared with Ifngr1+, Il21r+, and Il4r+ GC B cells in both WT BXD2 and BXD2-p19−/− mice (Fig. 7B, 7C, lower panels). Among these, the distribution of Rgs13+Aicda+, Calm1+Calm2+, and Pold4+Neil1+ cells in Il17ra+ GC B cells of BXD2-p19−/− mice were significantly lower than that in Il17ra+ GC B cells of WT BXD2 mice (Fig. 7B–D). These results also showed that other cytokines produced by Tfh cells cannot compensate for the loss of IL-17RA signals in BXD2-p19−/− mice to enable effective expression of these gene programs that are needed for DNA repair associated with AID-mediated CSR.
Ab-forming B cell development in the GC is a complex process involving B cell commitment, migration, migration arrest, cell cycle entry, AID induction, AID-induced DNA strand breaks, and excision repair, leading to differentiation of GC B cells into Ab-secreting plasma cells (57). The classical CXCR5+ LZ–oriented IL-21–producing Tfh cell model cannot fully explain the orchestration of Tfh cells in directing all of these processes. This is especially applicable during the CXCR4+ DZ stage of GC B cell development in which GC B cells undergo extensive proliferation, AID induction, and a series of AID-mediated effector functions needed for B cell affinity maturation (57, 58). We previously showed that the IL-17RA signal is important in mediating migration arrest to both CXCL12 and CXCL13 (5, 12–15). We also showed that IL-17 and IL-21 were expressed by distinct Tfh cell subpopulations, with IL-17 predominating in the CXCR5−ICOS+ subpopulation, whereas IL-21 expression was confined to the CXCR5+ subpopulation (5). The present study further suggests that IL-23–mediated induction of Tfh–IL-17 cells plays an essential role in shaping the autoantibody formation response through its effects on promoting CSR leading to IgG2b expression during the DZ Aicda+ GC program stage of B cell maturation lupus-prone BXD2 mice. The lack of such signaling in BXD2-p19−/− mice is associated with diminished autoantibody CSR to the IgG2b but not the IgG1 isotype.
Surprisingly, the diminished autoantibody production in BXD2-p19−/− mice was only partially identical to the effects of IL-17RA deficiency on the spontaneous GC responses. We previously showed that IL-17 was important to upregulate RGS13 to enable GC B cell migration arrest and close interaction with GC T cells (5, 12–15). Thus, IL-17 acts as a “glue” to hold the GCs together as well as to promote efficient T–B cell interactions, leading to increased AID expression and the resultant increased CSR and SHM (5, 12). In the absence of IL-17RA, efficient GC T–B cell interactions do not occur during GC development, and the development of GC cells and Ig CSR is abrogated. In the current study, we found that although the expression of Rgs13 in GC B cells was reduced and there was a significant reduction in the percentages and numbers of IL-17–expressing Tfh cells in BXD2-p19−/− mice, the percentages and numbers of GC B cells were comparable between WT and BXD2-p19−/− mice. Also, the administration of AdIL-23 to B6 mice did not enhance the size of GC B cells in both B6 and BXD2-p19−/− mice. These results, therefore, suggest that there are additional factors that can maintain the size and numbers of GCs in the absence of IL-23. These factors, however, are not affected by the absence of IL-17RA.
Analysis of the lineages of Tfh cells provided insights into the differential effects of IL-23 p19 deficiency versus IL-17RA deficiency on GC and autoantibody development. We have found that the absence of IL-23 diminished the percentages of Tfh–IL-17 cells and dramatically skewed Tfh cells into Tfh–IFN-γ cells, whereas the absence of IL-17RA did not promote Tfh–IFN-γ cell polarization.
One possible explanation for the enhanced Tfh–IFN-γ in BXD2-p19−/− mice is that in the absence of the p19 component, the dynamics of IL-12 p40/p35 heterodimer or p40/p40 homodimer may be altered. Constitutive higher levels of p40 expression has been shown to be critical for the IL-12 p40/p35 production (59). IL-12 p40/p40 homodimer, however, has been shown as an antagonist of IL-12 p40/p35 function and blocks IFN-γ production (60). As IFN-γ was enhanced in Tfh cells in BXD2-p19−/− mice, it is unlikely that there was an excessive IL-12 p40/p40 production. Although we do not rule out the possibility that deletion of p19 potentially enhances IL-12 p40/p35 production, resulting in a shift toward the IFN-γ pathway, the administration of excessive IL-23 using the Ad system leads to the induction of Th17 cytokine genes, but it also downregulated the expression of Th1, Th2, and Tfh cytokine and transcription factor genes. These results therefore suggest that IL-23 most likely influences the Tfh cell polarization program to regulate the Tfh cell lineage development.
The excessive skewing of the Tfh cells into other non–Tfh–IL-17 lineages of cells observed in the BXD2-p19−/− mice has provided the needed signals to initiate and maintain the size of the spontaneous GCs. They also suggest that there are unique IL-17–mediated B cell maturation responses in GC B cells that cannot be completely compensated for by other Tfh cytokines. Perturbation of IL-23 in the BXD2-p19−/− mouse provided an opportunity to interrogate the pro–IL-17– versus pro–IFN-γ–mediated GC B cell developmental responses and signatures in an autoimmune mouse. The LZ is more widely studied for centrocytes and T cell interactions because these GC B cells closely interact with the canonical CXCR5+ Tfh cells for commitment to the plasma cell fate (21, 61). However, our scRNA-seq differentially expressed gene and pseudotime results show that upregulation of Aicda is coordinated with the upregulation of Bcl6, Cxcr4, and Pax5, which have been shown to occur in the DZ (62). Interestingly, pseudotime analysis suggests that the Il17r expression in GC B cells is synchronized with the expression of Aicda and other transcriptome signatures of DZ centroblasts, including Pax5, Bcl6, Rgs13, Tcf3 (which encodes E2A), and Cxcr4. In contrast, Ifngr1 and Il21r are best synchronized with the transcriptome signature of LZ centrocytes/PBs, including Prdm1, Irf4, Xbp1, Id2, and Cxcr5. Pax5 directly promotes AID induction, and E2A has been shown to induce B cell CSR (63, 64). In contrast, Blimp1 and Id2 repress AID and Pax5 and thereby suppress CSR (65, 66). The synchronized expression pattern of Il17ra with Aicda and other GC program genes suggest that Tfh–IL-17 cells may operate in the DZ to facilitate AID and its downstream effector functions. This is consistent with our previous results showing that IL-17 was expressed predominantly in CXCR5−ICOS+ CD4 T cells, whereas IL-21 was expressed in CXCR5+ICOS+ CD4 T cells in BXD2 mice (5).
Coexpression analysis further shows that Il17r+ GC B cells expressed higher levels of genes related to genome instability and DNA excision repair compared with Ifngr1+, Il4ra+, and Il21r+ GC B cells. Although we did not find differences in pseudotime gene expression kinetics between WT and p19−/− GC B cells, the pseudotime gene trajectory analysis has enabled identification of novel DNA excision repair genes that might be related to DZ B cell excision repair, leading to IgG2b class switch. In the classical NHEJ excision repair, most studies thus far focused on the role of H2afx and the ligase complex of DNA ligase IV, XRCC4, XRCC5 (Ku80), and XRCC6 (Ku70) (28, 67). Although the current study does not negate the previous findings, the scRNA-seq gene expression and trajectory analysis reveal that there are numerous novel factors that work separately during the DZ or the LZ phase to enable the CSR response. In the histone variate family members, we have identified that the expression of H1fx, H2afv, H2afx, H2afz, and H3f3a were synchronized with the expression of Aicda, and these genes were downregulated in GC B cells of BXD2-p19−/− mice. Additional genes that were identified included Pold4, Neil1, Calm1, Calm2, and Xrcc1. Although the role of POLD4 in immunity has not been identified previously, a recent study identified that two patients with biallelic mutations affecting the POLD1 and POLD2 subunits of polymerase δ exhibited developmental defects and combined immunodeficiency (32). Further, NEIL1, a human H2TH family glycosylase, and Calm1 have been implicated in the repair of lesions in transcription bubbles (68) and histone H2AX phosphorylation (69), respectively. The present scRNA-seq data therefore suggest a novel group of complex factors are at play in AID-mediated CSR. Importantly, the current study also suggests that the T cell factor that facilitates these processes is IL-17. We propose that such IL-17–producing cells be considered as T–GC helper cells, an underappreciated GC T cell subpopulation that facilitates the selection and maturation of centroblasts and can be distinguished from the well-established Tfh cells that operate in the LZ to regulate a GC response.
Interestingly, the effect of IL-23 deficiency is unique to the development of IgG2b-producing B cells. Because IgG1 transcripts and circulating Ab levels were not affected in BXD2-p19−/− mice, these results suggest that IgG1 production either does not require IL-23 and Tfh–IL-17 or was rescued by the elevation of other Tfh cytokines. This is consistent with an earlier model in which there is an hierarchy of cytokines for inducing isotype switching such that IFN-γ is dominant over IL-4 for the survival of IgG2a versus IgG1 and IgE switched GC B cells (70). In contrast, CSR to IgG2b has been shown to be TGF-β dependent but IFN-γ and IL-4 independent (70, 71). This is also consistent with previous observations that a high dose of IL-17 enhanced IgG2b but not IgG1 CSR (72) and suggests that IL-17R signaling may be uniquely synchronized with AID and genes essential for NHEJ, leading to IgG2b formation.
Our results suggest a unique effect of IL-23 in orchestrating IgG2b CSR in the BXD2 mouse model of lupus. The results therefore underscore the necessity for a coordinated upregulation of not only AID but also the DNA manipulation machinery that occurs in an environment that contains IL-23. The findings of the current study are also relevant to the prominent Phase II clinical trial of ustekinumab (Stelara) for systemic lupus erythematosus patients with active disease (73). Ustekinumab is a human mAb that binds to the p40 subunit shared by both IL-12 and IL-23 and will not lead to a skewed B cell response to Tfh–IFN-γ cells as demonstrated in the p19-deficient condition described in this study. Based on the results generated in this work, we anticipate that ustekinumab will diminish the IL-23–mediated Tfh–IL-17 cells and their AID-mediated effector functions in the DZ and will also diminish the effects of IFN-γ in assisting canonical Tfh cells in the LZ. However, ustekinumab treatment also has been shown to induce subacute cutaneous lupus (74). It will be important to determine if such effects are relevant to the undesirable skewing toward Th2 responses by IL-12/IL-23 blockade.
We thank Drs. Michael Crowley and David K. Crossman at the UAB Heflin Genomics Core Laboratory for carrying out the RNA sequencing assay and assembling the Cell Ranger pipeline.
This work was supported by the National Institutes of Health (R01-AI-071110, R01 AI134023, and P30-AR-048311 to J.D.M. and R01-AI-083705 to H.-C.H.), U.S. Department of Veterans Affairs grants (I01BX004049 and 1I01BX000600 to J.D.M.), and the Lupus Research Alliance (a distinguished innovator award to J.D.M. and an alliance target identification award to H.-C.H.).
The single-cell RNA sequencing data presented in this article have been submitted to the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE145922) under accession number GSE145922.
The online version of this article contains supplemental material.
Abbreviations used in this article:
activation-induced cytidine deaminase
contiguous sequence of overlapping DNA fragments
Kyoto Encyclopedia of Genes and Genomes
nucleotide excision repair-like
nonhomologous end joining
periarteriolar lymphatic sheath
DNA polymerase δ
regulator of G-protein signaling
single-cell RNA sequencing
follicular T helper
University of Alabama at Birmingham
The authors have no financial conflicts of interest.