Patient ancestry significantly contributes to molecular heterogeneity of systemic lupus erythematosus

Abstract

Gene expression signatures can stratify patients with heterogeneous diseases, such as systemic lupus erythematosus (SLE), yet understanding the contributions of ancestral background to this heterogeneity is not well understood. We hypothesized that ancestry would significantly influence gene expression signatures and measured 34 gene modules in 1566 SLE patients of African ancestry (AA), European ancestry (EA), or Native American ancestry (NAA). Healthy subject ancestry-specific gene expression provided the transcriptomic background upon which the SLE patient signatures were built. Although standard therapy affected every gene signature and significantly increased myeloid cell signatures, logistic regression analysis determined that ancestral background significantly changed 23 of 34 gene signatures. Additionally, the strongest association to gene expression changes was found with autoantibodies, and this also had etiology in ancestry: the AA predisposition to have both RNP and dsDNA autoantibodies compared with EA predisposition to have only anti-dsDNA. A machine learning approach was used to determine a gene signature characteristic to distinguish AA SLE and was most influenced by genes characteristic of the perturbed B cell axis in AA SLE patients.

Keywords: Autoimmunity; Rheumatology.

Figure 1. Individual SLE patients manifest varied patterns of signatures for 34 cell and process modules.

(A) GSVA was carried out on 17 female HC to determine the mean and SD of control GSVA scores for 34 cell type and process modules. HC mean scores ± 1 SD were used to determine a normal range for GSVA scores. SLE female patient (GSE88884 ILL1 and ILL2 data set cohorts; n = 1566) GSVA scores were determined and compared with HC values to determine whether patients had increased (+1), decreased (–1), or normal (zero) values. GSVA enrichment gene symbols for each module are in Supplemental Table 5. (B and C) Percentage of patients within each ancestry (AA, n = 216; NAA, n = 232; EA, n = 1118) with > 1 (B) or < 1 (C) SD GSVA scores for each cell type and process module. Fisher’s exact P < 0.05 are indicated by different color asterisk: black asterisks for comparisons between all 3, red asterisks between NAA and AA/EA, orange asterisks between NAA and EA, light blue asterisks between AA and EA, and dark blue asterisks between AA and NAA/EA. Exact P values and percentages are listed in Supplemental Table 6. (D) WGCNA was carried out on data set GSE88884 ILL1 and ILL2 cohorts separately. Pearson correlation r values to ancestry were determined for each module and listed if P < 0.05.

Figure 2. Gene expression differences in SLE patients are similar to ancestral gene expression differences in healthy controls.

(A) Limma DE analysis was carried out between HC AA and EA for 2 separate data sets (Supplemental Table 9). Increased (Up in AA) and decreased (Up in EA) transcripts were compared with 4 SLE cohorts of AA DE to EA. Overlap P values were all below 1 × 10^–22 for OR above 1. (B) GSVA for the 34 cell and process modules was carried out on healthy AA and EA subjects from 2 separate data sets. Welch’s t test was used to determine significant differences between ancestral GSVA scores; the mean and CI for the 10 GSVA scores significantly different (P < 0.05) between ancestries are shown.

Figure 3. Autoantibodies and complement levels were associated with gene expression profiles.

The mean difference in GSVA enrichment scores is shown for manifestations with significant (P < 0.05, Sidak multiple comparisons test) differences in enrichment scores as compared with all other manifestations. Asterisks indicate that, for dsDNA autoantibodies and low C, patients were compared with patients without either dsDNA autoantibodies (IU < 30) or low C (C3 > 0.8 g/L and C4 > 0.1 g/L). All patients in these analyses were positive for ANA. Number of patients with each SLEDAI component manifestation are shown in parentheses.

Figure 4. The higher number and different types of autoantibodies in AA SLE patients led to higher plasma cell, IGS, cell cycle, Treg, and myeloid-secreted signatures.

(A) Percentage of patients with different numbers of 5 autoantibodies (RNP, Sm, SSA, SSB, and dsDNA) by ancestry. (B) Comparison of plasma cell and IGS GSVA scores by the number of autoantibodies. (C) GSVA enrichment scores for all 34 cell and process modules were compared, for each autoantibody group, with patients of the same ancestry with 0 of 5 autoantibodies. Tukey’s multiple comparisons test was used to determine significant differences; 8 cell and process module signatures had significant differences (P < 0.05) between autoantibody⁺ and autoantibody^– groups. (D) GSVA enrichment scores had significant differences between autoantibody groups. (B and D) Dots represent single patient scores, and data are presented as mean ± SD. Numbers of patients in each group are shown in parentheses. The black dotted lines represent the mean ± 1 SD of the HC for GSVA scores. Tukey’s multiple comparisons test was used to determine if significant differences existed between GSVA scores for plasma cells and IFN signatures for each group, and P < 0.05 are shown. RNP⁺ and/or dsDNA⁺ autoantibody groups could also have Sm, SSA, or SSB autoantibodies in any combination.

Figure 5. Association of corticosteroid (CS) use and immunosuppressive therapy with changes in gene expression profiles.

Female SLE patients (1566 patients; GSE88884) were separated by ancestry and GSVA scores for each cell type or process module in patients receiving each therapy and were compared with GSVA scores for each cell type or process module in patients taking all other therapies. The patient numbers are in parentheses. Sidak multiple comparisons test was used to determine significant differences between therapies. The mean difference in GSVA score related to the treatment is shown for therapies with P < 0.05. Two EA patients were receiving cyclophosphamide and are included in the immunosuppressive (IS) calculation for EA.

Figure 6. Stepwise logistic regression analysis determined the importance of ancestry, SoC drugs, and SLEDAI components to the gene expression profile.

(A–D) CIRCOS visualization of odds ratios (OR) using stepwise logistic regression analysis for ancestry (A), serology (B), SoC drug (C), and time from onset of disease, age > 50, and SLE manifestation (D) to GSVA categories with P < 0.05 (P values, OR, and CI in Supplemental Tables 18–21). The thickness of the lines from the 26 variables to the GSVA categories represent the magnitude of the ORs. An interval graph was used to assign thickness of the lines where OR < 2, 1 pt; 2 ≥ OR < 3, 5pt; 3 ≥ OR < 10, 10pt; OR ≥ 10, 20pt. Red lines indicate OR above 1, and blue lines indicate OR below 1. OR between 0 and 1 are represented as 1/odds ratio to accurately reflect the magnitude of the negative relationship to the GSVA enrichment score.

Figure 7. A machine learning approach predicted AA from EA SLE patients and demonstrated the perturbed B cell axis in AA SLE.

(A) SLE patients were classified as AA using logistic regression, generalized linear models (GLM), and support vector machine (SVM) classifiers. ROC curve for logistic regression and the 2 different machine learning models in GSE88884 (ILL1 and ILL2 combined). (B) Top 25 gene predictors determined by SVM model.