Gene Expression Trend Pattern Analysis in Peripheral Blood From Patients With Preclinical Systemic Sclerosis

Abstract

In this study, we identified significant blood gene expression trend patterns associated with specific immune cell types in early-stage differentiation in systemic sclerosis (SSc) progression. SSc is a heterogeneous autoimmune connective tissue disease, which is mainly classified as diffuse cutaneous SSc (dcSSc) and limited cutaneous SSc (lcSSc) depending on the extent of skin fibrosis [1]. Although SSc is a chronic immune disorder like allergic rhinitis (AR), the antigen type in SSc is a self-antigen, in contrast to a foreign antigen in AR. We previously reported blood immune gene signatures associated with immune cell frequency change and clinical symptoms after nasal allergen challenge in patients with AR, which was used to identify immune genes with significantly different trend patterns after immunotherapy [2, 3]. Investigating time series gene expression from peripheral blood is useful for understanding systemic immune response changes reflecting the progression or pathophysiology of the diseases [2, 4].

We hypothesize that trend pattern analysis of time series blood gene expression, namely regrouping gene expression clusters into simplified trend patterns (up, down, and steady) by a given ratio threshold, can provide instinctual interpretable insight to understand the pathophysiological systemic immune response that occurs during the development of the main SSc manifestations, such as fibrosis of the skin and/or inner organs. Specifically, such an analysis may provide a viewpoint to understand key immune cell frequency or activation changes at critical crossroads during the course of the disease, particularly in early immune disease progression.

To test our hypothesis, we used public data, GSE224849, which Bellocchi et al. recently registered to report a global gene expression study of patients with preclinical systemic sclerosis (PreSSc) [5]. PreSSc patients with the early signs of SSc were classified by the LeRoy & Medsger classification criteria—the presence of Raynaud's phenomenon plus SSc-specific autoantibodies and/or SSc-specific nailfold video-capilaroscopic change without any other sign of definite SSc and/or fibrosis [1, 5]. The RNA-sequencing data were generated using blood from two PreSSc patient groups with different disease progression stages at follow-up visits.

The blood samples from 33 PreSSc patients were collected at baseline and follow-up visits (4 years later). Fourteen patients developed lcSSc (Evolving-PreSSc), a severe form of SSc with skin features such as puffy fingers and/or skin fibrosis, whereas 19 patients remained stable (Stable-PreSSc) at follow-up: Evolving-PreSSc was classified at follow-up visit as the patients with definite SSc, herein lcSSc, by reaching the minimum score of 9 based on the 2013 ACR/EULAR classification criteria for SSc. Evolving-PreSSc was characterized as having mainly puffy fingers and/or partly telangiectasia (spider veins) or sclerodactyly: Patients with a definite SSc but without skin fibrosis yet with puffy fingers were also categorized in the lcSSc group [5]. Sixteen healthy controls were characterized who had similar age, ethnicity, and geographical residence to patients. To investigate the gene expression trends across time, we had a set of assumptions:

Based on Assumption 1, we used the gene expression measurements in healthy controls as the measurements at pre-onset (hypothetical time-point) in PreSSc, which allowed us to have three-time points with the same start point in both PreSSc subgroups (13 evolving PreSSc patients—we excluded one patient with missing data—and 19 stable PreSSc patients) (Figure 1A): Pre-onset, baseline post-onset, and follow-up post-onset; two periods, period 1 (from pre-onset to baseline post-onset) and period 2 (from baseline post-onset to follow-up post-onset).

Based on Assumption 2, a ratio threshold—a distinguishable relative change—provides a rule to decide what change is negligible small naturally-occurring fluctuations or oscillation, which was used for selecting genes and regrouping clusters into trends (Figure 1C,F). We used three commonly-applied ratio thresholds: 1.2 (low), 1.5 (medium), and 2.0 (high).

For preprocessing, genes (raw data) with low counts across all samples were filtered out using the edgeR R package (version 4.0.16). Subsequently, in each PreSSc subgroup, the remaining 20 587 genes (TPM normalized data) were evaluated given a ratio threshold: for example, under low-ratio settings, genes with a maximum to minimum ratio of mean expression across time greater or equal to 1.2 were selected (Figure 1C). We applied ANOVA with permutation (lmPerm R package, version 2.1.0; 999 999 permutations) to the filtered genes (a total of 10 580 genes in Stable-PreSSc and 10 421 genes in Evolving-PreSSc) to identify those showing differential expression over time: In Stable-PreSSc there were 5469 significant (False discovery rate, FDR < 0.1) genes at low threshold ratio; 1699 genes at medium-ratio; and 265 genes at high-ratio. In Evolving-PreSSc (with a smaller sample size) there were 17 significant genes at low-ratio; 824 genes at medium-ratio; and 123 genes at high-ratio (Figure 1D).

Hierarchical clustering was performed on the pooled significant genes (5591 genes in Stable-PreSSc; 827 genes in Evolving-PreSSc) in each group using the dtwclust R package (version 6.0.0) after standardizing the mean expression of each gene (mean = 0, SD = 1). Based on the similarity of their expression patterns, genes clustered into 24 groups (Figure 1G). Then, after evaluating the trend pattern of mean values at each time in each cluster, the clusters were regrouped into nine interpretable trends. The nine trend patterns were based on the permutation with repetition of expression trends over the two time periods (n^r = 9, when n = 3: up, down, and steady between two successive time points; r = 2: period 1 and 2) (Figure 1E,G). The ratio thresholds helped to determine whether a trend in a particular time period was negligible (slope 0, steady trend) or significant (upward or downward trend) compared with the other time period (Figure 1F). Most of the overlapped significant genes between Stable- and Evolving-PreSScs showed consistent expression trends in period 1 but varied in period 2 between PreSSc groups, which may imply the association between differential gene expression at period 2 and the differential disease progression between Stable- and Evolving-PreSScs (Figure 2A).

Significant genes (FDR < 0.1) in a given ratio threshold with specific expression trends were applied to cell enrichment analyses using Enrichr (blood cell types, gene set library: HuBMAP_ASCTplusB_augmented_2022) (Figure 2A,B). Stable-PreSSc showed three clusters were significantly associated with immune cells; Evolving-PreSSc, four clusters (Figure 2B). These showed enriched and intriguing results between the two groups besides the down-regulation of cytotoxic/NK cells in Evolving-PreSSc that Bellocchi et al. [5] reported with modular analysis.

In Stable-PreSSc, although CD14-positive monocytes, classical monocytes, were associated with both C1 up-up and C3 up-down patterns, CD14-low-CD16-positive monocytes, nonclassical monocytes, were only associated with C3 up-down pattern (Figure 2B). Monocytes are early responders to pathogens and maintain vascular homeostasis beyond simply being macrophage precursors; particularly nonclassical monocytes that recognize and clear dying endothelial cells to maintain vascular homeostasis can contribute to the inflammation associated with chronic diseases [6]. Intriguingly, peripheral blood nonclassical monocytes were overrepresented in SSc compared with healthy control and characterized by increased prostaglandin E synthase gene expression [7]. Nonclassical monocytes and M2 macrophages have recently been emphasized as key drivers of inflammatory and fibrotic manifestations in SSc [7, 8]. Thus, the normalization of CD14-low-CD16-positive monocytes activated by sensing connective or vascular tissue damage is likely associated with stabilizing the disease as Stable-PreSSc showed. CD4-positive, CD8-positive, and regulatory T cells were significantly associated with the C8 down-steady pattern, still, CD4-positive T cells and neutrophils had the C3 up-down pattern. The monocyte subtype changes and downregulated T cells recalled the significantly different proportion of low-dose aspirin administration in two groups (100% Stable-PreSSc vs. 64% Evolving-PreSSc), the only significant difference in the demographics comparison [5]. Aspirin enhances the anti-inflammatory and anti-regulatory functional activities of monocytes while it downregulates CD16 and CD40 expression [9]. Dose-dependent effects of aspirin on immune cells and endothelial cells are anti-inflammatory but varied [10, 11]. Further studies on prescription timing and dose-dependent effects of aspirin are needed to elucidate the potential benefits in the treatment of SSc (Figure 2C).

In Evolving-PreSSc, B cells were associated with C2 up-steady or C3 up-down patterns: B cells were activated or increased in period 1, and the upregulation remained steady or normalized in period 2. B cell depletion therapies in the early stage of SSc ameliorated skin fibrosis [12]. Reduction of apoptotic signaling in B cell transcriptome in patients with systemic lupus erythematosus compared with healthy controls was associated with puffy fingers [13]. B cell activation and their disease-specific autoantibodies may play a key role in the development of manifestations. In contrast, CD8-positive T cells and Natural Killer (NK) cells were associated with C8 down-steady or C9 down-down patterns. These cell types have important anti-fibrotic effects on the immune system. Although their role may be controversial in fibrosis depending on the target, CD8-positive T cells control fibrosis by inducing apoptosis of myofibroblasts, a key cell type in fibrosis [14]. Depletion of CD8-positive T cells has been shown to exacerbate fibrosis whereas depletion of CD4-positive T cells reduced fibrosis in a mouse model of renal fibrosis [15]. Noticeably, the downward trend was associated with CD16-negative-CD56-bright NK cells exerting immunoregulatory and immunosuppressive effects on both innate and adaptive immune cells as potent cytokine and chemokine producers and cytotoxic cells against autologous activated CD4-positive T cells [16]. The opposing patterns between the upregulation of B cells and the downregulation of CD8-positive T cells and CD16-negative-CD56-bright NK cells may indicate a pro-fibrotic and pro-inflammatory environment exacerbating puffy fingers and skin fibrosis in Evolving-PreSSc (Figure 2B,C).

In both Stable- and Evolving-PreSScs, the gene expression trend associated with CD4-positive-CD45RA-positive effector memory T cells (CD4-positive TEMRA), a CD4-positive memory T cell subtype associated with senescence, were significantly identified but the trend patterns were different: C8 down-steady in Stable-PreSSc but C3 up-down in Evolving-PreSSc (Figure 2B). The upregulated CD4-positive TEMRA in period 1 in Evolving-PreSSc likely contributed to the disease development while considering the relationship between various diseases including autoimmunity and senescent T cells including CD4-positive TEMRA [17].

In addition, although effector CD8-positive αβ T cells had a downward trend in period 1 in both PreSScs, the cause of the trend may be different between PreSScs while considering the other lymphocyte patterns in period 1: downregulation of CD4-positive and CD8-positive T cells and regulatory T cells in Stable-PreSSc; but upregulation of B cells and CD4-positive memory T cells whereas downregulation of CD8-positive T cells and NK cells in Evolving-PreSSc. The downregulation of CD8-positive T cells including effector CD8-positive αβ T cells in Evolving-PreSSc may be associated with type 2 immunity rather than immune inhibition. A relative predominance of type 2 responses (pro-fibrotic response) and of type 1 and 17 responses were associated with, respectively, early and later in the disease course of SSc [18].

The utility of blood gene expression trend pattern analysis was demonstrated by identifying significantly distinguishable alternations of specific immune cell types in early-stage differentiation in SSc progression. Further research is needed to clarify the immune cell-frequency changes or -regulation in peripheral blood while considering T cell senescence, leukocyte migration to target cells, and their relationship with progress in chronic inflammation. More robust conclusions can be drawn from larger cohorts with more longitudinal measures including real pre-onset measurement in the skin/organs and blood, considering immune cell frequency, gene interactions, sophisticated threshold setting, and a single cell-level approach, such as single-cell RNA-sequencing.

Y.W.K. designed the study and performed bioinformatics and statistical analysis; S.J.T. and A.S. supervised the project; Y.W.K. wrote the first draft of the manuscript; and all authors reviewed and contributed to the final version of the manuscript.

The authors declare no conflicts of interest.