Single-cell BCR and transcriptome analysis reveals peripheral immune signatures in patients with thyroid-associated ophthalmopathy

Study design

Blood samples from 6 TAO patients and 3 NCs were used to conduct scRNA-seq and scBCR-seq. The baseline information of these samples was listed in Table 1. 15 inactive TAO patients, 7 active TAO patients, and 10 NCs were used for flowcytometry analysis. The basic characteristics of these samples were listed in Table 2. All patients were diagnosed based on Bartley criteria [25]. The patients with other autoimmune diseases, accepted glucocorticoid and/or immunosuppressor treatments were excluded. Healthy adults without autoimmune disease history or infection within the last 3 months were recruited as the NCs and were matched by age and sex. Disease severity was evaluated based on the NOSPECS classification [26].

Peripheral blood monocytes (PBMCs) isolation

The PBMCs were isolated by using a density gradient centrifugation method. Briefly, peripheral blood samples were transferred to a 15mL centrifuge tube and diluted with 1xPBS of equal volume. Equal volume of human peripheral blood lymphocyte separation (Ficoll-Paque, P8900, Solarbio, China) was added to the diluted samples. After mixing, the samples were centrifuged at 700g at 4° C for 30min. The PBMCs located between 1xPBS and Ficoll layer were transferred to a new centrifuge tube, and 1xPBS was added to clean the cells. An appropriate volume (200ul) was taken into 1.5mL centrifuge tube, and a triploid volume (600ul) of erythrocyte lysate was added, gently upside down mixed and placed ice lysis for 8 min. After erythrocyte lysis, incubated at 4° C for 4~5min, centrifuged at 300g at 4° C for 10 min immediately after completion and removing the supernatant, and added 1 ml of re-suspension (precooled 1x PBS+0.04% BSA without calcium and magnesium) to fully re-suspended the cells. Then test with cell counter (Countstar Rigel S2), adjusted the concentration according to the test result. The target concentration was 700-1200 cell/μL, the cell activity was greater than 90%, and the agglomeration was less than 15%. Once the final cell concentration and activity were reached, the cells were placed on ice and 10X Genomics single cell transcriptome chip on-board experiment was carried out within 30 minutes.

scRNA-seq library preparation and sequencing

We prepared single cell RNA-seq libraries with Chromium Next GEM Single Cell 3ʹ Reagent Kits v3.1 (1000268) on the Chromium. We prepared single cell suspensions from cultured cell lines and single cells were suspended in PBS containing 0.04% BSA. The cell suspension was loaded onto the Chromium Next GEM Chip G and ran the Chromium Controller to generate single-cell gel beads (Beckman Coulter-SPRI select Reagent Kit-B23318, USA) in the emulsion (GEMs) according to the manufacturer’s recommendation. Chromium Next GEM Chip K Single Cell Kit (48 rxns, 1000286) was used to capture cells. Captured cells were lysed and the released RNA was barcoded through reverse transcription in individual GEMs. Barcoded, full-length cDNA was generated and libraries were constructed according to the performer’s protocol. Chromium Next GEM Single Cell 5' Reagent Kits v2 (Dual Index) (16 rxns, PN1000263), Library Construction Kit (16 rxns, PN1000190) and Dual Index Plate TT Set A (96 rxns, PN3000431) were used to construct libraries. The quality of libraries was assessed by Agilent-High Sensitivity DNA Kit (5067-4626). Sequencing was performed on the Illumina NovaSeq 6000 with a sequencing depth of at least 50,000 reads per cell and 150 bp (PE150) paired-end reads (performed by Biomarker Technologies Corporation, Beijing, China).

Processing raw data and assay from scRNA-seq of 10× genomics

We performed alignment to this amended reference using 10x cellranger 6.0.0, which employs the STAR sequence aligner [27]. The reference genome was the human genome GRCh38. Overall, 93007 cells passed the quality control (Agilent-High Sensitivity DNA Kit -5067-4626). To remove the cells with low quality, cells with gene number over 500 and the ratio of mitochondria lower than 0.20 were maintained, and genes with at least one feature count in more than three cells were used for the following analysis. We determined gene expression counts using unique molecular identifiers (UMIs) for each cell barcode-gene combination. Following alignment, we filtered cell barcodes to identify those which contain cells using the approach implemented in cellranger 3.0.2, and only these barcodes were considered for downstream analysis including clustering and cell type identification, differential expression analysis by Seurat (version 4.0.1).

Clustering and cell type identification

We applied Louvain community detection [28] to the nearest neighbor graph constructed in PCA space to define a cluster partition. To infer cell types, we trained a neural network classifier to predict cell ontology classes given single-cell RNA-seq mRNA abundance profiles. First, we use singleR (4.0.4) which contains 7 data sets to automate the identification of cell types. In addition to these annotations, we manually added cell state annotations to SingleR to provide a level of granularity below cell ontology classes. After completing the annotations, we further corrected the cell types using the CellMarker database (http://bio-bigdata.hrbmu.edu.cn/CellMarker1.0/).

Gene ontology enrichment analysis, cell cycle scoring

We used Enrichr [29] to perform gene set enrichment analysis against the Gene Ontology Biological Process 2018 version gene set collection. We also used the MSigDB Hallmark gene sets [30], for which we computed enrichment scores using Fisher’s exact test. In both cases, we corrected for multiple hypothesis testing using the Benjamini-Hochberg procedure. We estimated cell cycle activity by scoring the expression of a set of S phase–associated and G2/M phase–associated genes, as shown previously [31] and as implemented in Seurat 4.0.1 [28].

Trajectory inference

We used the Monocle 2 (version 2.18.0) to place cells onto pseudotime trajectories (https://github.com/cole-trapnell-lab/monocle-release), which applies reversed graph embedding (RGE), a recently developed machine learning strategy, to accurately reconstruct complex single-cell trajectories. Monocle 2 requires the genes without a priori information that characterize the biological process, the number of cell fates or branch points in the trajectory, or the design of the experiment. Monocle 2 outperforms not only its previous version but also more recently developed methods, producing more accurate, robust trajectories. The Monocle 2 analysis is based on the following three steps. First, genetic screening: monocle looks for genes that change in “interesting” (not just noisy) ways, and uses those genes to construct the data. Second, lower dimensions: once the genes for cell sequencing are selected, Monocle reduces the dimensions of the data. Third, pseudotime cell analysis: by projecting expression data into a lower dimensional space, the differentiation locus between cells was constructed.

ScBCR-seq library preparation and sequencing

GEM generation and barcoding

GEMs are generated by combining barcoded Single Cell 5’ Gel Beads (Chromium Next GEM Single Cell 5' Reagent Kits v2 (Dual Index), 16 rxns, PN100026), a Master Mix with cells, and Partitioning Oil on a microfluidic chip. To achieve single cell resolution, the cells are delivered at a limiting dilution, such that the majority (~90–99%) of generated GEMs contains no cell, while the remainder largely contain a single cell. Immediately following generation of a GEM, the Single Cell 5’ Gel Bead is dissolved and any co-partitioned cell is lysed.

Upon dissolution of the Single Cell 5’ Gel Bead in a GEM, oligonucleotides containing (i) an Illumina® R1 sequence (read 1 sequencing primer), (ii) a 16 nt 10x Barcode, (iii) a 10 nt Unique Molecular Identifier (UMI), and (iv) a 13 nt template switch oligo (TSO) are released and mixed with cell lysate and a Master Mix that contains reverse transcription (RT) reagents and poly(dT) primers.

Incubation of the GEMs then produces barcoded, full-length cDNA from poly-adenylated mRNA. After incubation, the GEMs are broken and the pooled post GEM-RT reaction mixtures are recovered.

GEM-RT clean up

The pooled post GEM-RT reaction mixture contains barcoded first-strand cDNA from poly-adenylated mRNA, as well as leftover biochemical reagents and primers. Silane magnetic beads are used to purify the cDNA from this mixture. Barcoded, full-length V(D)J segments can then be directly enriched from the purified post Silane magnetic beads which are used to remove leftover biochemical reagents and primers from the post GEM reaction mixture. Full-length, barcoded cDNA is then amplified by PCR to generate sufficient mass for library construction.

cDNA amplification

Barcoded, full-length cDNA is amplified via PCR with primers against common 5’ and 3’ ends added during GEM-RT. This amplification reaction generates sufficient material to construct multiple libraries from B cell-enriched libraries and 5’ gene expression libraries.

Target enrichment

Barcoded, full-length V(D)J segments are enriched from amplified cDNA via PCR amplification with primers specific to Ig constant regions prior to library construction.

Enriched library construction

Enzymatic Fragmentation and Size Selection are used to generate variable length fragments that collectively span the V(D)J segments of the enriched Ig transcripts prior to library construction. R1 (read 1 primer sequence) is added to the molecules during GEM-RT incubation. P5 is added during Target Enrichment. P7, a sample index and R2 (read 2 primer sequence) are added during library construction via End Repair, A-tailing, Adaptor Ligation and PCR. The final libraries contain the P5 and P7 priming sites used in Illumina® bridge amplification.

Sequencing libraries

The Single Cell V(D)J Reagent Kit (Chromium Single Cell Human TCR Amplification Kit, 16 rxns, PN1000252; Chromium Single Cell Human BCR Amplification Kit, 16 rxns, PN1000253) protocol produces V(D)J enriched and 5’ gene expression Illumina-ready sequencing libraries. A library comprises standard Illumina paired-end constructs which begin and end with P5 and P7. For V(D)J enriched libraries, Read 1 encodes the 16 bp 10x™ Barcode, 10 bp UMI, and 13 bp Switch Oligo, as well as the 5’ end of an enriched transcript. For 5’ gene expression libraries, Read 1 encodes the 16 bp 10x Barcode and 10 bp UMI. Due to Enzymatic Fragmentation, for both libraries Read 2 encodes a random internal fragment of the corresponding insert. Sample index sequences are incorporated as the i7 index read.

5’ gene expression data processing

Processing raw data and assay

We performed alignment to this amended reference using 10x cellranger 6.0.0, which employs the STAR sequence aligner. The reference genome was the human genome GRCh38 or the mouse genome mm10. Overall, 36,000 cells passed the quality control. To remove the cells with low quality, cells with gene number over 500 and the ratio of mitochondria lower than 0.20 were maintained, and genes with at least one feature count in more than three cells were used for the following analysis. We determined gene expression counts using unique molecular identifiers (UMIs) for each cell barcode-gene combination. Following alignment, we filtered cell barcodes to identify those which contain cells using the approach implemented in Cellranger 3.0.2, and only these barcodes were considered for downstream analysis including clustering and cell type identification, differential expression analysis by Seurat (version 4.0.1).

Dimensionality reduction

To enable unsupervised clustering and cell type identification, we perform dimensionality reduction with principal component analysis (PCA) on the combined set of samples for each tissue. To visualize the data, we further reduced the dimensionality of all 26,231 cells using Seurat and used t-SNE to project the cells into 2D space. The steps include: 1) Using the LogNormalize method of the “Normalization” function of the Seurat software to calculate the expression value of genes; 2) PCA (Principal component analysis) analysis was performed using the normalized expression value, within all the PCs, the top 10 PCs were used to do clustering and t-SNE analysis; 3) To find clusters, selecting weighted Shared Nearest Neighbor (SNN) graph-based clustering method. Marker genes for each cluster were identified by the “bimod” (Likelihood-ratio test) with default parameters via the FindAllMarkers function in Seurat (4.0.1). Filter Fold Change>1.5 and FDR<0.1 from the calculation results of FindMarker, and then sort to top10 gene as the marker gene.

V(D)J data processing

V(D)J assembly and annotation

Cell Ranger for V(D)J (https://support.10xgenomics.com/single-cell-vdj/software/overview/welcome) is a set of analysis pipelines that process Chromium single cell 5′ RNA-seq output to assemble, quantify, and annotate paired V(D)J transcript sequences. Cell Ranger workflow starts by demultiplexing the Illumina sequencer's base call files (BCLs) for each flow cell directory into FASTQ files.

Cell barcode and UMI correction

Cell Barcode must be on static list of known cell barcode sequences. The barcode may be 1 mismatch away from the list if the mismatch occurs at a low-quality position, it is then corrected. UMIs that are 1 mismatch away from a higher-count UMI are corrected to that UMI if they share a cell barcode.

Read trimming

Trim known adapter and primer sequences from the 5′ and 3′ ends of reads using the cut adapt tool. This tool uses Smith Waterman alignment and allows for a small number of differences from the expected primer sequences.

Read filtering

Cell Ranger aligns reads to all the V(D)J gene segments included in the reference. Read-pairs that exceed a specified alignment score and include at least one 15bp exact match against at least one of the reference segments are included in the set of reads to be assembled. These mappings are not in full alignments and are only used for filtering reads before assembly.

Cell calling

Even though we see many putative cell barcodes in the data, only a fraction of them correspond to droplets that truly contained a cell. The remaining droplets generate background reads. Cell calling algorithm is to select the barcodes corresponding to droplets that contained cells using a 2-component mixture model. Cell calling is performed independently of V(D)J read filtering and assembly.

Flow cytometry analysis

PBMCs were isolated using Ficoll (GE Healthcare, USA) gradient centrifugation, and were suspended in flow cytometry staining buffer (BD Pharmingen, USA) at a final concentration of 10⁷ cells/ml. After incubation with Fixable Viability Stain, BD Horizon™ Brilliant Stain Buffer, and Fc-blocking antibodies (BD Pharmingen, USA), cells were labeled with antibodies specific for surface markers at 4° C for 30 min. Cells were surface-stained with anti-CD4-FITC, anti-CD3-AF700, anti-CD8-Percp.cy5.5, anti-CD11b-BV510, anti-CD14-BV605, anti-CD16-BV421, anti-CD19-APC, anti-CD45-BV650, anti-CD56-PE (Biolegend, USA). To conduct intracellular cytokine staining, the cells were washed, fixed, permeabilized using the BD Cytofix/Cytoperm™ Fixation/Permeabilization Solution Kit, and then stained using anti-IL-10-APC-R7 antibodies (BD Pharmingen, USA). The results were analyzed using FlowJo software version 10.0.7.

Statistical analysis

The results were shown as the mean ± standard error of the mean depending on normality. Comparisons between two groups were performed using Student’s t tests (normal distribution). Spearman rank correlation analysis was used to determine the statistical correlation between variables. Statistical significance was accepted at p < 0.05. All statistical analyses were performed and images were generated using IBM SPSS Statistics 26 and GraphPad Prism (version 7.0), respectively.

Single-cell transcriptional profiling of peripheral immune cells in patients with TAO

To understand the heterogeneity of single-cell transcriptomes at different stages of the disease, we categorized TAO patients into two groups according to the CAS criteria: those with notable orbital inflammation were designated as active group, while those with minimal clinical symptoms were classified as the inactive group [32]. PBMCs were isolated from the whole blood of each subject for further study.

We performed scRNA-seq on a total of nine samples using the 10x Genomics Chromium platform (Supplementary Figure 1A). After initial quality control, a total of 93,007 PBMCs were obtained for subsequent analysis, comprising 23,051 cells from active TAO, 34797 cells from inactive TAO, and 35,159 cells from normal controls (NCs) (Supplementary Table 1). Graph-based clustering analysis identified 17 distinct clusters. Following manual confirmation of the cell type annotations, we successfully categorized the clusters into five major PBMC cell types. These included various populations of CD4⁺ T cell (CD4⁺), NK (CD56⁺), CD8⁺ T cells (CD8A⁺), B cells (CD19⁺), and myeloid (CD11b⁺) (Figure 1A, 1B and Supplementary Figure 1B and Supplementary Table 2). The presence of these five cell types identified by scRNA-seq was validated by flow cytometry analysis (Figure 1C). Notably, 12 of the 17 cell clusters (70.6%) contained at least 100 cells from multiple individual samples, with a median of 9 individual samples represented within each cluster (Supplementary Figure 1C).

Figure 1. Single-cell profiling of blood immune microenvironment in TAO patients. (A) A total of 5 main PBMC cell types were identified in PBMC by automated annotation using SingleR software. Color scheme denoting clusters was labeled with inferred cell types. (B) The t-SNE feature plots showed expression of marker genes of the five cell types. (C) FACS plots validated that CD4⁺ T cells, CD8⁺ T cells, B cells, and NK cells were present in the peripheral blood. Blood sample from one normal control was used for analysis. (D) Bar plots showing the percentage of five PBMC cell types from TAO patients and NCs. Error bars show SEM (* p<0.05, ** p<0.01, *** p<0.001).

We next sought to examine the variations in the prevalence of these populations according to the disease stage. We calculated the percentage of each cluster identified through scRNA analysis relative to the total number of clusters for each patient sample (Figure 1D and Supplementary Table 3). Our findings revealed a remarkable diversity in the immune profiles across patients, with active TAO showing significant infiltration of B cells (p=0.047) and myeloid cells infiltration (p=0.008). In contrast, NK cells were predominantly observed in inactive TAO (p=0.022) and NCs (p=0.059). Although the difference did not reach statistical significance, both inactive TAO and NCs samples showed an enrichment of CD4⁺ T cells. The proportion of CD8⁺ T cells exhibited minimal changes across there sample groups. These findings emphasized the alterations in PBMC composition in individuals with active TAO, underscoring the importance of gaining a deeper understanding of the diversity of cell subtypes and states that related to disease phases.

Regulatory B cells involve in immune regulation during TAO

To comprehensively understand the heterogeneity of B cells and their role in TAO, we conducted further analysis to identify their distinct subpopulations based on RNA-seq data. Using graph-based clustering, we identified five subpopulations comprising a total of 9493 B cells (Figure 2A). Furthermore, we thoroughly assessed the expression of a comprehensive panel of genes associated with B cell lineage to gain deeper insights into their distinct features. The follicular B cells were characterized by the expression of IGHD (IgD) and CD23 (also known as FCER2). The marginal zone B population exhibited high expression of JCHAIN (IgM) and CD21 [33]. Bregs were identified based on the expression of CD1d, CD5, CD19, CD24, CD38, and CD27 [14]. Additionally, germinal center B cells were recognized as a cluster displaying high levels of CD20 (MS4A1), CD38, and FCRL3 [34]. Lastly, a small number of B cells were identified as multilymphoid progenitor cells, expressing CD45RA (PTPRC), CD34, and CD10 (MME) (Figure 2B). The top ten different expression genes (DEGs) for each B cell subpopulation provided additional evidence for our annotation and supported the distinct identity and functional characteristics of each B cell subcluster (Figure 2C).

Figure 2. scRNA-seq revealed Bregs with immunoregulatory properties in TAO. (A) t-SNE plot showing the annotated B cells captured from Figure 1B. Each dot denotes an individual B cell, color denotes cluster origin. (B) The t-SNE feature plots showed expression of marker genes within individually identified B cell populations. (C) Heatmap of DEGs of each B cell subpopulation provides phenotypic information for individual B cell populations. Expression values are scaled between minimum and maximum expression for each gene across all clusters. (D) Prevalence of each B cell subpopulation within all samples from TAO patients with different activity and normal controls. (E) Representative FACS dot plots showed the proportion of Bregs from a patient with active TAO, a patient with inactive TAO, and a normal control (left). Scatter plot showing the frequency of Bregs in active TAO, inactive TAO, and NC group (right). (F) Trajectory analysis of B cell subsets in NCs and patients with TAO. Upper: distribution of the B cell subpopulations on each of the branches. Lower: the relative proportion of B cell subpopulations in each state. (G) Significantly activated motifs in the differentiation process colored by cell clusters. (H) Heatmap showing upregulated DEGs in active TAO patients’ Bregs. The upregulated GO terms enriched in the DEGs are listed with p-value and gene numbers. (I) Violin plots showing upregulated DEGs of Bregs. The diamond inside the violin plot represents the mean. (J) Heatmap showing the distribution of interaction pairs across the cell types.

We subsequently examined the B cell composition in each sample during different phases of the disease. All B cell clusters were observed in each PBMC sample (Supplementary Figure 2A). Although there was no significant difference between the relative proportions of Bregs for the three groups, the proportions of Bregs were lower in active TAO when compared to inactive TAO and normal control groups (p=0.078 and p=0.065, respectively) (Figure 2D). Despite the observed differences not reaching statistical significance, this result indicated a potential trend in the data. To validate these findings, we conducted flow cytometry analysis and used CD19⁺IL-10⁺ as the marker for identifying Bregs cell population. IL-10 had been widely recognized as the hallmark of Breg due to its ability to suppress pro-inflammatory responses [35]. This analysis confirmed that Bregs from active TAO exhibited a reduced frequency compared to those in inactive TAO and the NC group (Figure 2E).

Trajectory analysis was performed to elucidate the dynamic relationships among different B cell subgroups. The analysis uncovered 3 branchpoints and 5 distinct states. State 1 was characterized by a simultaneous distribution of marginal zone B cells and FCRL3 high GC B cells. State 2 consisted of a small fraction of follicular B cells. State 3 predominantly comprised follicular B cells. State 4 and state 5 were primarily composed of Bregs. Additionally, all states exhibited diffused multilymphoid progenitor cells (Figure 2F and Supplementary Figure 2B). In line with previous findings, we observed a reduced cell density of Bregs in state 4 and state 5 in active TAO. According to the SCENIC analysis, a distinct set of genes that regulated by different transcription factors was found to be specifically expressed in Bregs while being almost absent in other cell types (Figure 2G). In Bregs, the genes regulated by RORA, CEBPD, PRDM1, FOSL2, and EOMES displayed a notably higher expression level. RORA, a crucial regulator in cell differentiation and immunity [36], had been recently implicated in B cell development [37]. CEBPD, known to modulate various biological processes, including cell differentiation, proliferation, and cell death, played a pivotal role in regulating immune and inflammatory responses [38]. The DNA-binding protein Blimp-1, encoded by PRDM1, dynamically influences the generation, differentiation, and IL-10 production of Bregs [39]. FOSL2 (also known as FRA2) acts as a key enhancer in regulating early B cell proliferation and differentiation [40]. While EOMES had been recognized as essential for the differentiation of T cells and NK cells [41, 42], its role in B cell development remains to be elucidated.

The GO enrichment analysis revealed that Bregs exhibited upregulated genes associated with immunoregulation (AIF1, CSK, HLA-DRB5, CYBA), cell cycle progression, signal transduction, and apoptosis (ARPC1B, CORO1A, YWHAB, RHOA, CLIC1, S100A6, NME2) (Figure 2H, 2I). Notably, in active TAO, the upregulated genes of Bregs were enriched in inflammation-related GO terms and cellular process-related GO terms. These findings strongly suggested that Bregs may play a crucial role in modulating inflammation and immune responses during the TAO process. According to the CellPhoneDB analysis, Bregs exhibited a high frequency of communication with myeloid cells (Figure 2J). This finding indicates that Bregs and myeloid cells are likely to engage in coordinated immunoregulatory processes during TAO.

scBCR-seq analysis uncovered a notable increase in the diversity of BCRs in active TAO

To allow for a comprehensive analysis of both gene expression and BCR diversity in the B cell populations, scBCR-seq was performed for the PBMC samples used in scRNA-seq analysis. In the active TAO, the proportion of the top 10 clonotypes (0.17%) exhibited a significant decrease compared to in the inactive group (0.37%, p=0.010) and the NC group (0.41%, p=9.31E-05) (Figure 3A). However, analysis of the length of the complementarity determining region 3 (CDR3) did not show any noticeable changes (Figure 3B). When examining the diversity of CDR3, as quantified by the Chao1 index, a significant increase was observed in the active TAO group compared to both the NC group (p=0.016) and the inactive TAO group (p=0.002) (Figure 3C). A similar trend was noted in Breg cells (Figure 3D), with a significant increase compared to the NC group (p=0.011) but no statistically significant difference between the active TAO group and the NC group (p=0.213).

Figure 3. scBCR-seq revealed increased diversity of BCRs in active TAO. (A) The box plot showed the top 10 high-frequency B-cell clonotypes for each sample in active TAO, inactive TAO, and NC group. (B) The distribution of the length of CDR3 amino acid sequence in active TAO, inactive TAO, and NC group. (C) Box plots of Chao1 index for each sample were used to compare the diversity of CDR3of B cells (* p<0.05, ** p<0.01). (D) Box plots of Chao1 index for each sample was used to compare the diversity of CDR3 of Bregs (* p<0.05, ** p<0.01).

These findings suggested that while the frequency of the top 10 clonotypes was reduced in active TAO, there was an increased diversity of CDR3 sequences in this condition, which might have implications for the regulation of immune responses and contribute to the pathogenesis of active TAO.

Alteration in configuration of myeloid cells in active TAO

To assess the involvement of myeloid cells in TAO, we identified the major subtypes of myeloid cells and examined their configuration in active TAO. The analysis revealed five distinct myeloid cell populations, namely monocytes, megakaryocytes, dendritic cells (DCs), monocyte-derived DCs, and erythrocytes (Figure 4A and Supplementary Figure 3A). The identification of these cell types was based on differential gene expression and examination of established lineage markers (Figure 4B and Supplementary Figure 3B). A comparison between active TAO and inactive TAO unveiled notable differences in the proportions of monocytes and DCs. Active TAO patients displayed a significantly higher proportion of monocytes (p=0.006), while the proportion of DCs in active TAO was significantly lower compared to the inactive TAO group (p=0.003) (Figure 4C). These findings suggested that increased levels of monocytes and decreased populations of DCs might be characteristic of the active disease state in TAO patients. To further understand the clinical relevance of quantitative changes of Bregs, DCs, and monocytes in TAO, we examined correlations between the proportion of these cells and TAO disease activity and severity. As shown in Supplementary Figure 4, the proportion of monocytes were positively correlated with the CAS score of TAO patients (p=0.0004, r²=0.9668) (r² represented the coefficient of determination in Spearman rank correlation analysis). And the proportion of Bregs were negatively correlated with the NOSPECS score of TAO patients (p=0.0158, r²=0.8021). However, no significant correlations were found between the proportion of Bregs, DCs and the CAS score of TAO. Similarly, the proportion of DCs and monocytes showed no correlations with NOSPECS score of TAO (Supplementary Table 4 and Supplementary Figure 4).

Figure 4. Analysis of the myeloid landscape reveals decreased DCs and increased monocytes in active TAO. (A) Sub-clustering of myeloid cells identifies 5 cellular categories, including monocytes, megakaryocytes, monocyte-derived dendritic cells, dendritic cells, and erythrocytes. (B) tSNE feature plot representation of marker genes expression within individually identified myeloid cells. (C) Bar plots showing the percentage of myeloid subclusters in active TAO, inactive TAO, and NC group (* p<0.05, ** p<0.01).

DCs might involve in the inflammatory process during active TAO

DCs were subdivided into three distinct subclusters. Cluster 0 was designated as cDC1s, which was characterized by high expression of CDKN1C, FCGR3A, SMIM25, IFITM3, LST1, and MS4A7, and low expression of CD1C and CD141. Clusters 1 and 2 were identified as cDC2s, which exhibited elevated expression levels of inflammatory genes such as CD14, S100A8, GNLY, NKG7, IL32, and S100A12 (Figure 5A, 5B). The proportion of cDC1s was significantly higher in active TAO, whereas cDC2s were more prevalent in inactive TAO and NCs (Figure 5C). Given that cDC1s and cDC2s perform distinct roles in immune regulation – with cDC1 being specialized in inducing Th2 response and cDC2 subclusters exhibiting a more potent ability to induce Th1 response [42], we further examined their gene expression profiles. Regarding cDC1s, the genes associated with MHC II antigen presentation (HLA-DRB5, IFI30), immune activation (NAPSB), and developmental and neurogenetic diseases (NBPF14) were highly expressed in active TAO compared to NCs (Figure 5D, 5E). For cDC2s, higher levels of HLA-DRB5 and FCN1 were observed in active TAO compared to NCs (Figure 5D, 5E). Moreover, for both cDC1 and cDC2, most of the related genes, including HLA-DRB5, NAPSB, NBPF14, IFI30, and FCN1, also exhibited relatively higher expression levels in active TAO compared to inactive TAO (Supplementary Figure 5A and Figure 5E). Enrichment analysis revealed that for both DCs the upregulated genes in active TAO were primarily associated with various processes related to inflammation, immune activation pathways, and complement activation pathways (Figure 5F). Interestingly, the downregulated genes were predominantly associated with Toll-like receptor (TLR) cascade pathways which are essential for DCs activation and maturation [43, 44]. These results suggested that the activation of DCs might be impaired in active TAO, and the observed inflammation during this active phase could be closely related to these cells.

Figure 5. Activation of DCs may be associated with the inflammation during active TAO. (A) tSNE plot of subgroups of DCs. (B) Heatmap of DEGs between three different DC subgroups. (C) Bar chart showing the relative proportion of CD1C^-CD141^- DCs and CD1C⁺B DCs from active TAO, inactive TAO and NCs (* p<0.05). (D) Volcano plots showing DEGs of CD1C^-CD141^- DCs (left), and CD1C⁺B DCs (right) from active TAO patients and NCs. (E) Violin plots showing the expression of HLA-DRB5, IFI30, NAPSB, FCN1, MT-ATP6, FTL, SERPINA1, LTA4H, and NBPF14 in CD1C^-CD141^- DCs and CD1C⁺B DCs from active TAO patients and NCs. (F) Representative Reactome terms enriched by upregulated (left) and downregulated (right) DEGs of DCs in active TAO.

Monocytes play a role in active TAO by regulating inflammatory cytokine production

From a total of 12,959 monocytes, our clustering analysis successfully identified three distinct subpopulations. The first cluster was characterized as classical monocytes (CMOs), which exhibited high expression levels of CD14, S100A8, and S100A9. The second cluster was identified as non-classical monocytes (NMOs), characterized by their high expression of FCGR3A (also known as CD16). The third cluster was recognized as intermediate monocytes (IMOs), with both CD14 and FCGR3A highly expressed (Figure 6A, 6B and Supplementary Figure 5B). To validate these monocyte subclusters, flow cytometry analysis was conducted, confirming CMOs as CD14⁺⁺CD16^-, NMOs as CD14⁺CD16⁺⁺, and IMOs as CD14⁺⁺CD16⁺ (Figure 6C).

Figure 6. Monocytes displayed increasing inflammatory cytokine production and cos-stimulatory molecule expression in active TAO. (A) tSNE plot showing clusters of monocyte subsets. (B) Pseudotime analysis of monocytes recapitulate known lineage relationships, with CMOs (CD14⁺S100A8⁺S100A9⁺) branching into either NMOs (FCGR3A⁺) or IMOs (CD14⁺⁺FCGR3A⁺). (C) FACS plots validated that CMOs, IMOs, and NMOs were present in the peripheral blood. Blood sample from one normal control was used for analysis. (D) Volcano plots showing DEGs of CMOs, IMOs and NMOs in active TAO/NC comparison. (E) Inflammatory cytokines and costimulatory molecules (CD86 and TNF) increased in CMOs and NMOs in active TAO. (F) Representative GO terms enriched by upregulated DEGs of monocytes in active TAO.

Further comparisons were made to assess the proportions of the three monocyte subsets among different groups. Although not statistically significant, it was observed that NMOs and IMOs were more abundant in inactive TAO when compared to active TAO (Supplementary Figure 5C, 5D). Trajectory inference analysis revealed that monocytes followed a developmental trajectory starting from CMOs and branching towards NMOs and IMOs states (Figure 6B). These observations of developmental state alterations in monocyte cells during TAO suggest their potential roles in the disease.

The results of the DEGs analysis revealed that CMOs exhibited significantly higher expression of genes such as HLA-DRB5, NAPSB, and FCN1 in active TAO compared to other groups. Similarly, NMOs in active TAO showed relatively higher expression of genes such as HLA-DRB5 and SNORP3A (Figure 6D and Supplementary Figure 5E). Both CMOs and NMOs exhibited higher expression of TNF and CD86 in active TAO, but these genes were significantly lower in expression in IMOs compared to the normal control group (Figure 6E). Enrichment analysis of DEGs suggested that inflammation-associated pathways and MHC II antigen presentation were upregulated in CMOs and NMOs (Figure 6F). These data indicated that monocytes play a role in regulating pro-inflammatory effects during active TAO.

Crosstalk among Bregs and myeloid cells in TAO

Cell interactions among Bregs, DCs, and monocytes in the context of active TAO were further examined. This analysis identified a total of 1,352 pairs of significant interactions. Remarkably, Bregs were found to exhibit the highest number of ligands, whereas cDC1 subclusters showed the highest number of receptors. Notably, Bregs displayed an increased number of ligands in active TAO when compared to both inactive TAO and the normal control group (NCs) (Figure 7A).

Figure 7. CellChat analysis based on scRNA-seq showed the crosstalk between Bregs and myeloid cells in the immune network. (A) Heatmap showed the distribution of interaction pairs across the cell types. (B) Chord plots showed the interactions of ligand/receptor pairs between active TAO, inactive TAO, and NC group. (C) Violin plots showed the expression levels of LTA and TNF family cytokines in each immune cell type for both TAO and HC groups.

We conducted a detailed analysis of the expression patterns of corresponding receptors and ligands associated with inflammation and immune activation. Compared to inactive TAO, and normal control (NC) groups, it was observed that the CD22-PTPRC interaction between Bregs and cDC1s, cDC2s, IMOs, NMOs, or CMOs were significantly stronger in active TAO patients (Figure 7B). Moreover, we observed robust inhibitory interactions in active TAO, such as BTLA-TNFRSF14 between Bregs and cDC1s, cDC2s, or IMOs. This finding suggests a potentially more immunosuppressive role for Bregs in the context of active TAO. (Figure 7B).

Through examination of the expression levels of genes having potential role in immunosuppressive and inflammatory processes, such as CD22, PTPRC, BTLA, and TNFRSF14 in these cell subclusters, we found that CD22 expression was increased in Bregs during active TAO. Furthermore, the expression levels of BTLA and TNFRSF14 showed an overall increase in Bregs, DCs, and monocytes in active TAO when compared to inactive TAO and the normal control (NC) group (Figure 7C).

These findings suggest intricate regulatory interactions between Bregs, DCs, and monocytes during active TAO, which may contribute to the inflammatory and immune responses observed in the disease.

We know that due to limited understanding of the potential immune aberrations of TAO, it is difficult to develop new and effective therapies and interventions. In the above study, the immune components in peripheral blood of TAO patients and healthy subjects were analyzed by single-cell RNA sequencing and B-cell receptor profile (BCR), and it was found that the proportion of regulatory B cells (Bregs) and type 2 conventional dendritic cells (dc) decreased significantly during the active phase of TAO patients. In contrast, the proportion of type 1 dc is significantly increased, and the transition of B cells to the Breg phenotype may be impaired. In other words, the inflammatory process of TAO is related to the dysfunction of immune regulation. At the same time, abnormalities in the peripheral immune system are primarily driven by enhanced interactions between Breg, dc, and monocytes. There is a link between these identified immune aberrations and potential therapeutic targets, and these findings could provide new therapeutic approaches or interventions for TAO.

These results indicate that some immune cell populations change during TAO activity. The expression level of inflammatory genes is increased, and different inflammatory genes play different roles in immune regulation, that is to say, the change of immune cell population plays a positive role in immune regulation. In addition, there are complex regulatory interactions between different cells during TAO activity, which may contribute to the inflammatory and immune responses observed in disease. In addition, dynamic changes in the activity of key regulators suggest that TAO activation appears to involve inflammation and immune dysfunction. Overall, our findings provide comprehensive insights into the molecular regulation and cell remodeling of the active phase of TAO at the single-cell level, with positive clinical implications for suggesting the potential role of monocytes in disease through the observed developmental status of monocytes during TAO.

Data limitations and perspectives

The present study contains a small sample size and lack of functionally validated study of immunosuppressive function of Bregs in the evolution of TAO. However, our study highlighted significant changes in the proportion and transcriptional heterogeneity of B cells, DCs, and monocytes in individuals with active TAO, inactive TAO, and NCs. These findings provided valuable insights into the alterations in immune cell infiltration within the peripheral immune environment during active TAO. Notably, we identified potential inhibitory effects of two ligand-receptor pairs, CD22-PTPRC and BTLA-TNFRSF14, which may contribute to the decrease of Bregs and impaired inhibition of inflammation in active TAO. These potential pathogenic cells and molecules may account for the inflammatory changes of TAO, and, thus, may be identified as novel therapeutic targets for this disease.

Additionally, we described the immune characteristics of BCR repertoires in TAO, revealing a significant increase in BCR diversity during active TAO.

Our observations shed light on the molecular and cellular basis of the inflammatory reaction observed in active TAO. By gaining a better understanding of the immune dysregulation in this condition, our study lays the groundwork for the development of more effective immunotherapy strategies in the future. With further validation and exploration of the underlying mechanisms, these discoveries hold the potential to contribute significantly to the advancement of targeted and personalized treatments for individuals with TAO.