Article Text
PDF
Abstract
Objective To combine targeted transcriptomic and proteomic data in an unsupervised hierarchical clustering method to stratify patients with childhood-onset SLE (cSLE) into similar biological phenotypes, and study the immunological cellular landscape that characterises the clusters.
Methods Targeted whole blood gene expression and serum cytokines were determined in patients with cSLE, preselected on disease activity state (at diagnosis, Low Lupus Disease Activity State (LLDAS), flare). Unsupervised hierarchical clustering, agnostic to disease characteristics, was used to identify clusters with distinct biological phenotypes. Disease activity was scored by clinical SELENA-SLEDAI (Safety of Estrogens in Systemic Lupus Erythematosus National Assessment-Systemic Lupus Erythematosus Disease Activity Index). High-dimensional 40-colour flow cytometry was used to identify immune cell subsets.
Results Three unique clusters were identified, each characterised by a set of differentially expressed genes and cytokines, and by disease activity state: cluster 1 contained primarily patients in LLDAS, cluster 2 contained mainly treatment-naïve patients at diagnosis and cluster 3 contained a mixed group of patients, namely in LLDAS, at diagnosis and disease flare. The biological phenotypes did not reflect previous organ system involvement and over time, patients could move from one cluster to another. Healthy controls clustered together in cluster 1. Specific immune cell subsets, including CD11c+ B cells, conventional dendritic cells, plasmablasts and early effector CD4+ T cells, differed between the clusters.
Conclusion Using a targeted multiomic approach, we clustered patients into distinct biological phenotypes that are related to disease activity state but not to organ system involvement. This supports a new concept where choice of treatment and tapering strategies are not solely based on clinical phenotype but includes measuring novel biological parameters.
- Disease Activity
- Lupus Erythematosus, Systemic
- Autoimmune Diseases
- Cytokines
Data availability statement
Data are available upon reasonable request. De-identified data including R-scripts underlying this article are available to all scientists upon request.
This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See: http://creativecommons.org/licenses/by-nc/4.0/.
Statistics from Altmetric.com
WHAT IS ALREADY KNOWN ON THIS TOPIC
The quest for personalised therapies, which are more efficacious while avoiding unnecessary side-effects, has fuelled omic studies investigating biological pathways in primarily adults with SLE to find novel biomarkers for disease activity states. Children could benefit the most from better disease understanding as they have a more severe disease than adults with SLE.
WHAT THIS STUDY ADDS
Targeted multiomics, combined with unsupervised hierarchical clustering analysis performed in children with SLE, led to the identification of clusters of patients with distinct biological phenotypes that were associated with disease activity states.
HOW THIS STUDY MIGHT AFFECT RESEARCH, PRACTICE OR POLICY
Clinical decision-making in SLE and patient outcome will be improved by better understanding of the biological pathways underlying active disease and disease remission. Multiomic studies have shown their value in adult-onset SLE where the first steps are made towards stratification of patients to guide treatment choices. With the results of our study, we made a start to achieve the same for children with SLE.
Introduction
Childhood-onset SLE (cSLE) is a devastating relapsing remitting autoimmune disease characterised by significant heterogeneity between patients.1 Before the age of 30 years, 5%–10% of the patients with cSLE will have died as a result of severe disease.2–4 At 30 years of age, the majority of the patients will have developed irreversible damage, due to the disease or its treatment. Prednisone accounts for around 50% of the damage, which is a reason to limit or withdraw prednisone completely.2 3 Response to therapy and when to taper or withdraw medication cannot be clearly predicted by clinical disease phenotypes. Therefore, personalised treatment strategies are needed, especially those that not only use clinical characteristics but also novel biomarkers reflecting underlying immunological pathway activation.
The quest for personalised therapies that are more efficacious while avoiding unnecessary side-effects has fuelled studies investigating biological pathways in (c)SLE to stratify patients in distinct groups. These studies mostly used single-layer biological techniques (primarily transcriptomics, next to genetics and flow cytometry), which could differentiate patients from healthy controls (HCs) and were moderately successful in stratifying patients in clinically distinct groups.5–7 However, the overlap between patient groups was considerable and the clinical usefulness of the identified biological parameters remains to be determined.
The clinical complexity of cSLE may be impossible to dissect by using single-layer biological pathway analysis. Combining data from multiple biological layers to stratify patients with (c)SLE in groups with similar biological pathway activation has been done and is a promising new option to improve patient stratification.8–11 Type I interferons (IFNS) play a pivotal role in the pathogenesis of SLE and blocking their activity results in improvement of disease activity.12 13 Additionally, other pro-inflammatory and anti-inflammatory cytokines, like interleukin (IL)-6, IL-10, IL-8, IL-17 and chemokine (C-C motif) ligand 2 (CCL2), have been associated with SLE disease activity and proposed as potential biomarkers.14–16 Moreover, biological agents targeting cytokines and their receptors have shown their value for the treatment of SLE, thus underlining the contribution of these proteins to disease pathogenesis.13 17 Hence, combining transcriptomic with proteomic (eg, cytokine) data would be a promising option. Unsupervised computational models integrate high-dimensional data from multiple ‘-omes’ (eg, transcriptome and proteome) and are a method to find combinations of biological markers that stratify patients into homogeneous subgroups.
Here we combined targeted transcriptomic and proteomic data in an unsupervised hierarchical clustering method, agnostic to disease characteristics, to stratify patients into similar biological phenotypes. Moreover, we identified unique immune cell subsets that characterise these clusters. Lastly, we studied the relation of these biological phenotypes to the patient’s clinical characteristics.
Methods
Study population
In this exploratory study, two cohorts of patients with cSLE were studied. These patients were preselected in three equally sized groups based on disease activity state (treatment naïve at diagnosis, Low Lupus Disease Activity State (LLDAS), disease flare) (table 1). The discovery cohort consisted of 17 patients with cSLE. Blood samples from 12 of 17 patients with cSLE were collected at a single time point. For 5 of 17 patients, two to three blood samples per patient were used; these were taken at different time points when the patient was in a different disease activity state. The replication cohort consisted of 12 other patients with cSLE with blood samples taken at a single time point (online supplemental table 1). These patients were either in LLDAS or had active disease (clinical SELENA-SLEDAI (Safety of Estrogens in Systemic Lupus Erythematosus National Assessment-Systemic Lupus Erythematosus Disease Activity Index) ≥4). All patients fulfilled the Systemic Lupus International Collaborating Clinics classification criteria.1 Additionally, eight age-matched HCs, without symptoms of underlying viral infections or the use of any medications, were included (online supplemental table 1).
Supplemental material
Patient characteristics in the discovery cohort
Gene signatures
A selection of genes from four previously described gene modules was determined by reverse transcription-PCR.5 These gene modules have been associated with disease features and recently been translated into applicable gene signatures represented by the selected genes.18 The genes are listed in table 2.
Overview of genes and cytokines
Cytokine measurement
Fifteen cytokines were measured using an ELLA Simple Plex system (Protein simple, San Jose, California, USA) and four subsequent cytokines were measured using a conventional ELISA according to the manufacturer’s protocol (table 2). For further analysis, we used cytokines with a value above the lower limit of quantification in at least 20% of the patients to prevent skewing of the data. Twelve of the 19 cytokines that were measured fulfilled this requirement.
High-dimensional flow cytometry
Peripheral blood mononuclear cells from patients and HCs were stained with a 40-colour antibody panel as previously described19 and analysed using a five-laser Aurora spectral flow cytometer (Cytek Biosciences, California, USA). The unsupervised and statistical inference portions of the flow cytometry analysis were performed using the OMIQ data analysis software (www.omiq.ai). The unsupervised analysis method based on surface markers without any two-dimensional gating was used.20 The workflow included running flowAI to check for changes in channels over acquisition time, uniform manifold approximation and projection for dimensionality reduction, flowSOM for clustering and edgeR for statistical inference. For the statistical comparisons of abundance, the Flow Cytometry Standard files were subsampled to ensure the same number of events was included per group (either per biological phenotype or disease activity state).
Statistical analysis
For the unsupervised hierarchical clustering and principal component analysis (PCA), we loaded the patient’s clinical, gene expression and cytokine data into R (V.4.0.4). Only cytokines which were positive in >20% of the patients were used in the analysis. If a cytokine was included in the analysis and a sample had a value below detection limit, this value was set as the lower limit of quantification. Standardised scores (z-scores) were calculated based on the log10-transformed gene expression and cytokine levels. Unsupervised hierarchical clustering was performed using Ward’s Hierarchical Agglomerative Clustering Method (ward.d2) on the Euclidean distances of the z-scores. The optimal number of clusters for both cohorts was assigned with the NbClust (V.1.0.12) package in R. Subsequently, heatmaps were plotted using the R package Complexheatmap (V.2.10.0).21 PCA was performed and visualised based on the z-scores using the R function prcomp. The grouping of patients in the different clusters was agnostic to parameters reflecting clinical disease phenotype.
The Kruskal-Wallis test (three groups) was used to analyse comparisons between medians. Levene’s test was used for homogeneity of variance across the groups. Dunn’s multiple comparisons post hoc test was used to correct for multiple testing as indicated in the legend of each figure. Spearman’s correlation coefficient was calculated to assess correlation. Values of p <0.05 were considered statistically significant. Non-significant results are not shown in figures. The false discovery rate (FDR) was used to compare immune cell subsets between the biological phenotypes. A cut-off value of FDR <0.1 was used to indicate significant differences between the populations. Graphpad Prism V.8.0 (Graphpad Software, La Jolla, California, USA) and R statistical software were used for graph design and statistical analysis.
Results
Unsupervised hierarchical clustering identifies three unique biological phenotypes associated with disease activity state
Gene expression and cytokine levels were measured in 23 samples from 17 patients with cSLE. These data were combined and subjected to an unsupervised hierarchical clustering approach agnostic to any clinical or laboratory disease characteristics. This approach identified three clusters, further termed as biological phenotypes (figure 1A). Patients within each cluster had a unique and remarkably similar combination of differentially expressed genes and cytokine profiles. Patients in cluster 1 had a low gene expression and cytokine profile, and were termed low biological phenotype (LBP). Patients in cluster 2 had a high gene expression and high cytokine profile, and were termed high biological phenotype (HBP). Lastly, patients in cluster 3 had a high gene expression and a low cytokine profile and were termed mixed biological phenotype (MBP) (figure 1A and online supplemental figure 1).
Unsupervised hierarchical clustering identifies three unique biological phenotypes associated with disease activity state. (A) Unsupervised hierarchical clustering using Ward’s agglomerative method and passing the Euclidean distance between samples and using row-based log-transformed z-scores identified three clusters. Visit numbers per patient are depicted at the bottom. (B) Jaccard index of the three identified clusters after k-means bootstrapping (measures how many times cells consistently clustered together). Orange dashed bar indicates the best range for cluster stability. (C) Principal component analysis showing that the three clusters lie apart from each other. The first two components with their percentage of variance are shown in parentheses. (D) Correlation plot depicting the inter/intracorrelation between genes and cytokines (N=23). Blue indicates a positive and red indicates a negative correlation. Numbers indicate Spearman’s r, squares indicated a significant correlation (p<0.05). Red-blue colour indicates the z-scores. HBP, high biological phenotype; LBP, low biological phenotype; LLDAS, Low Lupus Disease Activity State; MBP, mixed biological phenotype; SELENA-SLEDAI, Safety of Estrogens in Systemic Lupus Erythematosus National Assessment-Systemic Lupus Erythematosus Disease Activity Index.
To confirm that our clustering approach was not skewed by the patients from whom we had multiple samples collected at different time points, we repeated the unsupervised hierarchical clustering with only data from the sample collected at the first time point of each patient (online supplemental figure 2). The clusters did not change, confirming that the biological phenotypes were linked to disease activity.
Next, to assess whether these biological phenotypes were patient or disease activity state dependent, we studied the five patients from whom multiple time points in different disease activity states were available. In three out of five patients, a change in the biological phenotype was directly associated with a change in disease activity state. The biological phenotypes from the other two patients did not change, but stayed in the MBP group, while the disease activity state changed from flare to LLDAS or vice versa. Notably, these two patients were the only patients in LLDAS in the MBP group, and of all patients in LLDAS, they were the only two to have a clinical SELENA-SLEDAI above 0 (figure 1A).
Cluster stability in the study cohort was assessed by resampling the dataset and computing the Jaccard index to the original cluster.22 A stable cluster should yield a Jaccard index of >0.75.22 After 1000 iterations, the LBP, HBP and MBP group showed a Jaccard index of 0.86, 0.84 and 0.91, respectively (figure 1B), indicating the stability of the three clusters. Next, a PCA showed a clear difference between the identified clusters (figure 1C). PC1 and PC2 combined indicated a proportion of variation of 62.64%. These findings underscored the robustness of the identified clusters.
Subsequently, we studied the association between the gene expression and cytokine data to explore the driving mechanisms underlying the identified biological phenotypes. While a strong correlation was observed between the tested genes as well as between the tested cytokines, a much lower correlation was present between these two datasets (figure 1D). Additionally, while the gene expression and cytokine data are discordant in the MBP group, gene expression and cytokine data in the other two groups are concordant (online supplemental figure 1). This supports the additive value of combined analysis of transcriptional and cytokine data in a single analysis and indicates that genes and cytokines have independent roles in the network that underlies the biological phenotypes.
Validation of biological phenotypes in an independent cSLE cohort
In order to validate the biological phenotypes, we applied our clustering approach on a second cohort of patients with cSLE (n=12), who were either in LLDAS or had active disease (clinical SLEDAI >4). From these patients, we had transcriptional data of all targeted genes and serum analysis of 7 of the 12 cytokines studied in the discovery cohort. Therefore, we adjusted the clustering algorithm in our discovery cohort regarding the number of cytokines. All patients had the same biological phenotype as they were initially assigned to (online supplemental figure 3A). To find out whether a specific cytokine may be dominant in the cluster analysis, we reanalysed the data by removing one cytokine at a time to see the effect on the clusters. In this analysis, any cytokine could be removed except for IL-8 and CCL2, which seem to be the main drivers of these clusters (data not shown). Applying our clustering method to the replication cohort showed a clear distinction between the patients in LLDAS and with active disease (online supplemental figure 3B). PCA, representing the gene expression and cytokine profiles, indicated that measurements in the replication cohort showed good agreement with those in the discovery cohort (online supplemental figure 3C). Together, these data confirmed the validity and robustness of the study outcome, showing identical clustering despite the relative low patient numbers in each cohort.
Unsupervised hierarchical clustering groups HCs within the low biological phenotype cluster of patients with cSLE
When data from HCs were added to our clustering approach, patients with cSLE in the discovery cohort remained in the same cluster (online supplemental figure 4A). Remarkably, all HCs grouped together with patients with cSLE in the LBP cluster and PCA showed a clear overlap between these two groups (online supplemental figure 4A,B). PCA of the combination of HCs and patients in the LBP group from the discovery and replication cohort also showed overlap (online supplemental figure 4C). Patients with cSLE, who clustered together with HCs, were primarily patients in LLDAS who were on medication, with the majority treated with both hydroxychloroquine (HCQ) and mycophenolate mofetil (MMF). Only one of these patients with cSLE was on low-dose prednisone (0.05 mg/kg/day) (table 1). Interestingly, all patients in LBP with LLDAS without prednisone remained in LLDAS for at least 1 year after stopping prednisone. These data illustrate that the biological phenotype of HCs with the used biological markers was not different from patients with cSLE with low disease activity in this explorative study.
Biological phenotypes reflect disease activity states
Most patients in the LBP group were in LLDAS (six of eight). The other two patients were included at diagnosis, with only involvement of the mucocutaneous domain next to serological activity (SELENA-SLEDAI 2 and 6). The median clinical SELENA-SLEDAI in this cluster was 0 (0–6) (figure 2A). The HBP group contained patients at diagnosis and one patient with a flare. The median clinical SELENA-SLEDAI in this cluster was 10 (4–15; p=0.011 (Kruskal-Wallis test with Dunn’s multiple comparisons post hoc)) (figure 2A). Lastly, for patients in the MBP group, which concerns a mix of patients in LLDAS, disease flare or at diagnosis, the median clinical SELENA-SLEDAI was 5.5 (0–14; p=0.026 (Kruskal-Wallis test with Dunn’s multiple comparisons post hoc)) (figure 2A). When disease duration was compared between the three clusters, patients in the HBP group had median disease duration of 0 (0–6.0) years, while patients in the MBP and LBP cluster had median disease duration of 1.6 (0–3.7) and 3.2 (0–4.1) years, respectively (figure 2B). This could suggest that disease duration could be the driving factor for the biological phenotypes. However, this is contradicted by the finding that two patients who were at diagnosis were also found in the LBP group and patients with a disease flare were in the HBP group. When the days in LLDAS were compared between LLDAS patients in the LBP and MBP group, we found that LLDAS patients in the LBP group were already in LLDAS for a much longer period of time (median 715 days) when compared with the two LLDAS patients in the MBP group (median time in LLDAS 198 days) (figure 2B). Together, these data underline that not only disease duration but particularly (biological) disease activity might be the driving factor for the biological disease phenotypes.
Biological phenotypes are associated with disease states and not influenced by the use of medication. (A) Clinical SELENA-SLEDAI score of patients with different biological phenotypes. Each dot represents individual patients. (B) Heatmap indicating the disease activity state, clinical SELENA-SLEDAI, prednisone usage in mg/kg/day, disease duration in years, days in LLDAS, hydroxychloroquine (HCQ) and mycophenolate mofetil (MMF) use. (C) Heatmap indicating organ domain involvement ever. Bar plot indicates the number of organ domains involved per patient in each cluster. (D) Bar plots indicating 10 laboratory parameters in patients with different biological immune phenotypes. Kruskal-Wallis test with Dunn’s multiple comparisons post hoc test was used for relation between the clusters. CRP, C reactive protein; ESR, erythrocyte sedimentation rate; HBP, high biological phenotype; LBP, low biological phenotype; LLDAS, Low Lupus Disease Activity State; MBP, mixed biological phenotype; SELENA-SLEDAI, Safety of Estrogens in Systemic Lupus Erythematosus National Assessment-Systemic Lupus Erythematosus Disease Activity Index.
Medication use does not explain the biological phenotypes identified by hierarchical clustering
When comparing medication use of the patients in the different clusters, HCQ and MMF usage were similar in the LBP and MBP groups, yet patients clustered separately (figure 2B). Notably, four out of five patients in the HBP group were treatment naïve at blood sampling (figure 2B). Interestingly, almost all (8 of 10) patients in the MBP group were treated with prednisone, where one out of eight and zero out of five were treated with prednisone in the LBP and HBP groups, respectively (figure 2B). In the MBP group, no correlation was found between prednisone dosage and gene expressions or cytokine levels (online supplemental figure 5 and table 1) suggesting that medication use does not explain the biological phenotypes identified by our clustering approach.
Biological phenotypes are not driven by organ system involvement or currently used laboratory parameters for SLE disease activity
To assess whether biological phenotypes reflect specific organ system involvement, we studied the organ systems that were ever affected in each patient. No difference was observed between the number of organ systems ever involved between the three clusters (figure 2C). Studying the autoantibody profiles of the patients in the clusters, no difference between patients with cSLE in the different clusters was found (table 1). Lastly, we studied whether routine laboratory parameters to measure SLE disease activity were different between the three biological phenotypes (figure 2D). Not surprisingly, anti-dsDNA, C reactive protein, erythrocyte sedimentation rate and IgG levels were highest in the HBP group, the group with newly diagnosed patients (figure 2D). However, C3, C4, thrombocyte, lymphocyte and neutrophil counts did not differ between any of the groups and importantly could not differentiate between patients in LBP and MBP. These data indicate that the identified biological phenotypes are not primarily driven by organ system involvement, and currently used laboratory parameters to measure cSLE disease activity are of limited influence on biological disease phenotype (LBP, MBP) during follow-up of the patient.
High-dimensional flow cytometry identifies immune cell subsets that differ between biological phenotypes
To investigate the cellular landscape that characterised the three identified biological phenotypes, we applied 40-colour spectral flow cytometry (figure 3A). FlowSOM analysis identified 32 distinct immune cell subsets (online supplemental figure 6). In line with previous findings,6 there was a clear difference in the immune cell subsets from HCs compared with the patients. In total, 14 of 32 populations differed significantly between these two groups, with central memory CD8+ T cells being the cell population that showed the highest fold difference (figure 3B,C and online supplemental table 2).
High-dimensional flow cytometry identifies immune subsets that differ between biological phenotypes. (A) Population cluster identification in high-dimensional 40-colour flow cytometry data using uniform manifold approximation and projection (UMAP) dimensionality reduction. Visualisation is performed on combined patients with childhood-onset (cSLE) and healthy control (HCs) data (n=31). (B) Volcano plot comparing patients/HCs. Positive fold changes indicate cellular populations increased in patients, while negative fold changes indicate decreases in patients. (C) UMAP plots of HCs (n = 8) and all patients with cSLE (n=23). Populations that are increased in each group are indicated by a dashed line. (D) Volcano plot comparing HBP (N=5) versus LBP (N=8). Positive fold changes indicate cellular populations increased in HBP compared with LBP. (E) Volcano plot comparing HBP (N=5) versus MBP (N=10). Positive fold changes indicate cellular populations increased in HBP compared with MBP. (F) Volcano plot comparing LBP (N=8) versus MBP (N=10). Positive fold changes indicate cellular populations increased in MBP compared with LBP. (G) Bar plots depict the percentage of cells in total population. The FDR is depicted as statistical analysis. DCs, dendritic cells; FDR, false discovery rate; HBP, high biological phenotype; LBP, low biological phenotype; MBP, mixed biological phenotype; NK, natural killer; pDCs, plasmacytoid DCs.
When we compared the difference between the biological phenotypes in patients with cSLE, the HBP group clearly differed from the other two biological phenotypes: 16 of 32 of the identified immune cell subsets differed between the HBP and LBP groups (figure 3D and online supplemental table 3), and 10 of 32 populations differed between the HBP and MBP groups (figure 3E,F and online supplemental table 4). When comparing the LBP and MBP groups, only 2 of 32 populations differed (figure 3F and online supplemental table 5). Of interest, CD11c+ CD16− dendritic cells (DCs) were increased in the LBP group, while a unique CD11c+ B cell population was increased in the MBP and HBP group (figure 3F). To assess associations with the biological phenotypes, we studied the correlation between the 32 immune cell subsets and the clinical SELENA-SLEDAI. Nine immune cell populations were significantly associated with the clinical SELENA-SLEDAI (online supplemental figure 7). Of these, IgG+ memory B cells and the central memory CD8+ T cells did not show any significant difference between the biological phenotypes (online supplemental figure 7). The seven other populations significantly differed between the identified biological phenotypes (figure 3G). Here, CD11c+ B cells, plasmablasts and early effector CD4+ T cells were increased in the HBP group compared with the LBP group. On the other hand, conventional DCs, TCRγδ+ CCR7− CD45RA−, plasmacytoid DCs and mature natural killer CD159+ cells were decreased in the HBP group compared with the LBP group (figure 3G).
Taken together, multiple immune cell subsets differed significantly between the biological phenotypes and were associated with disease activity. This cellular analysis underscores the connection between immune cell subsets and targeted transcriptomic/proteomic data, indicating the additional relevance of these cells in defining biological phenotypes.
Discussion
Our study illustrates that combining targeted transcriptomic and proteomic data is a reliable method to cluster patients in similar clinical and biological phenotypes that represent a specific disease activity state rather than organ system involvement. Moreover, high-dimensional flow cytometry revealed immune cell subsets that were characteristic for these phenotypes. Thus, these biological phenotypes can facilitate patient stratification for optimal treatment choices and provide information to improve interventional clinical trial design.
Transcriptomic studies in SLE have been able to discriminate between patients with low and high disease activity. Yet, gene signatures are often expressed in binary values and remain relatively stable over time, and therefore are not sensitive to changes that occur during clinical follow-up of patients.23 Multiomics have led to ground-breaking findings in the field of oncology indicating their usefulness for other fields like autoimmunity.24 One study, in patients with adult-onset SLE, combined different gene modules and various plasma soluble mediators and identified seven unique molecular profiles.11 However, these profiles did not distinguish between patients with different disease activity scores. Furthermore, this study included multiple gene modules from which the majority had no relation with SLE disease activity.5 In contrast, we chose a targeted approach and selected genes and cytokines that previously were associated with disease activity5 18 expecting a higher clinical relevance from these genes and cytokines. Indeed, this approach enabled a more precise stratification of patients in groups with similar biological phenotypes that are related to disease activity state.
For the first time, we were able to link almost all patients with cSLE in LLDAS to a specific biological disease phenotype, namely the LBP group. Even more, HCs clustered with these LLDAS patients in the LBP group. This shows that patients on medication can achieve comparable combinations of gene expression and cytokine profiles as HCs and additionally increases our understanding of LLDAS being associated with reduced risk of flares, less damage development and corticosteroid sparing.25 26 Hence, we speculate that for more accurate characterising of LLDAS, analysis of these novel biological parameters could be used. This concept is in line with findings in patients with ulcerative colitis in which histological remission at the gut level, in addition to clinical remission, is associated with corticoid sparing and less hospitalisation for disease flares.27
Medication usage can be an important confounding factor in analysing biological profiles in patients. For example, high-dose corticosteroids have been reported to have an inhibitory effect on multiple proinflammatory cytokines and the type I IFN signature.28 29 In our study, patients who used prednisone (primarily in the MBP group) still had high levels of type I IFN-induced genes and prednisone dose did not correlate with gene expression or protein levels. Moreover, between the different biological phenotypes, no major differences in HCQ and MMF usage were observed. Remarkably, almost all patients in the LBP group, who have low gene and cytokine profiles, were in LLDAS and used HCQ and disease-modifying antirheumatic drugs but not prednisone. Interestingly, these patients remained in LLDAS for at least 1 year after stopping prednisone. This underscores that continued prednisone use is not needed to remain in LLDAS.
As expected, pro-inflammatory cytokines are upregulated in the HBP group compared with the LBP group. An interesting finding is the upregulation of transforming growth factor beta (TGF-β) in HCs, patients in the LBP and also in patients in the HBP group. This latter finding is somewhat contradictory to literature showing that low concentrations of TGF-β are associated with high disease activity.30 However, we also know that TGF-β can have dual roles with the activation of resting monocytes but also the inhibition of activated macrophages.31 32 This in line with the theory that often the effects of TGF-β are context dependent and more research is needed to understand the biological role of TGF-β.33 Additionally, when analysing whether one specific cytokine was dominant in driving the clusters, we found that IL-8 and CCL2 were the main drivers of the clusters in our explorative study. Although these findings are suggesting that the analyte number could be reduced, additional studies with larger datasets should first confirm this.
In SLE, gene signatures and cytokines have been associated with involvement of certain organ systems.23 34 35 Currently, organ system involvement is seen as an important trait that guides treatment strategies. Here we show that biological phenotypes are able to change and are more a reflection of the patients’ disease activity state than that of the specific organ system(s) involved. This supports a new concept where choice of treatment and tapering strategies are not solely based on clinical phenotype but includes measuring these novel biological parameters. This concept is supported by our recent work showing that patients with cSLE, irrespective of organ system involvement, reached LLDAS in a median time of 6 months. This was achieved by a treatment strategy focusing on limited use of corticosteroids and early introduction of immunosuppressives, specifically when corticosteroids could not be tapered within 3 months.36
High-dimensional flow cytometry confirmed previous findings that immune cell subsets differ between HCs and patients with SLE.10 A caveat in the previous study in cSLE that focused on immune cell subset differences is the high number of patients with low disease activity (SLEDAI ≤4).6 In our study, we preselected patients in equally sized groups with different disease activity states to overcome this caveat. Interestingly, distinct immune cell populations differed between the patients with different biological phenotypes. This indeed shows the interconnection between the biological phenotypes based on a combination of transcriptional and cytokine data and immune cells in peripheral blood, which hints towards the idea that these immune cells could also be used to reflect the identified biological phenotypes. In the ideal setting, measuring a small number of cell populations could predict a specific clinical disease course and guide treatment strategy, as has been implemented in oncology.37 Possible candidates for this purpose are the CD11c+ B cell and CD11c+ CD16- DC populations identified in our study, as these cells were able to discriminate between patients in the LBP and MBP group, which may be a difficult task based on clinical characteristics. The CD11c+ B cell population has previously been associated with disease activity and proposed to be driven by an extrafollicular B cell activation route in SLE.38 39 Likewise, decreased blood levels of the CD11c+ CD16- DC population have been associated with active disease in SLE.40
Better understanding of biological pathways underlying active disease, LLDAS and remission will facilitate tailored treatment strategies and improve patient outcome. Children would benefit the most as they have a more severe disease than adults with SLE.3 Multiomic data have shown their value in adult-onset SLE where the first steps are made towards stratification of patients to guide treatment choices.41 Our explorative study in children with SLE for the first time links disease activity states to biological phenotypes. We identified three clusters, in which the MBP group is the most interesting group, as from a clinical point of view, we would not have clustered these patients together. The next step is to extend patient numbers and perform machine learning techniques on combined targeted gene expression, cytokine, flow cytometry and clinical data to facilitate further understanding of underlying biological mechanisms and make a start towards stratification of patients to guide treatment choices in cSLE.
Our study has limitations. First, we have a low number of patients in our analysis. Expanding this study in a larger set of patients with different ethnic backgrounds will be our next step. The validation of the results by using a PCA, bootstrapping and a replication cohort makes us confident that expanding the cohort will produce similar results. Second, some of the analyses performed are based on multiple time points from the same individuals. Although we performed extra analyses that suggest that any underlying correlations from repeated samples from the same individual did not appear to substantively affect the clustering, we recognise that repeated observations from the same person are correlated. Third, we have only used a combination of (targeted) transcriptomics and proteomics, while there are several other omic techniques that may also reflect the biological pathways that play a role in the pathogenesis of SLE.42
In conclusion, here we show that integrating targeted transcriptomic and proteomic data in an unsupervised hierarchical clustering model leads to the identification of biological phenotypes that reliably distinguish patients with cSLE with different disease activity states, which could not be achieved by currently used clinical characteristics. This is a first step towards developing personalised treatment strategies in children with lupus based on clinical phenotype and new biomarkers. Future studies should translate and simplify these tools by reducing the number of analytes that need to be measured, which will facilitate implementation in daily clinical practice.
Data availability statement
Data are available upon reasonable request. De-identified data including R-scripts underlying this article are available to all scientists upon request.
Ethics statements
Patient consent for publication
Ethics approval
This study involves human participants and was approved by the Medical Ethics Review Committee of the Erasmus University Medical Center, Rotterdam, the Netherlands (MEC-2019-0412/MEC-2005-0137). Written informed consent was obtained from all participants in compliance with the Declaration of Helsinki.
Acknowledgments
The authors thank all the patients with cSLE and HCs for taking part in this study.
References
Supplementary materials
Supplementary Data
This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.
Footnotes
MAV and SK are joint senior authors.
Contributors All authors were involved in drafting the article or revising it critically for important intellectual content, and all authors approved the final version to be published. Study conception and design—MJW, YMM, PDK, SK and MAV. Experimental work—MJW, SJvT, YMM, CGVH-M, HdW, AWL, PDK and MAV. Acquisition of clinical data—MJW, AWL, MJG, AM, MV and SK. Analysis and interpretation—MJW, SJvT, YMM, PDK, MAV and SK. Writing (draft manuscript)—MJW, MAV and SK. Revision of the manuscript—MJW, SJvT, YMM, HdW, CGVH-M, AWL, MJG, AM, MV, PDK, MAV and SK. MAV and SK are responsible for the overall content as the guarantors.
Funding This work was made possible by the support of the Sophia Children’s Hospital Fund (B18-04), NVLE (Dutch patient organization for Lupus, APS, Scleroderma and MCTD) (BP12-1-261) and Dutch Arthritis Society (CO-19-001). The research was performed within the framework of the Erasmus Postgraduate School Molecular Medicine.
Competing interests None declared.
Patient and public involvement Patients and/or the public were not involved in the design, or conduct, or reporting, or dissemination plans of this research.
Provenance and peer review Not commissioned; externally peer reviewed.
© Author(s) (or their employer(s)) 2023. Re-use permitted under CC BY-NC. No commercial re-use. See rights and permissions. Published by BMJ.
Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.