Single-cell transcriptomic analysis in two patients with rare systemic autoinflammatory diseases treated with anti-TNF therapy

Systemic autoinflammatory diseases (SAIDs) are a group of rare diseases characterized by recurrent or continuous inflammation, typically accompanied by genetic variants. Good responses to anti-TNF therapy were observed in SAIDs patients. However, the mechanisms underlying the disease flare and the response to TNF blocking therapy have not been fully elucidated. Here, single-cell RNA sequencing technology was used to describe the transcriptomic profile of PBMCs and PMNs in two SAID patients both before and after anti-TNF treatment. Interferon responses were involved in the disease flare. After anti-TNF therapy, clinical symptoms were alleviated while TNF and IL-1 were unexpectedly increased, indicating that these inflammatory cytokines are not positively correlated with disease activity. Trajectory analysis showed that inhibition of macrophage differentiation, rather than reduction of the inflammatory cytokines, as the potential mechanism of anti-TNF treatment response in SAIDs.

1 Introduction

Systemic autoinflammatory diseases (SAIDs) are a group of rare diseases caused by defects or dysregulation of the innate immune system, characterized by recurrent or continuous inflammation (e.g. fever, serosal, synovial, or cutaneous inflammation) and by the lack of a primary pathogenic role of the adaptive immune system (autoreactive T-cells or autoantibodies production) (1–4). This group of diseases includes more than 40 monogenic, as well as several polygenic diseases (5, 6). Classic monogenic SAIDs include familial Mediterranean fever (FMF, caused by MEFV gene variants), TNF receptor-associated periodic fever syndrome (TRAPS, by autosomal dominant variants in the TNFRSF1A gene), cryopyrin-associated periodic syndromes (CAPS)/NLRP3-associated autoinflammatory disease (NLRP3-AID, by gain-of-function variants in the NLRP3 gene), and hyper-IgD syndrome (HIDS)/mevalonate kinase deficiency (MKD, by autosomal recessive variants in the MVK gene) (6). Patients with high disease activity can present with recurrent attacks of fever and other inflammatory symptoms, and a successful treatment can reduce disease activity, as manifested by reduced frequency or severity of disease flare or no flare. Regarding to pathogenesis, these classic SAIDs are classified as IL-1 mediated SAIDs, as the variants can directly or indirectly promote the assembly of inflammasomes that proteolytically process the inactive pro-IL-1β into active IL-1β. This can occur by increased intracellular sensor/pattern-recognition receptors (PRR) function (e.g. CAPS, FMF), generation of intracellular stress (e.g. TRAPS, HIDS), or loss of a negative regulator (e.g. Deficiency of IL-1 receptor antagonist, DIRA) (6). A lot of studies focus on the pathogenic role of the inflammasome – IL-1 axis in these SAIDs, and many patients have good responses to IL-1-blocking agents (e.g. anakinra, canakinumab, and rilonacept), although some patients with FMF, TRAPS or HIDS are less responsive to IL-1 inhibition (7, 8; 9–15). Little is known about other pathways in the pathogenesis of SAIDs, and the treatment response to other biological agents, e.g. TNF or IL-6 blocking agents. Some studies have reported increased NFκB activation and defective autophagy as potential pathogenic mechanisms (16–18). In our previous study, we observed a good response to TNF inhibitors in CAPS patients (19). Whether there are other pathways involved in the pathogenesis of SAIDs, and what the underlying mechanisms are in response to TNF blocking therapy, remain to be elucidated.

Single-cell RNA sequencing (scRNAseq) is a powerful technique to detect transcriptomics of single cells, which can further reveal cellular heterogeneity, gene regulatory networks, cell-cell interactions, and differentiation trajectories in tissues. This technique has been widely used in oncology, microbiology, neurology, and immunology (20). Yet, it has never been applied to SAIDs research. In this study, we are the first to apply this technique to classic SAIDs, aiming to identify the potential mechanisms of disease flare other than the IL-1 pathway and the response to anti-TNF therapy, and to provide the first transcriptomic resource of SAIDs spectrum diseases.

2 Results

2.1 Single cell profiling of PBMCs and PMNs in SAID patients

This study enrolled 2 patients with classic SAIDs. One patient was a 20-year-old Chinese woman diagnosed as CAPS, with NLRP3 D303G variation, and another patient was a 16-year-old Chinese man diagnosed as TRAPS, with TNFRSF1A V202D variation. Detailed information can be found in Supplemental Data S1. Both patients had high disease activity during the first visit (pre-treatment). The CAPS patient presented urticarial rash, arthritis (interphalangeal joint and knees), lymphadenopathy, hearing loss, and vision loss, while the TRAPS patient was in the intermittent phase of disease episodes and did not show significant symptoms. Both patients received Etanercept treatment, and during the 3-month follow-up visit (on-treatment), both patients reported significant remission of clinical symptoms and no disease relapse.

Peripheral blood mononuclear cells (PBMCs) and polymorphonuclear neutrophils (PMNs) were extracted from pre-treatment (CAPS_NT, TRAPS_NT) and on-treatment (CAPS_TNFi, TRAPS_TNFi) blood samples of both patients, as well as the blood of one healthy donor (HD). Then the cells were subjected to 3’-scRNAseq (10X Genomics) to study the transcriptomic profiles (Figure 1A). As the average detected genes (N_gene) in different cell types varied dramatically (e.g. neutrophils and platelets had very low N_gene, while monocytes and plasma cells had higher N_gene, Figure 1F), we used different quality filtering criteria for each cluster (see Materials and methods) and acquired 30,377 cells for the downstream analysis. To identify the common clusters of all 5 samples, Harmony algorithm (21) was used to remove the batch effect, and identified 17 clusters by unsupervised clustering (Figures 1B, C), including two B-cell clusters (memory, naïve), plasma cells, five T-cell clusters (memory CD4, naïve CD4, memory CD8, naïve CD8, effector CD8), two NK-cell subtypes (CD16⁺ and CD16^-), two monocyte subtypes (CD14⁺, CD16⁺), DCs, pDCs, neutrophils, platelets, and progenitors (Figure 1B). Representative markers of each cell type were shown in Figure 1D. T-cells, CD14 monocytes, CD16 NK-cells, and neutrophils were the more frequent cell types in the blood (Figure 1E).

FIGURE 1

Figure 1 Single Cell Profiling of PBMCs and PMNs in SAID Patients. (A) Study design, comparing transcriptomics of PBMCs and PMNs from pre-treatment (CAPS_NT, TRAPS_NT) and Etanercept treated (CAPS_TNFi, TRAPS_TNFi) blood samples of both patients and one healthy donor (HD). (B, C) UMAP plot, colored by 17 subtypes (B) identified by unsupervised clustering or 5 sample origins (C). Batch effect removed by Harmony algorithm. (D) Expression of representative markers of each cell type (y-axis) in 17 clusters (x-axis). Dot size represents the percentage of cells in which the gene is detected. Color indicates the centered mean expression. (E) Fraction of 5 sample origins in 17 cell types. Total cell count of each sample was normalized to 100. (F) Boxplot displaying the detected gene number (nFeature_RNA) and percentage of mitochondria genes (percent.mt) in each cell type.

2.2 Gene expression profile alters significantly after anti-TNF therapy in CAPS patient

After identifying common clusters, we investigated the differences among these samples in each cell type. Unaligned uniform manifold approximation and projection (UMAP) was plotted based on PCA without Harmony integration and it intuitively showed marked changes in the gene expression profile of CAPS after treatment, compared to other samples (Figure 2A). This may suggest that the transcriptome is more associated with disease attack rather than disease activity, as both patients had high activity but only the CAPS patient was at disease flare. To make the comparison more quantitative, we calculated the Pearson Correlation of average gene expression of each cluster between groups (Figure 2B). Note that low N_gene and cell counts can usually generate lower Pearson Correlation value due to a higher variance, like neutrophils, platelets, DCs/pDCs and CD16^- NK cells. In general, CAPS_TNFi shows a lower correlation value compared to other samples, which confirmed our observation in Figure 2A. B cells, effector and memory CD8 T cells, NK cells and monocytes had relatively higher variance than naïve or CD4 T cells, indicating that these can be the main cell types that respond to anti-TNF therapy in the CAPS patient.

FIGURE 2

Figure 2 Gene Expression Profile Alters Significantly after Anti-TNF Therapy in CAPS Patient. (A) UMAP plot, colored by 17 subtypes (left) or 5 sample origins (right). Batch effect not corrected. (B) Pearson correlation of all detected gene expression between samples in each cell type. (C) Heatmap showing top 10 up- and down-regulated GO pathways in CAPS_TNFi.

To further look into which pathways contribute to the major changes in CAPS_TNFi, Gene Set Enrichment Analysis (GSEA) was performed based on Gene Ontology (GO) database, and identified common top 10 up- and down-regulated pathways in some representative cell types (Figure 2C). Surprisingly, many up-regulated pathways after anti-TNF therapy in CAPS were related to stress response and innate immune response-activating signal transduction, including pattern recognition receptor (PRR) signaling pathways like toll−like receptor (TLR) and NOD2 signaling pathways. Meanwhile, some metabolic pathways related to aerobic respiration were down-regulated (Figure 2C). To confirm this, SCENIC analysis (22) was performed which inferred the gene regulatory network and predicted regulon (transcription factor, TF) activities, and found that indeed many inflammatory-related TFs are activated in CAPS_TNFi, including canonical NFκB pathway TFs (NFKB1, REL) and AP-1 complex (FOS, JUN) (23) (Figures 3A, B). We reasoned that the upregulation of inflammatory pathways could be a compensatory adjustment to anti-TNF therapy, although this requires further validation.

FIGURE 3

Figure 3 Gene Regulatory Network (GRN) Predicted by SCENIC Analysis. (A) Heatmap showing top 6 most activated TFs in CAPS_TNFi predicted by SCENIC. (B) Gene regulatory network (GRN) of CD16⁺ monocytes, predicted by SCENIC. The gene expression (round nodes) or regulon activity (square node) in CAPS_NT (left) or CAPS_TNFi (right) is shown as node color (z-score).

Although gene expression changes were less pronounced in the TRAPS patient after anti-TNF treatment than in CAPS, we also performed a comprehensive analysis of differentially expressed genes (DEGs), up/downregulated TF activities and GO pathways in the TRAPS patient, which are documented in Supplementary Data S2. The same analysis was performed comparing CAPS_NT/TRAPS_NT to healthy donors to identify the potential pathogenic genes in SAID (Supplementary Data S2). However, this may be difficult due to the small sample size, as some DEGs may be individual-specific (e.g. genes on the Y chromosome, ribosomal proteins) and apparently not associated with the disease. These data may be more useful when more SAID data becomes available in the future.

2.3 Disease flare is unnecessary to be positively correlated with serum inflammatory cytokines

One important clinical question is which markers are elevated during disease attack. Normally, we evaluate clinical manifestations and laboratory tests including white blood cell (WBC) count, erythrocyte sedimentation rate (ESR), and C-reactive protein (CRP), etc (24). Some inflammatory cytokines, like TNF, IL-1β and IL-6 are considered to be the cause of systemic inflammation, although IL-6 has both pro- and anti-inflammatory functions (25). Thus, we asked how these inflammatory cytokines changed after anti-TNF therapy. First, we checked which cell types were the main source of each cytokine (Figure 4A), and found that TNF was mainly secreted by CD16⁺ monocytes, IL-1β mainly by monocytes, DCs and neutrophils, and IL-6 by a small subset of B cells. Next, the expression levels of these cytokines were examined in the relevant clusters. We noticed that in general they were even more highly expressed after anti-TNF treatment (Figure 4B). This appears counterintuitive, while similar to what we found in Figure 2C. Thus, we asked if the increase of inflammatory cytokines after anti-TNF treatment occurred commonly, by chance, or was just a technical effect of scRNAseq. Therefore, we tried to validate this using other methods. As some patients were tested for serum TNF and IL-6 during their visit to the hospital, we checked if these cytokines were increased after anti-TNF treatment at the protein level. We reviewed our SAIDs cohort from 2015 (26), until Jan 2021, and among total of 228 patients, 34 were received anti-TNF therapy (either Etanercept, Infliximab, Golimumab or Adalimumab), and 8 patients were tested for TNF and IL-6 both before and after anti-TNF therapy. The serum TNF levels were increased while IL-6 decreased after the treatment (Figure 4C). Although it is still not clear what role IL-6 plays during the treatment, we confirmed that TNF was increased after anti-TNF therapy, by both scRNAseq and serum protein test. Taken together, the inflammatory cytokines, at least TNF and IL-1β can be up-regulated after anti-TNF therapy and are not necessarily positively correlated with disease flare.

FIGURE 4

Figure 4 Disease Flare Is Unnecessary to Be Positively Correlated with Inflammatory Cytokines. (A) Expression of inflammatory markers TNF, IL-1 and IL-6 of each cell type (y-axis) in 17 clusters (x-axis). Dot size represents the percentage of cells in which the gene is detected. Color indicates the centered mean expression. (B) Violin plots displaying the expression of inflammatory markers in their main source cell types, comparing 5 sample origins. (C) Serum TNF and IL-6 levels in SAID patients, comparing before and after anti-TNF therapy. (D) Bar plot showing the number of up- and down-regulated differentially expressed genes (DEGs) of CAPS_NT of each cell type. (E) Heatmap showing top 10 up-regulated DEGs in CAPS_NT.

Next, we asked what changes reflect the disease attack, if not TNF or IL-1β. As we know that among the 5 samples, the CAPS_NT sample was taken during disease attack while other samples were not. Therefore, the DEGs of CAPS_NT in each cell type were checked (Figures 4D, E). We noted that some interferon-induced genes, including IFI, IFIT and IFITM families, which are involved in broad-spectrum antiviral functions (27) were upregulated in CAPS_NT, remarkably in the CD16⁺ monocyte subtype (Figure 4E). Consistently, by SCENIC analysis we also observed that several interferon regulatory factors, e.g. IRF2, 5 and 7 were activated in the CAPS_NT sample in CD16 monocytes (Figure 3B). These data indicated that the interferon response might play a role in the disease flare.

We further investigated whether the increased interferon response in CAPS_NT was caused by elevated stimuli levels (IFN-α, IFN-γ), or enhanced sensitivity in these monocytes. First, IFN-α and IFN-γ were checked at both RNA and protein levels. In the single-cell RNA-seq data, we did not see increased IFN expression in non-treated samples (Figure S1A). At the protein level, we acquired several serum samples of non-treated SAID patients (CAPS = 4, TRAPS = 1) from the biobank in our hospital as well as 3 healthy donors, and tested IFN-α and IFN-γ levels by ELISA. Consistent with scRNAseq data, there was no increased serum IFN in non-treated SAID samples (Figure S1B). We then checked the IFN signaling pathways in CD16 monocytes. The IFN signaling gene lists from the Reactome database (28) were used and calculated the enrichment score using the AUCell package (22). We found that indeed there was higher activation of IFN signaling pathways in CAPS_NT (Figure S1C). These data supported that the increased interferon response in CAPS_NT was probably caused by enhanced sensitivity to the IFN stimuli in CD16 monocytes, rather than elevated ligand levels.

2.4 Inhibition of macrophage differentiation is a potential mechanism of anti-TNF therapy

From the previous analysis, we also noticed that complement components (C1Qs) were more upregulated in pre-treatment than on-treatment samples in CD16⁺ monocytes (Figure 4E), which is one of the key features of macrophages (29). Interestingly, CD16⁺ monocytes indeed resemble macrophages in many aspects, as they have low expression of classic monocyte markers like CD14 and CCR2, and highly express CX3CR1 which explains why they migrate and adhere more than CD16⁻ monocytes to fractalkine-secreting endothelium, and specialize in complement and FcR-mediated phagocytosis and anti-viral responses (30–32). The lower expression in complement after TNF blockade gave us a hint that the treatment might change the differentiation state of CD16⁺ monocytes. To validate this hypothesis, pseudo-time trajectory analysis was performed within monocyte clusters using Palantir and scVelo/Velocyto (33–35) which could predict the direction of differentiation and “RNA velocity” — the time derivative of the gene expression state — based on the relative abundance of nascent (unspliced) and mature (spliced) mRNA (Figure 5A, B). We found that before treatment, the CD16⁺ monocytes were highly dynamic and differentiating into the terminal state, while after TNF blocking, the differentiation was blocked, depicted by the arrows with shorter lengths than untreated (Figure 5B). Indeed, TNF is known to be involved in macrophage differentiation (36). This analysis shows that inhibition of macrophage differentiation, rather than reducing the overall level of inflammatory cytokines, could be a potential mechanism of anti-TNF treatment response in SAIDs.

FIGURE 5

Figure 5 Differentiation Trajectory Analysis of Monocytes. (A, B) Differentiation trajectory of NT/TNFi monocytes, predicted by Velocyto/ScVelo, based on tSNE embeddings calculated by Palantir. Differentiation direction shown by arrows, in all samples (A) or split by each sample (B).

3 Discussion

In this study, we for the first time described the transcriptomic profile of PBMCs and PMNs in SAIDs patients, and also revealed transcriptional changes after anti-TNF therapy. We observed that after anti-TNF therapy, some pathways related to stress response and innate immune response-activating signal transduction were up-regulated, while some metabolic pathways related to aerobic respiration were down-regulated in CAPS, whereas there were no significant changes in TRAPS. This may suggest that the transcriptome is more associated with disease attack rather than disease activity, and may partially explain the unexpected good response to IL-1 blockade compared with anti-TNF therapy observed in TRAPS patients in previous studies (37, 38). We also identified that inflammatory cytokines such as TNF and IL-1β were up-regulated after anti-TNF therapy and were not necessarily positively correlated with disease activity. Interestingly, some interferon-induced genes were up-regulated in the CAPS pre-treatment sample, which indicated that the interferon response can be correlated with disease flare. Further analysis revealed that one of the potential mechanisms of anti-TNF treatment response in SAIDs may be the inhibition of macrophage differentiation rather than decreasing the overall level of inflammatory cytokines. As macrophages are a major cell type involved in inflammation, and their activation can be uncontrolled in SAIDs due to the genetic defect, inhibition of macrophage differentiation may help reduce the frequency of disease flare.

As SAIDs are rare diseases, in this exploratory study where we only enrolled 2 patients (one CAPS and one TRAPS) to test the feasibility of using scRNAseq for SAIDs research. As the number of samples was limited, we could not gain robust conclusions from this study. However, we can see that after anti-TNF therapy, we observed significant transcriptional changes after treatment in CAPS patient, indicating that scRNAseq can be a useful technique for the future study of SAIDs and treatment response. Meanwhile, we are also aware of the shortcomings of this study design, which can be improved in the future: i) Some differentiated immune cells in the peripheral tissues (e.g. myeloid cells in skin lesions or synovial fluid) may have important pro-inflammatory functions. Therefore, we should consider not only extracting circulating PBMCs/PMNs but also collecting samples from peripheral tissues. ii) Some long-term effects caused by genetic variation or treatment may be the result of epigenetic modifications that cannot be revealed by scRNAseq, which is more reflective of immediate changes in cell state. Novel techniques such as single cell assay for transposase-accessible chromatin (ATAC) together with RNA sequencing may help us better understand such mechanisms. iii) In the study design we included neutrophils as they have important pro-inflammatory functions. However, we then realized that they had very few transcripts detected and less usable sequencing reads, as they had relatively low RNA content and relatively high levels of RNases with other inhibitory compounds. This brought great challenges to the subsequent analysis as we could not get high-quality information from these cells. In future studies, we should consider either supplementing RNase inhibitor to preserve the RNA or only keeping other high-quality cells without neutrophils.

In conclusion, these results may indicate that interferon response, rather than TNF or IL-1, may be involved in the disease flare of the SAIDs. Trajectory analysis showed that inhibition of macrophage differentiation as the potential mechanism of anti-TNF treatment response in SAIDs. Our work provided new insights into the potential mechanisms of disease flare and anti-TNF treatment response, new transcriptomic resources of SAID spectrum diseases, and expanded the application of single-cell RNA sequencing technology to new disease areas.

3.1 Limitations of the study

As an exploratory study, the main limitation of this study is the small sample size, making it impossible to generalize in this situation. In addition, there is a lack of appropriate mouse models to validate our conclusions. The aim of the further research should be to expand the study cohort and apply multi-omics technology (e.g. single-cell RNA-seq, scATAC-seq and sc-proteomics) to various tissue types (e.g. PBMC, skin lesion or synovial fluid) to gain further insight into the pathogenesis of SAIDs, and to validate the findings by genetically engineered mouse models (GEMM) that mimic the clinical phenotypes.

4 Materials and methods

4.1 Patients

This study enrolled 2 patients with classic SAIDs. One patient was a 20-year-old Chinese woman diagnosed as CAPS, with NLRP3 D303G variation, and another patient was a 16-year-old Chinese man diagnosed as TRAPS, with TNFRSF1A V202D variation. Detailed information can be found in Supplemental Data S1. Both patients were in the cohort of adult SAID patients in Peking Union Medical College Hospital (PUMCH), which was described in the previous publication (26). This study was approved by the Institutional Review Board of Peking Union Medical College Hospital and performed according to the Declaration of Helsinki. Informed consents were obtained from all participants.

4.2 Single-cell RNA sequencing

4.2.1 Sample preparation

The processing of blood samples, including PBMC and neutrophil isolation, library preparation, and single-cell RNA sequencing were done by CapitalBio Technology, Beijing. In brief, ACCUSPIN™ System-Histopaque^®-1077 and Histopaque^®-1119 reagent were used to isolate mononuclear cells and granulocytes, respectively, following the protocol on Sigma-Aldrich website. Using single cell 3’ Library and Gel Bead Kit V3 (10x Genomics, 1000075) and Chromium Single Cell B Chip Kit (10x Genomics, 1000074), the cell suspension (300-600 living cells per microliter determined by Count Star) was loaded onto the Chromium single cell controller (10x Genomics) to generate single-cell gel beads in the emulsion according to the manufacturer’s protocol. In short, single cells were suspended in PBS containing 0.04% BSA. About 6,000 cells were added to each channel, and the target cell will be recovered was estimated to be about 3,000 cells. Captured cells were lysed and the released RNA were barcoded through reverse transcription in individual GEMs. Reverse transcription was performed on a S1000TM Touch Thermal Cycler (Bio Rad) at 53°C for 45 min, followed by 85°C for 5 min, and hold at 4°C. The cDNA was generated and then amplified, and quality assessed using an Agilent 4200 (performed by CapitalBio Technology, Beijing).

4.2.2 Data processing and cell clustering

Raw sequencing data (fastq files) were mapped to the human genome (build GRCh38) using CellRanger software (10x Genomics, version 3.0.2). Raw gene expression matrices generated per sample were analyzed with the Seurat package in R (39). To achieve clean cell clustering results, we divided the cell filtering process into two major steps: primary clustering and fine adjustment. Primary clustering: 5 samples (CAPS_NT, CAPS_TNFi, TRAPS_NT, TRAPS_TNFi, HD) were merged together and cells were filtered by nFeature_RNA (genes detected) > 200 and percent.mt (percentage of mitochondria genes) < 12.5%. High variable genes were selected by FindVariableFeatures and auto-scaled by ScaleData function using default parameters, and a principal component analysis (PCA) was performed for all datasets using the default RunPCA function in the Seurat package. Batch effect correction of each sample was done using the Harmony algorithm (21) based on PCA space, followed by FindNeighbors and FindClusters function (dims = 30, resolution = 0.5) in the Seurat package for unsupervised clustering. In total 20 clusters were found. As plasma cells were not correctly identified by unsupervised clustering, they were manually annotated. Fine adjustment: scDblFinder package was used to predict potential doublets in the datasets (40). As neutrophils and platelets naturally have much fewer transcripts than other cell types, and DCs are often misclassified as doublets, we divided all cell types into 3 groups and use different filter criteria for each group. Group1: including Naïve CD4, Naïve CD8, Memory CD4, Memory CD8, Effector CD8, CD16⁺ NK, CD16^- NK, CD14 Mono, Naïve B, and Memory B, these cells were filtered by nFeature_RNA > 1200 and kept only singlets by scDblFinder; Group 2: including Neutrophil and Platelet, these cells were filtered by nFeature_RNA < 800 and kept only singlets by scDblFinder; Group 3, including FCGR3A Mono, DC, and pDC, these cells were filtered only by nFeature_RNA > 1200 regardless of scDblFinder prediction. In addition, 1 RBC cluster, 3 doublet clusters and 1 mitotic cluster were removed as they were not informative. All plasma cells were kept manually. After filtering, we ran the same pipeline as primary clustering mentioned above with slightly changes of several parameter (dims = 20, resolution = 0.4). Final clustering results were shown in Figure 1B.

4.2.3 Differential expression analysis and data visualization

To visualize cells on a 2D plot, Uniform manifold approximation and projection (UMAP) was done by RunUMAP function (dims = 20) in Seurat. Differentially expressed genes (DEGs) were identified by the Wilcoxon Rank Sum test using the FindMarkers function. Gene expression levels or gene set enrichment scores (AUCell score) were shown in t-score or z-score for heatmaps or waterfall plots. UMAP or tSNE plots were done using DimPlot or FeaturePlot functions in Seurat. Heatmaps, modified stacked violin plots, boxplots, and bar plots of cluster proportion were generated using customized codes in R, and these functions were integrated into the R package “SeuratExtend” which is available on Github (https://github.com/huayc09/SeuratExtend).

4.2.4 SCENIC and gene set enrichment analysis

To carry out the transcription factor network inference, SCENIC workflow was performed using Nextflow pipeline (22), and regulon activity of each cell was evaluated using AUCell score with Bioconductor package AUCell. For functional/pathway analysis, gene set lists were collected from databases including Gene Ontology (GO) and Reactome. For Gene Set Enrichment Analysis (GSEA), the enrichment of given gene sets of each cell was evaluated using AUCell package as well. Gene regulatory networks (Figure 3) were plotted using Cytoscape software (41).

4.2.5 Trajectory analysis

For trajectory analysis (Figure 5), monocyte subsets of CAPS and TRAPS samples were extracted and integrated by Harmony algorithm. Python package Palantir (34) was used to calculate diffusion map and diffusion components based on Harmony space, then cells were visualized by tSNE. To predict the differentiation direction, we conducted a Velocyto pipeline (33) using the *.bam file and barcode information generated by CellRanger, and used ScVelo in Python (35) for better visualization.

4.3 ELISA

The plasma levels of IFN-α and IFN-γ were detected according to the instructions of the commercial ELISA kits (EXCELL Bio, China).

4.4 Statistical analysis

Differentially expressed genes (DEGs) were identified by the Wilcoxon Rank Sum test using the FindMarkers function of Seurat package in R. All remaining statistical analysis were performed by the ggpubr package in R with default parameters.

Data availability statement

The raw sequence data reported in this paper have been deposited in the Genome Sequence Archive (42) in National Genomics Data Center (43), China National Center for Bioinformation / Beijing Institute of Genomics, Chinese Academy of Sciences (GSA-Human: HRA003572) that are publicly accessible at https://ngdc.cncb.ac.cn/gsa-human. All software is freely or commercially available. To ensure data accessibility to non-bioinformaticians, we made the an interactive web tool generated by ShinyCell (44) at https://yichao-hua.shinyapps.io/shen_lab_said_anti-tnfa/.

Ethics statement

The studies involving human participants were reviewed and approved by Institutional Review Board of Peking Union Medical College Hospital. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin. Written informed consent was obtained from the individual(s), and minor(s)’ legal guardian/next of kin, for the publication of any potentially identifiable images or data included in this article.

Author contributions

MS and YH conceived the study. YH performed the bioinformatics analyses. NW performed the Elisa experiment. JM summarized the medical history of SAID cohort. YH and MS wrote the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by the Natural Science Foundation of Beijing (Grant No.7192170), the National Key Research and Development Program of China (Grant No.2016YFC0901500; 2016YFC0901501).

Acknowledgments

The authors thank CapitalBio Technology for the technical support of single-cell RNA sequencing, Peking Union Medical College Hospital Laboratory Department for testing serum cytokines, Dr. Junbin Qian (Zhejiang University) for the suggestions of single cell analysis, Dr. Kathryn Jacobs (VIB-KU Leuven) for the help with language modification.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1091336/full#supplementary-material

Supplementary Data Sheet 2 | Differentially expressed genes, regulon activities and GSEA

References

1. Centola M, Aksentijevich I, Kastner DL. The hereditary periodic fever syndromes: molecular analysis of a new family of inflammatory diseases. Hum Mol Genet Oxford Univ Press (OUP) (1998) 7(10):1581–8. doi: 10.1093/hmg/7.10.1581

CrossRef Full Text | Google Scholar

2. Manthiram K, Zhou Q, Aksentijevich I, Kastner DL. The monogenic autoinflammatory diseases define new pathways in human innate immunity and inflammation. Nat Immunol (2017) 18(8):832–42. doi: 10.1038/ni.3777

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Ben-Chetrit E, Gattorno M, Gul A, Kastner DL, Lachmann HJ, Touitou I, et al. Consensus proposal for taxonomy and definition of the autoinflammatory diseases (AIDs): a Delphi study. Ann rheumatic diseases BMJ (2018) 77(11):1558–65. doi: 10.1136/annrheumdis-2017-212515

CrossRef Full Text | Google Scholar

4. Nigrovic PA, Lee PY, Hoffman HM. Monogenic autoinflammatory disorders: Conceptual overview, phenotype, and clinical approach. J Allergy Clin Immunol (2020) 146(5):925–37. doi: 10.1016/j.jaci.2020.08.017

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Doria A, Zen M, Bettio S, Gatto M, Bassi N, Nalotto L, et al. Autoinflammation and autoimmunity: bridging the divide. Autoimmun Rev (2012) 12(1):22–30. doi: 10.1016/j.autrev.2012.07.018

PubMed Abstract | CrossRef Full Text | Google Scholar

6. de Jesus AA, Canna SW, Liu Y, Goldbach-Mansky R. Molecular mechanisms in genetically defined autoinflammatory diseases: disorders of amplified danger signaling. Annu Rev Immunol (2015) 33:823–74. doi: 10.1146/annurev-immunol-032414-112227

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Hoffman HM, Throne ML, Amar NJ, Sebai M, Kivitz AJ, Kavanaugh A, et al. Efficacy and safety of rilonacept (interleukin-1 trap) in patients with cryopyrin-associated periodic syndromes: results from two sequential placebo-controlled studies. Arthritis rheumatism (2008) 58(8):2443–52. doi: 10.1002/art.23687

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Hoffman HM, Throne ML, Amar NJ, Cartwright RC, Kivitz AJ, Soo Y, et al. “Long-term efficacy and safety profile of rilonacept in the treatment of cryopryin-associated periodic syndromes: results of a 72-week open-label extension study. Clin Ther (2012) 34(10):2091–103. doi: 10.1016/j.clinthera.2012.09.009

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Kuemmerle-Deschner JB, Ramos E, Blank N, Roesler J, Felix SD, Jung T, et al. Canakinumab (ACZ885, a fully human IgG1 anti-IL-1β mAb) induces sustained remission in pediatric patients with cryopyrin-associated periodic syndrome (CAPS). Arthritis Res Ther (2011) 13(1):R34. doi: 10.1186/ar3266

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Kuemmerle-Deschner JB, Tyrrell PN, Koetter I, Wittkowski H, Bialkowski A, Tzaribachev N, et al. Efficacy and safety of anakinra therapy in pediatric and adult patients with the autoinflammatory muckle-wells syndrome. Arthritis rheumatism (2011) 63(3):840–9. doi: 10.1002/art.30149

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Dinarello CA, van der Meer JWM. Treating inflammation by blocking interleukin-1 in humans. Semin Immunol (2013) 25(6):469–84. doi: 10.1016/j.smim.2013.10.008

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Jesus AA, Goldbach-Mansky R. IL-1 blockade in autoinflammatory syndromes. Annu Rev Med (2014) 65:223–44. doi: 10.1146/annurev-med-061512-150641

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Cantarini L, Lopalco G, Cattalini M, Vitale A, Galeazzi M, Rigante D, et al. Interleukin-1: Ariadne’s thread in autoinflammatory and autoimmune disorders. Israel Med Assoc journal: IMAJ (2015) 17(2):93–7.

Google Scholar

14. Schett G, Dayer J-M, Manger B. Interleukin-1 function and role in rheumatic disease. Nat Rev Rheumatol (2016) 12(1):14–24. doi: 10.1038/nrrheum.2016.166

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Szekanecz Z, McInnes IB, Schett G, Szamosi S, Benkő S, Szűcs G. Autoinflammation and autoimmunity across rheumatic and musculoskeletal diseases. Nat Rev Rheumatol Springer Sci Business Media LLC (2021) 17(10):585–95. doi: 10.1038/s41584-021-00652-9

CrossRef Full Text | Google Scholar

16. Nedjai B, Hitman GA, Yousaf N, Chernajovsky Y, Stjernberg-Salmela S, Pettersson T, et al. Abnormal tumor necrosis factor receptor I cell surface expression and NF-kappaB activation in tumor necrosis factor receptor-associated periodic syndrome. Arthritis rheumatism (2008) 58(1):273–83. doi: 10.1002/art.23123

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Bachetti T, Chiesa S, Castagnola P, Bani D, Di Zanni E, Omenetti A, et al. Autophagy contributes to inflammation in patients with TNFR-associated periodic syndrome (TRAPS). Ann rheumatic diseases BMJ (2013) 72(6):1044–52. doi: 10.1136/annrheumdis-2012-201952

CrossRef Full Text | Google Scholar

18. Hua Y, Shen M, McDonald C, Yao Q. Autophagy dysfunction in autoinflammatory diseases. J Autoimmun (2018) 88:11–20. doi: 10.1016/j.jaut.2017.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Wu N, Wu D, Zhao M, Miao J, Yu W, Wang Y, et al. Clinical benefits of TNF-α inhibitors in Chinese adult patients with NLRP3-associated autoinflammatory disease. J Internal Med (2021) 290(4):878–85. doi: 10.1111/joim.13334

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Tang X, Huang Y, Lei J, Luo H, Zhu X. The single-cell sequencing: new developments and medical applications. Cell biosci (2019) 9(1):53. doi: 10.1186/s13578-019-0314-y

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with harmony. Nat Methods (2019) 16(12):1289–96. doi: 10.1038/s41592-019-0619-0

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Aibar S, González-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods (2017) 14(11):1083–6. doi: 10.1038/nmeth.4463

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zenz R, Eferl R, Scheinecker C, Redlich K, Smolen J, Schonthaler HB, et al. Activator protein 1 (Fos/Jun) functions in inflammatory bone and skin disease. Arthritis Res Ther (2008) 10(1):201. doi: 10.1186/ar2338

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Singh G. C-reactive protein and erythrocyte sedimentation rate: Continuing role for erythrocyte sedimentation rate. Adv Biol Chem (2014) 04(01):5–9. doi: 10.4236/abc.2014.41002

CrossRef Full Text | Google Scholar

25. Kany S, Vollrath JT, Relja B. Cytokines in inflammatory disease. Int J Mol Sci (2019) 20(23):6008. doi: 10.3390/ijms20236008

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Hua Y, Wu D, Shen M, Yu K, Zhang W, Zeng X, et al. Phenotypes and genotypes of Chinese adult patients with systemic autoinflammatory diseases. Semin Arthritis rheumatism (2019) 49(3):446–52. doi: 10.1016/j.semarthrit.2019.05.002

CrossRef Full Text | Google Scholar

27. Diamond MS, Farzan M. The broad-spectrum antiviral functions of IFIT and IFITM proteins,” Nature reviews. Immunology (2013) 13(1):46–57. doi: 10.1038/nri3344

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Gillespie M, Jassal B, Stephan R, Milacic M, Rothfels K, Senff-Ribeiro A, et al. The reactome pathway knowledgebase 2022. Nucleic Acids Res (2022) 50(D1):D687–92. doi: 10.1093/nar/gkab1028

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Hartung HP, Hadding U. Synthesis of complement by macrophages and modulation of their functions through complement activation. Springer Semin immunopathol (1983) 6(4):283–326. doi: 10.1007/bf02116277

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Davies LC, Jenkins SJ, Allen JE, Taylor PR. Tissue-resident macrophages. Nat Immunol (2013) 14(10):986–95. doi: 10.1038/ni.2705

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Italiani P, Boraschi D. From monocytes to M1/M2 macrophages: Phenotypical vs. functional differentiation. Front Immunol (2014) 5:514. doi: 10.3389/fimmu.2014.00514

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Kapellos TS, Bonaguro L, Gemünd I, Reusch N, Saglam A, Hinkley ER, et al. Human monocyte subsets and phenotypes in major chronic inflammatory diseases. Front Immunol (2019) 10:2035. doi: 10.3389/fimmu.2019.02035

PubMed Abstract | CrossRef Full Text | Google Scholar

33. La Manno G, Lange M, Peidli S, Wolf FA, Theis FJ. RNA Velocity of single cells. Nature (2018) 560(7719):494–8. doi: 10.1038/s41586-018-0414-6

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Setty M, Kiseliovas V, Levine J, Gayoso A, Mazutis L, Pe’er D. Characterization of cell fate probabilities in single-cell data with palantir. Nat Biotechnol (2019) 37(4):451–60. doi: 10.1038/s41587-019-0068-4

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Bergen V, Lange M, Peidli S, Wolf FA, Theis FJ. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat Biotechnol (2020) 38(12):1408–14. doi: 10.1038/s41587-020-0591-3

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Witsell AL, Schook LB. Tumor necrosis factor alpha is an autocrine growth regulator during macrophage differentiation. Proc Natl Acad Sci United States America (1992) 89(10):4754–8. doi: 10.1073/pnas.89.10.4754

CrossRef Full Text | Google Scholar

37. Gattorno M, Pelagatti MA, Meini A, Obici L, Barcellona R, Federici S, et al. Persistent efficacy of anakinra in patients with tumor necrosis factor receptor-associated periodic syndrome. Arthritis rheumatism (2008) 58(5):1516–20. doi: 10.1002/art.23475

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Holzinger D, Kessel C, Omenetti A, Gattorno M. From bench to bedside and back again: translational research in autoinflammation. Nat Rev Rheumatol (2015) 11(10):573–85. doi: 10.1038/nrrheum.2015.79

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM 3rd, et al. Comprehensive integration of single-cell data. Cell (2019) 177(7):1888–1902.e21. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Germain P-L, Lun A, Garcia Meixide C, Macnair W, Robinson MD. Doublet identification in single-cell sequencing data using scDblFinder. F1000Research (2021) 10:979. doi: 10.12688/f1000research.73600.2

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res (2003) 13(11):2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Chen T, Chen X, Zhang S, Zhu J, Tang B, Wang A, et al. The genome sequence archive family: Toward explosive data growth and diverse data types. Genomics Proteomics Bioinf (2021) 19(4):578–83. doi: 10.1016/j.gpb.2021.08.001

CrossRef Full Text | Google Scholar

43. CNCB-NGDC Members and Partners. Database resources of the national genomics data center, China national center for bioinformation in 2022. Nucleic Acids Res (2022) 50(D1):D27–38. doi: 10.1093/nar/gkab951

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Ouyang JF, Kamaraj US, Cao EY, Rackham OJL. ShinyCell: Simple and sharable visualisation of single-cell gene expression data. Bioinformatics (2021). doi: 10.1093/bioinformatics/btab209

PubMed Abstract | CrossRef Full Text | Google Scholar