Single-cell transcriptome landscape of circulating CD4+ T cell populations in autoimmune diseases

Abstract

CD4⁺ T cells are key mediators of various autoimmune diseases; however, their role in disease progression remains unclear due to cellular heterogeneity. Here, we evaluated CD4⁺ T cell subpopulations using decomposition-based transcriptome characterization and canonical clustering strategies. This approach identified 12 independent gene programs governing whole CD4⁺ T cell heterogeneity, which can explain the ambiguity of canonical clustering. In addition, we performed a meta-analysis using public single-cell datasets of over 1.8 million peripheral CD4⁺ T cells from 953 individuals by projecting cells onto the reference and cataloging cell frequency and qualitative alterations of the populations in 20 diseases. The analyses revealed that the 12 transcriptional programs were useful in characterizing each autoimmune disease and predicting its clinical status. Moreover, genetic variants associated with autoimmune diseases showed disease-specific enrichment within the 12 gene programs. The results collectively provide a landscape of single-cell transcriptomes of CD4⁺ T cell subpopulations involved in autoimmune disease.

Keywords: CD4(+) T cells; GWAS; autoimmune diseases; immunogenomics; single-cell RNA-seq.

Graphical abstract

Figure 1

Global profiling of CD4⁺ T cells (A) Sample collection strategy. HC, healthy control; MG, myasthenia gravis; MS, multiple sclerosis; SLE, systemic lupus erythematosus. (B and C) Cluster layer 1 (L1) and layer 2 (L2) on Uniform Manifold Approximation and Projection (UMAP) embeddings. (D) Dot plot depicting signature genes' mean expression levels and percentage of cells expressing them across clusters. Marker genes for the plot were manually selected. See also Figure S1C for automatically extracted marker genes. (E) Expression correlation of clusters with bulk RNA-seq for sorted CD4⁺ T cell fractions generated by the DICE project. (F) UMAP plot showing clonotype size. The color indicates the number of cells sharing the same T cell receptor (TCR) sequences. (G and H) Clonotype size distributions across clusters. (I and J) TCR similarity networks in autoimmune patients, healthy donors (I), and all donors (J). TCR similarity was calculated for each sample, and only edges where overlapping clonotypes were detected in ≥2 (I) or ≥4 (J) samples are depicted as robust overlaps. The edge color indicates the average TCR similarity of all samples. (K) Degree centrality of TCR networks. Significance across clusters was calculated by one-way analysis of variance and, after multiple test corrections by false discovery rate, Tcm (Tfh) and Treg Act were retained as significant cell types. Then, pairwise Tukey-HSD post hoc tests were performed. ∗p_adj < 0.05 in comparison with HCs. (L) TCR-intrinsic regulatory potential (TiRP) score distributions on UMAP plot. Mean scores for each cluster L2 were shown. (M) TiRP score distribution across cluster L2. The dot shows the mean, and the confidence interval (CI) shows 95% CI of the bootstrap distribution of means (n = 1,000). Adjusted p values of significant clusters, Tnaive MX1 1.29 × 10⁻², Tem (Th1) 5.82 × 10⁻¹¹, Temra (Th1) 4.84 × 10⁻²², Treg Naive 2.14 × 10⁻¹³, Treg Act 2.14 × 10⁻¹³, Treg Eff 3.82 × 10⁻⁵ (two-sided Mann-Whitney U test was performed for one cluster vs. the other clusters iteratively).

Figure 2

NMF captured 12 CD4⁺ T cell features (A) Schematic view of NMF and NMF projection (NMFproj). (B) Matrixplot showing the mean scaled NMF feature weight for each cluster L2 population. The explained variance (Evar) is also shown on the right. The NMF feature weight is scaled by the maximum value for each feature for visualization. (C) NMF cell feature value on UMAP plots. (D) Gene features for each component. The top 10 genes for each feature were selected. The 12 gene features are annotated using top genes and previous reports as NMF0 Cytotoxic-Feature (F), NMF1 Treg-F, NMF2 Th17-F, NMF3 Naive-F, NMF4 Activation-F (Act-F), NMF5 TregEff/Th2-F, NMF6 Tfh-F, NMF7 IFN-F, NMF8 Central Memory-F, NMF9 Thymic emigrant-F, NMF10 Tissue-F, and NMF11 Th1-F. The Reactome pathway with the smallest p value for each gene feature is shown below the heatmap (Figure S3B; Table S5).

Figure 3

Pan-autoimmune meta-analysis of peripheral CD4⁺ T cells (A) Strategy for meta-analysis of peripheral CD4⁺ T cells across diseases. First, CD4⁺ T cells were extracted from peripheral blood mononuclear cell (PBMC) scRNA-seq datasets using Azimuth. Extracted CD4⁺ T cells were mapped on our reference using Symphony with a batch correction. Mapped cells were used to assess cell frequency and NMF cell features for each cluster. (B) Bar plots showing the number of samples (left) and the number of CD4⁺ T cells (right) enrolled in the meta-analysis. The dashed line in the left plot indicates a sample size of 10. Celiac, celiac disease; Kawasaki, Kawasaki disease; T1D, type 1 diabetes; AD, Alzheimer’s disease; PPP, palmoplantar pustulosis; ParkinsonDis., Parkinson’s disease; BD, Behçet’s disease; CD, Crohn’s disease; Flu, influenza; pSS, primary Sjögren syndrome; RA, rheumatoid arthritis; UC, ulcerative colitis. (C) Dot plot showing changes in cell frequency at cluster L2 resolution. Dot colors show coefficients, and sizes show the significance of the generalized linear model (GLM) (STAR Methods). Detailed statistics can be found in Table S8. Only significant dots (p_adj < 0.05) are shown. (D–F) Principal-component analysis (PCA) plots of samples based on cell frequencies. Sample distributions for each disease state (D), loading vectors for each cell type (E), and sample characteristics in healthy donors (F) are shown. (G) Chord diagram showing the top 100 significant associations with positive coefficients between NMF features and cells in each condition, calculated by GLM (STAR Methods). Detailed statistics are shown in Table S9. The thickness of edges indicates the coefficient of GLM, and colors indicate conditions such as diseases, gender, and age. (H) Strategy for predicting autoimmune states from CD4⁺ T cell profiles using machine learning framework. As the input parameters, one model took only cell frequency, age, and gender (without NMFproj), while the other took cell frequency, NMF cell features in Tcm (Th0) and Tnaive, age, and gender (with NMFproj). (I) Receiver operating characteristic curves of logistic regression models trained by cell frequencies (top left), by NMFproj values in Tnaive and Tcm (Th0) (top right), and both cell frequencies and NMFproj values (bottom). SLE, COVID-19, and MS patients were trained on 159, 116, and 35 patients with the same number of healthy subjects, and evaluated on 40, 89, and 17 patients and the same number of healthy subjects from independent datasets. Numbers in parentheses indicate the area under the curve (AUC).

Figure 4

Partitioned heritability of autoimmune diseases by CD4⁺ T cell features (A) Bar plot showing maximum –log₁₀(q_ E score) among NMF gene features. Partitioned heritability of GWASs was measured using the sc-linker framework. Enrichment of each category is the following: autoimmune diseases, p = 8.51 × 10⁻¹²; inflammatory traits, p = 0.131; and blood cell count, p = 7.59 × 10⁻⁸ (two-sided Mann-Whitney U test). (B) Dot plot showing enrichment of partitioned heritability of autoimmune diseases across NMF gene features. The dashed boxes indicate the factor with the highest E score for each disease. Duplicated traits were removed for the visualization. Full statistics are shown in Table S11.

Figure 5

Relationship between genetic factors and phenotypic changes in CD4⁺ T cells (A) Model of genetic effect on phenotypic changes in CD4⁺ T cells. CD4⁺ T cell changes are observed as qualitative (NMFproj cell features) and quantitative (cell-type frequencies) changes. (B) Scatterplot showing the genetic effect on cell frequencies (x axis) and NMF features (y axis). Sc-linker weight per cell was calculated by dot products of sc-linker outcome (NMF) and NMF cell features. For cell frequencies and NMF cell features, coefficients of GLM output for each cluster L2 population were used. Spearman’s correlation of sc-linker weight and cell frequency/NMF cell feature changes were calculated. For the correlation with sc-linker and NMF cell feature changes, we used the maximum R among NMF features for the visualization. COVID19-A, very severe respiratory symptom; COVID19-B, hospitalized; COVID19-C, SARS-CoV-2 infection. (C) Individual sc-linker weights, cell frequency changes (coefficient for each cluster L2), and NMF cell feature changes in the factor with the highest magnitude (E score) (coefficient for each cluster L2) of MS, MG, and SLE were visualized on the UMAP embeddings (left panel). For the coefficient of the NMF cell feature changes, only one representative factor with the highest E score for each disease is shown. The bar plot of Spearman’s correlation of cell frequency and NMFproj changes with partitioned heritability is shown in the right panel. The colors of the bars, except for cell frequency, indicate the E score calculated using sc-linker.