key: cord-0019201-einb03a7
authors: Kingsmore, Kathryn M.; Bachali, Prathyusha; Catalina, Michelle D.; Daamen, Andrea R.; Heuer, Sarah E.; Robl, Robert D.; Grammer, Amrie C.; Lipsky, Peter E.
title: Altered expression of genes controlling metabolism characterizes the tissue response to immune injury in lupus
date: 2021-07-20
journal: Sci Rep
DOI: 10.1038/s41598-021-93034-w
sha: 8320bd3df863745468cef94fb31c6fdaeb6e003f
doc_id: 19201
cord_uid: einb03a7

To compare lupus pathogenesis in disparate tissues, we analyzed gene expression profiles of human discoid lupus erythematosus (DLE) and lupus nephritis (LN). We found common increases in myeloid cell-defining gene sets and decreases in genes controlling glucose and lipid metabolism in lupus-affected skin and kidney. Regression models in DLE indicated increased glycolysis was correlated with keratinocyte, endothelial, and inflammatory cell transcripts, and decreased tricarboxylic (TCA) cycle genes were correlated with the keratinocyte signature. In LN, regression models demonstrated decreased glycolysis and TCA cycle genes were correlated with increased endothelial or decreased kidney cell transcripts, respectively. Less severe glomerular LN exhibited similar alterations in metabolism and tissue cell transcripts before monocyte/myeloid cell infiltration in some patients. Additionally, changes to mitochondrial and peroxisomal transcripts were associated with specific cells rather than global signal changes. Examination of murine LN gene expression demonstrated metabolic changes were not driven by acute exposure to type I interferon and could be restored after immunosuppression. Finally, expression of HAVCR1, a tubule damage marker, was negatively correlated with the TCA cycle signature in LN models. These results indicate that altered metabolic dysfunction is a common, reversible change in lupus-affected tissues and appears to reflect damage downstream of immunologic processes.

. Dysregulation of metabolic gene signatures is common among lupus-affected tissues. (a) Comparison of DEGs among DLE, class III/IV LN GL, and class III/IV LN TI. (b) MCODE protein-protein interactions of common UP and DOWN DEGs were generated with Cytoscape using the STRING and ClusterMaker2 plugins and annotated with BIG-C functional categories (odds ratio (OR) > 1, p < 0.05) in Adobe Illustrator. Overlap p-value was calculated using Fisher's exact test. GSVA of signatures for (c) glycolysis, (d) the PPP, (e) the TCA cycle, (f) OXPHOS, (g) FAAO, (h) FABO, and (i) AA metabolism in lupus tissues and controls (CTLs). Each point represents an individual sample. Significant differences in enrichment of the metabolic signatures in each lupus tissue as compared to CTL was determined by Welch's t-test with Bonferroni correction. Numbers below each tissue indicate the number of lupus patients with enrichment scores 1 SD less than (< 1SD) or greater than (> 1SD) the CTL mean. For all calculations, the following sample numbers were used: DLE [CTL (n = 8), DLE (n = 9)], LN GL [CTL (n = 14), Cl III/IV (n = 22)], and LN TI [CTL (n = 15), Cl III/IV (n = 22) ]. **, p < 0.01; ***, p < 0.001; ****, p < 0.0001.

To determine whether cellular changes accounted for the decreased metabolic signatures, we examined enrichment of immune and non-hematopoietic cell signatures in lupus tissues. Increased tissue enrichment of immune cell signatures was variable, with DLE and LN GL demonstrating more enrichment of inflammatory cell signatures as compared to LN TI (Fig. 2a, Supplemental Fig. S1-S2 ). There was evidence of increased plasmacytoid dendritic cell (pDC), dendritic cell (DC), or monocyte/myeloid cell (monocyte/MC) signatures in some patients from each tissue. Non-hematopoietic cell gene expression was additionally altered. The endothelial cell (EC) signature was increased in LN GL. Keratinocyte transcripts were increased in 4/9 DLE patients, whereas melanocyte transcripts were decreased in 5/9 patients. General kidney cell transcripts www.nature.com/scientificreports/ were decreased in the majority of LN patients. Genes indicative of podocytes were decreased, whereas genes indicative of mesangial cells were increased in LN GL. Finally, the proximal tubule signature was significantly decreased and the collecting duct signature was significantly increased in LN TI. Although the monocyte/MC signature was consistently increased among lupus-affected tissues, its nature in each tissue varied. Linear regression of the monocyte/MC signature and FCN1 expression indicated DLE and LN GL contained inflammatory monocyte-derived macrophages 23 , but the correlation in LN TI was weak (Fig. 2b,  Supplemental Fig. S3a ). Tissue-resident macrophage populations were evaluated in DLE by expression of F13A1, FOLR2, SEPP1, and TXNIP 24 and in LN by co-expression of CD74 and CD81, which is characteristic of renal tissue-resident macrophages 25 . The tissue-resident macrophage signature was present in LN TI as evidenced by monocyte/MC correlation with both markers and co-expression of CD74 and CD81 (Fig. 2b,c, Supplemental  Fig. S3b) . However, the tissue-resident macrophage signature was not identified in DLE as shown by lack of correlation with F13A1, SEPP1, and TXNIP. Similarly, the tissue-resident macrophage signature was not found in LN GL as shown by lack of correlation with CD81 and no co-expression between CD74 and CD81. This suggests monocyte/MCs in DLE and LN GL are predominantly infiltrating inflammatory monocyte-derived macrophages, whereas LN TI is populated by tissue-resident macrophages.

Class II LN GL is molecularly similar to class III/IV LN GL. To elucidate whether the same metabolic signature changes were present in less severe LN, we expanded our analysis to incorporate WHO class II LN samples, where less immune cells have been observed histologically 26 . Unexpectedly, we found that genes controlling glycolysis, the TCA cycle, FAO, and AA metabolism were decreased in class II LN GL ( Fig. 3a-g) , and gene expression of pDCs and B cells was increased in class II ( Fig. 3h-o) . It is notable that gene expression of non-hematopoietic cell populations was altered in class II LN GL, including a decrease in kidney cell genes and an increase in EC genes ( Fig. 3p-t) . These abnormalities in expression of metabolic genes and non-hematopoietic cell genes were noted even in patients in which an increase in monocyte/MC genes was not detected (Supplemental Fig. S4 ). Altogether, even though inflammatory cell gene signatures were less pronounced in some patients relative to class III/IV, class II LN GL was molecularly similar to class III/IV LN GL (Fig. 3u) .

As in the glomerulus, class II LN TI was not statistically different from class III/IV LN TI (Supplemental Fig. S5 ); however, class II LN TI did not display decreased metabolic signatures. Class II LN TI had less indication of tissue cell damage and inflammatory infiltrate than class II LN GL (Supplemental Fig. S6 ).

Stepwise regression, which identifies the independent variables that best explain the dependent variable 27 , was employed to analyze cellular signature associations with metabolism gene signature changes. To improve precision, highly collinear cellular signatures in DLE were combined for stepwise analysis (Methods). In DLE, the inflammatory cell, EC, and keratinocyte signatures were positively correlated with the glycolysis signature, as indicated by a positive regression coefficient (Supplemental Fig. S7 , Supplemental Data S4). Conversely, the GC B cell and fibroblast signatures were negatively correlated with other metabolic signatures. The keratinocyte signature was also negatively correlated with the TCA cycle signature.

We then implemented both stepwise regression and classification and regression tree (CART) analysis, which partitions data by independent variables 28 , to determine the cellular signatures that were most associated with the metabolic signature changes in all classes of LN and controls. In LN GL, all metabolic signatures except for the PPP signature exhibited some dependence on the kidney cell signature by either stepwise regression or CART (Fig. 4 , Supplemental Data S5). The stepwise regression coefficients for the kidney cell signature were positive and of the largest magnitude for the TCA cycle, FAAO, FABO, and AA metabolism signatures; a negative regression coefficient was noted for the EC signature and glycolysis. Moreover, CART identified the EC signature as the strongest contributor to the glycolysis signature, with a positive EC GSVA score predicting a lower glycolysis GSVA score. Importantly, the EC signature classifier alone separated 23 of 30 LN GL samples from controls. It is notable that by CART analysis upregulation of the monocyte/MC signature tended to mitigate the decreased glycolysis signature related to ECs, suggesting that monocyte/MC genes might contribute to glycolysis in opposite directions from EC genes. The mocyte/MC signature exhibited a negative stepwise regression coefficient with the OXPHOS signature, which was also positively associated with the kidney cell signature by CART. Together, the stepwise regression and CART models indicate decreased OXPHOS is reflective of decreased kidney cell and increased monocyte/MC signatures, although increased pDC genes may also contribute.

The contribution of kidney-specific cell signatures to most metabolic changes was also observed in the TI (Supplemental Fig. S8 , Supplemental Data S6). In contrast to LN GL, EC genes did not contribute to the decreased glycolysis signature in LN TI, but instead the kidney cell signature was predicted by stepwise regression and CART. Numerous cell types were found to contribute to the overall TCA cycle signature. Stepwise regression demonstrated a positive contribution from proximal tubule, kidney cell, and granulocyte transcripts, and a negative contribution from fibroblast and pDC cell signatures. CART identified the proximal tubule signature as the strongest contributor, with a positive correlation to the TCA cycle signature, and confirmed the involvement of the fibroblast signature as well. Monocyte/MC transcripts contributed most to the OXPHOS signature, although the regression coefficient was negative, implying that the presence of monocyte/MC contributed to a decreased OXPHOS signature. The FABO signature, which has previously been shown to be defective in affected TI of multiple etiologies 11, 12, 29 , FAAO, and AA metabolism signatures also demonstrated a positive relationship with expression of kidney cell and/or proximal tubule genes.

We sought to confirm these findings by referencing data from single-cell RNA-sequencing (scRNA-seq) of LN biopsies 30 . In LN, tissue-resident macrophages exhibited decreased OXPHOS genes (Supplemental Fig. S9 ), which aligns with the negative relationship that stepwise regression and CART predicted between the monocyte/ Fig. S8a and e) . Conversely, CD4 T cells had both increased and decreased expression of glycolysis genes, making their metabolism difficult to interpret. The kidney epithelial cell cluster, which was comprised of general kidney cell as well as tubule genes, exhibited significant negative differential expression of two glycolysis genes (G6PC, GAPDH) and one TCA cycle gene (PDK4), but four OXPHOS genes were significantly over-expressed (MT-CO2, MT-CO3, MT-ND5, NDUFS3). This supports the idea that kidney cells have high oxidative metabolism, and their absence or dysfunction could contribute to decreased OXPHOS signatures. However, insufficient genes were detected by scRNA-seq to confirm the current findings determined with bulk RNA analysis.

Mitochondrial and peroxisomal signature changes and local hypoxia contribute to changes in metabolic gene expression in specific cells. As mitochondria and peroxisomes are the primary organelles responsible for metabolic processes such as OXPHOS and FAO 31 , we sought to examine whether there were detectable changes to organelle-specific gene expression that could explain the altered metabolic state. www.nature.com/scientificreports/ There was no evidence of global mitochondrial gene expression changes, although mitochondrial genes were decreased in approximately half of DLE patients, and mitochondrial transcription was increased in 15/30 LN GL patients ( Fig. 5a-e) . Conversely, the apoptotic mitochondrial changes signature was increased in both DLE and LN GL (Fig. 5f ). Genes associated with peroxisome biogenesis were decreased in some lupus patients in each tissue; in contrast, genes associated with peroxisomal fission were increased in some class III/IV LN TI patients ( Fig. 5g-h) . Stepwise regression analysis revealed a positive association between the apoptotic mitochondrial signature and the granulocyte and inflammatory cell signatures in DLE, and positive associations between the apoptotic mitochondrial signature and MC signatures in LN ( Fig. 5i -k, Supplemental Data S7-S9). Positive associations were observed between the peroxisome biogenesis signature and the kidney cell and proximal tubule signatures in LN GL and LN TI, respectively. Since hypoxia has been cited as a driver of kidney disease in the TI 32 , and a hypoxic microenvironment may result in decreased oxidative metabolism, we examined the contribution of HIF1A to metabolic gene signatures. Some lupus patients had increased expression of HIF1A by GSVA (Fig. 5l) . Although HIF1A expression did not significantly affect metabolic signatures in DLE, in LN TI, the HIF1A gene signature had a positive correlation with the glycolysis signature, whereas in LN GL, the HIF1A gene signature was positively associated with the PPP signature ( Fig. 5m -o, Supplemental Data S10-S12). This suggests hypoxia contributes to specific metabolic alterations in the kidney.

Metabolic gene expression changes occur independent of acute IFN stimulation. The IFN gene signature (IGS), a known hallmark of lupus 33, 34 and lupus-affected tissues 35 , has been implicated in metabolic alteration of MCs [36] [37] [38] [39] . To determine the functional relationship between IFN stimulation and metabolic alteration, we examined metabolic gene expression longitudinally in the IFNα-accelerated NZB/W murine model of LN, where the IGS increases at early timepoints following injection of IFNα adenovirus and then increases again when kidney disease develops (Fig. 6a ). To determine whether metabolic signatures were changed with the elevated IGS, we examined metabolic signatures after IFNα administration. No significant changes in metabolic signatures were observed at week 1 (W1) after IFNα administration ( Fig. 6b-h) , the peak of early IGS expression. In contrast, the TCA cycle, OXPHOS, FAO, and AA metabolism gene signatures were decreased at W7 post-IFNα administration, when LN develops 40 . There were negative correlations between the IGS and metabolic signatures for the TCA cycle, OXPHOS, FAO, and AA metabolism, suggesting increased IFN could contribute to decreased metabolism. However, these is a clear separation of the W1 and W7/W9 datapoints in the correlation analysis for TCA cycle, OXPHOS and AA metabolism signatures, in which W1 datapoints show positive GSVA scores for both the IGS and metabolic signatures. Thus, based upon the lack of transcriptional metabolic changes following early IFN exposure in the mouse, it does not appear that acute exposure to type I IFN explains the metabolic changes in lupus tissues, as transcriptional changes to metabolism occur when disease develops and the IFN signature reoccurs, but not at W1, W2, or W3

To determine the robustness of the observed metabolic and cellular gene expression changes and determine whether dysregulation could be reversed by immunosuppressive therapy, we analyzed gene expression in pre-and post-treatment kidneys of lupus mice (Supplemental Data S1). Although some metabolic changes had previously been identified in the kidneys of untreated and treated NZB/W and NZM2410 mice 41 , no analysis of the nature of the affected cells was carried out.

Metabolic gene expression was significantly altered in NZM2410, NZB/W, IFNα-accelerated NZB/W (GSE72410), and NZW/BXSB mice. Treatment of NZM2410 mice with BAFF-R-Ig and proteinuric NZB/W mice with a combination of cyclophosphamide (CTX) + CTLA4-Ig + anti-CD154 restored TCA cycle, FAAO, FABO, and AA metabolism gene expression ( Fig. 7a, b) . Notably, in the NZB/W kidneys when combination therapy was discontinued and LN relapsed, some metabolic abnormalities recurred. In IFNα-accelerated NZB/W mice, both the glycolysis and TCA cycle signatures were significantly decreased with the onset of LN (IFN W7), and CTX treatment significantly restored TCA cycle gene expression to baseline pre-disease levels (Fig. 7c) . In MRL/lpr mice, the same trends were observed, although the data did not achieve statistical significance (Fig. 7d) . NZW/BXSB mice had decreased metabolic signatures for all processes except glycolysis, the PPP, and OXPHOS (Fig. 7e) .

GSVA of cellular changes in murine LN models demonstrated similar results to human LN, although inflammatory infiltrate and changes in the EC and podocyte signatures were less robust (Supplemental Figs. S10-S15). BAFF-R-Ig, combination therapy, and CTX-treatment restored kidney cell and proximal tubule gene signatures, and combination therapy decreased inflammatory cell gene signatures. Correlation analysis in all murine LN models showed strong positive correlations between metabolism and proximal tubule transcripts, whereas negative correlations were seen between metabolism and inflammatory cell signatures (Supplemental Fig. S16 ). This suggests that treatment may allow for functional metabolic recovery of non-hematopoietic cells in the kidney in the presence or absence of changes to the inflammatory cell signal.

Metabolic changes correlate with expression of genes indicating tubular damage. Finally, we examined whether changes in genes controlling metabolism occurred synchronously with changes in expression of HAVCR1 (KIM1) and LCN2, which are known markers of tubular damage 42, 43 . Although HAVCR1 expression was increased in class III/IV human LN TI patients, there were no significant changes to LCN2 in any class of LN TI (Fig. 8a, b) www.nature.com/scientificreports/ the kidney cell, proximal tubule, and TCA cycle signatures in both human and murine LN, although there were model-dependent differences in magnitude (Fig. 8c-f, Supplemental Fig. S17 ).

Multi-pronged bioinformatic analyses of gene expression data from human lupus tissues revealed that despite intra-tissue heterogeneity metabolic dysfunction was present in all tissues. Immune effector cells have high metabolic needs 13, 14, 44, 45 and, therefore, we initially hypothesized immune infiltration may be responsible for the observed lupus tissue-wide metabolic dysregulation. Although kidney-infiltrating CD8 T cells in murine LN are functionally exhausted with defective mitochondria 16 , anergic/activated T cell markers were not found in these human LN samples, and regression models indicated T cells contributed minimally to changes in renal metabolic gene expression. Similarly, monocyte/MCs, which were increased in some patients from all tissues, might be expected to contribute to enhanced glucose metabolism -either glycolysis (M1 macrophages) or OXPHOS (M2 macrophages) 44 . Indeed, monocyte/MC signatures were inversely correlated with OXPHOS in both LN tissues, and positively correlated with glycolysis in LN GL, suggesting they are likely M1 in nature and may contribute to the altered metabolic landscape of intact tissues. Although gene expression revealed differing origins of the renal MC populations, as they reflect monocyte-derived macrophages in LN GL and tissue-resident macrophages in LN TI, the consistent MC presence aligns with their prominent role in tissue damage 46 . Observed increases in the monocyte/MC signature and strong inverse correlations between MCs and metabolism may reflect the role of MCs in tissue damage, even when T and B cells are not yet abundant.

To examine whether metabolic abnormalities represented primarily tissue cell defects as opposed to changes in inflammatory cell metabolism, we analyzed gene expression in class II LN, in which less inflammatory infiltrate is evident histologically 26 . Even though class II LN samples had evidence of increased immune/inflammatory cell signatures that coincided with changes to metabolic signatures, examples were observed in which the changes in metabolic signatures were found in the absence of a monocyte/MC signature, suggesting alterations in metabolic signatures may be initiated immediately following immune complex (IC) deposition and complement activation. Subsequent monocyte/MC and other inflammatory cell activation/infiltration may then contribute to further damage of tissue cells. Indeed, previous studies report that changes in kidney gene expression can occur following early IC deposition, but before microscopic detection of inflammation 47 , consistent with our transcriptomic results in class II LN.

Our findings of altered metabolism in lupus tissues align with changes seen in other forms of tissue pathology. We observed a positive association between the keratinocyte and glycolysis signatures, and upregulation of glycolysis has been observed in keratinocytes during cutaneous infection 48 . Increased glycolysis 49, 50 and decreased PPAR signaling, TCA cycle, and OXPHOS have been reported in dermal fibroblasts subjected to radiation 50 , indicating dermal fibroblasts have the potential to contribute to inflammation-induced alterations in metabolism. However, in the current study, stepwise regression indicated that the fibroblast signature was negatively associated with the glycolysis, PPP, and OXPHOS signatures, whereas associations with FAO did not achieve statistical significance. The negative correlation between the fibroblast and PPP signatures contrasts with the observed upregulation of the PPP in cultured fibroblasts and keratinocytes that were exposed to UV-induced oxidative stress 8 . This suggests that in vivo fibroblasts are altered by signals different than that mediated by UV light, as expected since fibroblasts are deep in the dermis and shielded from such ambient exposure.

Metabolic dysregulation is also common in kidney disease 12, 15, 29, 51, 52 . In non-diabetic chronic kidney disease, TCA cycle abnormalities measured in urine metabolites coincided with changes to kidney gene expression 51 , supporting our conclusions that metabolic dysregulation primarily reflects altered renal cell function as opposed to changes in immune cell metabolism. Moreover, defects in FAO have been correlated with fibrosis progression in TI disease 11 . Both human and murine models with TI fibrosis exhibited decreased expression of FAO enzymes and resultant increases in lipid deposition, which was reversed by correcting the metabolic abnormalities 11 . We similarly found that FAO signatures aresubstantially decreased in the TI; however, regression models indicated that altered FAO transcripts were most associated with decreased kidney cell signatures. Although fibroblasts were not predicted by the models, fibrosis or fibroblast enrichment could contribute to decreased kidney cell and proximal tubule transcripts in the TI.

In the glomerulus, ECs appear to play an additional role in disease. Endothelial activation in LN has been suggested previously 53 , and EC transcripts were increased in 83% of all LN GL samples. Increased EC transcripts may reflect altered EC function, potentially resulting from cytokine/growth factor stimulation and/or hypoxia-induced cellular damage. Abnormal angiogenesis has been reported in diabetic nephropathy resulting in leaky vessels 54 , but the function of glomerular ECs in LN is less well-defined. Indeed, glomerular ECs have been found to be dependent upon podocyte stimulation for differentiation 55 , whereas other studies suggest EC damage precedes podocyte injury 56 . Our findings from class II LN support the latter, as signature changes to ECs occurred in the absence of podocyte changes in some patients, suggesting that ECs are early participants in LN.

The relationship between the EC signature and metabolic gene expression changes implied an alteration in EC physiology in LN. Although healthy ECs are highly glycolytic 57 , all regression techniques in LN GL indicated an inverse correlation between the EC and glycolysis signatures. Consistent with our findings, in diabetic nephropathy stalled glycolytic flux has been observed in ECs 57 . These data suggest that glomerular ECs are metabolically altered, perhaps because of IC and complement stimulation and/or cytokine exposure, making them less capable of maintaining normal function. Indeed, in the GL, the EC signature had a positive regression coefficient with the FABO signature, and quiescent ECs have been reported to upregulate FABO 58 , supporting the conclusion that ECs are functionally deranged in LN GL.

Analysis of metabolism-associated genes in cell clusters derived from scRNA-seq of LN biopsies 30 further supports our finding that changes in metabolism are most closely related to kidney cell gene expression, with www.nature.com/scientificreports/ minor contributions from resident or infiltrating immune cells. Tissue-resident macrophages exhibited decreased OXPHOS genes, whereas CD4 T cell metabolism was unclear. Notably, the kidney epithelial cell cluster reported many metabolism-associated genes, suggesting decreased glycolysis and TCA cycle, but increased OXPHOS. Proximal tubules, which have the most mitochondria of any kidney epithelial cell, are known to be dependent on oxidative metabolism 59 , and this further supports the idea that diminished OXPHOS in bulk RNA in part reflects decreased kidney epithelial cell transcripts. Additionally, because scRNA-seq looks at expression of individual cells as opposed to the bulk environment, the detected kidney epithelial cells are likely the residual functionally normal ones. However, because of technical issues including cell yield and read depth, there are difficulties in determining the status of individual cell metabolism from the scRNA-seq data. Altogether, scRNA-seq appeared to be less effective than deconvolution of bulk RNA to detect important but subtle changes in cellular metabolism.

To determine whether defective organelles were responsible for metabolic alteration, we examined gene expression specific for both mitochondrial and peroxisomal function. We observed changes to the mitochondrial and peroxisomal gene signals in some patients, and correlation analysis suggested these changes were associated with specific cell types. Notably, the peroxisome biogenesis signature was positively associated with signatures for ECs, kidney cells, and proximal tubules. Moreover, as the kidney is highly susceptible to hypoxia 60 , we investigated the propensity for hypoxia to contribute to altered metabolism. Although GSVA demonstrated no significant increases in HIF1A expression, there was an association between HIF1A and the PPP and glycolysis signatures in LN GL and LN TI, respectively, suggesting that hypoxia can have specific effects on metabolism, and may contribute to some of the metabolic changes observed in the tissues.

The relationship between the IGS and metabolic signature changes in IFNα-accelerated LN mice indicated that acute type I IFN exposure could not explain the observed metabolic changes. The mouse studies were particularly informative because IFNα exposure was regulated. Although it has been demonstrated that type I IFN stimulation can increase OXPHOS and FAO in DCs and MCs 36, 38 , inhibit isocitrate dehydrogenase (part of the TCA cycle) in macrophages 39 , and alter oxidative metabolism in other cells 61 , IFNα overexpression did not change metabolic signatures at early timepoints. However, there was an inverse relationship between the IGS and metabolic defects after LN onset, when the IGS recurred. These results support the conclusion that downregulation of metabolic pathways is unlikely to be explained by the known actions of type I IFN alone, but rather during LN progression, decreased metabolic signatures may be parallel reflections of continued IGS exposure and inflammation.

Metabolic gene expression was altered in four murine LN models and immunosuppressive treatment, not known to directly affect cellular metabolism, restored metabolic gene expression. Although combination therapy diminished inflammatory cell abundance in the kidneys of NZB/W mice, we also observed restoration of metabolic and kidney cell gene expression after treatment in models with little inflammatory infiltrate or those in which inflammatory cells were not significantly decreased with treatment. This suggests that although inflammatory cells play a critical role in mediating tissue cell damage and metabolic dysfunction, damage is not related only to local inflammatory cells, as intensive therapy with CTX restores tissue cell defects without significant changes to inflammatory cells. Importantly, these results demonstrate that metabolic abnormalities in tissue cells are reversible with immunosuppressive therapy and restored metabolic gene expression might be considered a goal of effective lupus treatment. Moreover, monitoring tissue metabolism may be especially important in situations where anti-metabolic drugs, such as metformin or 2DG, are employed. Whereas these drugs may be promising for correction of individual immune cell defects, they have potential consequences for already metabolically deranged tissue cells.

Additionally, these studies reveal subtle differences in pathology of glomerular and tubulointerstitial involvement in LN. In both kidney regions, we observed decreased resident non-hematopoietic cell signatures and increased monocyte/MC signatures. However, monocyte/MCs in the glomerulus were likely monocyte-derived, whereas those in the tubulointerstitium appeared more like tissue-resident macrophages. Moreover, tubulointerstitial diseases in class III/IV LN was characterized by less inflammatory infiltrate than was observed in the glomerulus, and although metabolic signatures were comparably decreased in all classes of glomerular LN, metabolic signature changes in class II tubulointerstitial LN were less consistently regulated. These results align with studies that show tubulointerstitial damage occurs later in LN and predicts end stage renal disease 62 . Poorer outcomes in class III/IV could be related to persistence of abnormalities or inhibition of repair mechanism that might contribute to progressive renal disease. Nonetheless, it is noteworthy that even the more modest immunologic damage in class II LN was associated with marked changes in metabolic signatures.

Detectable changes to immune cell, EC, kidney cell, and metabolic gene signatures in all classes of LN GL is notable. We found that gene expression may detect cellular changes with greater sensitivity than immunohistochemistry, when little or no inflammatory infiltrate is observed histologically. Gene expression may provide an advance to current classification or diagnostic techniques, as gene expression changes are detectable before discernable immunohistochemical changes, and transcriptomic analysis of metabolism can elucidate potential functional rather than merely histopathologic changes.

Our study is not without limitations. We performed post-hoc analysis of bulk gene expression in three lupusaffected tissues comprising limited numbers of lupus patients. Moreover, 37.5% of LN patients were being treated with immunosuppresives 53 , that may have affected gene signatures. Additionally, the gene signature we identified as reflecting general kidney cells may be more specific for tubule cells, despite the strong representation in the glomerulus. Furthermore, regression analyses provided an estimate of the cellular variables that are most associated with each metabolic signature, but accuracy can be limited by sample size, and there is a chance of overfitting. Future work with larger cohorts currently not available would be necessary to validate these results. Additionally, direct assessment of functional metabolism may be necessary to assay how metabolic changes at the gene expression level reflect changes in protein content and cellular function.

In conclusion, prominent alterations in cellular metabolism signatures are characteristic of lupus tissue pathology. Systems bioinformatics and assessment with regression modeling techniques revealed that the monocyte/ www.nature.com/scientificreports/ MC signature, including both monocyte-derived macrophages and tissue-resident macrophages, was increased in many lupus patients, kidney cell signatures were decreased in LN, the EC signature was increased in LN GL, and these cell signature changes were associated with altered metabolism signatures. Moreover, apoptotic mitochondrial gene changes were associated with MC genes in DLE and LN GL. In murine LN, metabolic dysregulation correlated with tubular damage marker expression and metabolic gene changes were reverted to normal by immunosuppressive therapy. Altogether, altered metabolism may serve as a promising biomarker or therapeutic target for lupus tissue disease, especially as metabolic gene expression changes precede expression of the renal damage biomarker LCN2 in human LN, and coincide with changes to HAVCR1. Indeed, urinalysis has been used to measure metabolite biomarkers in the kidney 51, 52 and metabolism transcripts could be used to estimate degree of kidney cell damage and assess treatment efficacy. Although treatment strategies aimed at metabolic restoration are not straightforward, the current findings support the conclusion that immunosuppressive therapy can restore metabolic function, and, thereby, may ameliorate damage in specific lupus-affected tissues.

Human and mouse gene expression datasets. Raw data from publicly available human and murine lupus datasets were derived from the Gene Expression Omnibus (GEO) repository. All datasets are summarized in Supplemental Data S1. GSE72535 comprises microarray analysis of lesional skin biopsies from human patients with DLE with no systemic involvement with Cutaneous Lupus Activity and Severity Index (CLASI) ≥ 2 63 . Some DLE patients were treated with various therapies including corticosteroids, immunomodulators, and hydroxychloroquine 63 .

GSE32591 consists of microarray analysis of human renal biopsies that were originally derived from the European Renal cDNA Bank (ERCB) 53 . LN patients from the ERCB (n = 32) had an average age of 35.1 ± 2.4 years, average proteinuria of 2.9 ± 0.6 g/day, and eGFR MDRD of 63.7 ± 5.4 ml/min/1.73m 2 . LN patients were treated with various therapies, including steroids and immunosuppressants 53 . Two patients from this group were excluded from our analyses because they had non-inflammatory class V LN.

GSE86423 consists of samples from the IFNα-accelerated NZB/W LN model. Female NZB/W mice were injected with an adenovirus vector expressing recombinant murine IFNα (5 × 10 9 particles) at 9 weeks 40 . Kidney gene expression was measured in mice at 0, 1, 2, 3, 4, 5, 7, and 9 weeks post IFNα injection 40 .

GSE32583 and GSE49898 consist of samples from three murine lupus models. Kidney gene expression was measured in NZM2410 mice including 6-8 week pre-disease mice, 22-30 week diseased mice, and treated mice in remission (Tx + 15w) 41, 53 . NZM2410 mice in the Tx + 15w group were treated with adenovirus expressing BAFF-R-Ig at 22 weeks, and then sacrificed at 30-35 weeks or 55 weeks 41 . Kidney gene expression in NZB/W mice was measured at 16w, 23w, 36w, and after treatment 41, 53 . Both 23w and 36w mice shown in the main figures had confirmed proteinuria. Some NZB/W mice with proteinuria (> 300 mg/dl) at two timepoints were treated with combination therapy -one dose of 50 mg/kg CTX, six doses of 100 µg CTLA-4-Ig, and one dose 250 µg of anti-CD40L 41 . Mice were determined to be in remission if they achieved proteinuria of ≤ 30 mg/dl at two timepoints 41 . One group was sacrificed 3-4 weeks after remission (Tx Rem. + 3-4w) and another was sacrificed > 5-14 weeks after remission (Tx Rem.+ > 5w) 41 . The latter group had confirmed histologic relapse 41 . Kidney gene expression from NZW/BXSB mice was measured at 17w (prenephritic mice) or 18-21w mice with confirmed proteinuria 53 .

GSE72410 comprises samples from the IFNα-accelerated NZB/W LN model. NZB/W mice were treated with adenovirus-expressing murine IFNα at 14-15w (1.2 × 108 IFNα rAd5-CMV) 64 . Kidney gene expression was measured in 17w naïve NZB/W (naïve), as well as IFNα-accelerated NZB/W mice: 17w mice 21 days post IFNα injection (IFN W3), and in mice treated with vehicle or CTX for four weeks beginning three weeks post IFNα njection (IFN W7 + Veh or IFN W7 + CTX) 64 .

GSE153021 consists of samples from the MRL/lpr LN model. MRL/lpr mice were treated with vehicle, prednisone, mycophenolate mofetil (MMF), FK506, or all three (Multi-target) for eight weeks beginning at week 8 or age-matched wildtype MRL/MpJ mice treated with vehicle 65 .

Quality control and data normalization. Microarray data were processed as previously described 66 using free, open source programs (GEOquery, affy, affycoretools, limma, and simpleaffy). Unnormalized arrays were inspected for visual artifacts or poor RNA hybridization using QC plots. Datasets were annotated using their native chip definition files (CDFs). Probes missing gene annotation data were discarded. Raw data (CEL files) from the Affymetrix platform were background corrected and normalized using guanine cytosine robust multiarray average (GCRMA) or robust multichip average (RMA) algorithms, whereas raw data files from Illumina chip were read and normalized using neqc (limma R package).

RNA-seq data (GSE72410 and GSE153021) was processed from FASTQ files as previously described by Daamen et al. 67 and also described below. SRA files were downloaded and converted into FASTQ format. Read ends and adapters were trimmed with Trimmomatic (v0.38) using a sliding window, ilmnclip, and headcrop filters. The reads were head cropped at 6 bp and adapters were removed before read alignment. Reads were mapped to the mouse reference genome m17 using STAR, and the .sam files were converted to sorted .bam files using Sambamba. The mouse reference genome was downloaded from GENCODE. Read counts were summarized using the featureCounts function of the Subread package (v1.61). The DESeq2 workflow was used to filter RNA-seq genes with low expression (i.e. genes with very few reads). The filtered raw counts were normalized using DESeq method and then log2 transformed.

Principal component analysis was used to inspect the raw data files from each dataset for outliers. All log2 transformed data was formatted into R expression set objects (E-sets). 73, 74 . The keratinocyte signature was derived from the keratinocyte-specific genes of Gazel et al. 75 . Kidney-specific lists generated from the Human Protein Atlas and single-cell data were additionally modified to incorporate genes found by both transcriptomics and immunohistochemistry 76 . The mesangial cell signature was derived from PanglaoDB 77 . In all tissues, the hematopoietic, EC, and fibroblast gene sets were evaluated. For non-hematopoietic gene signatures, those relevant to each tissue were reported (DLE-keratinocyte, melanocyte; LN GL-kidney cell, podocyte, mesangial cell; LN TI-kidney cell, proximal tubule, Loop of Henle (LoH) cell, distal tubule, collecting duct cell). Mitochondrial and peroxisomal gene sets were derived from literature mining and the BIG-C, and also compared to signatures in the MSIG database. The Apoptotic Mitochondrial Changes signature (M7482) was accessed on the MSIG database 78, 79 and is derived from GO_0008637. GO_Mitochondrial_Fission (M12786) and GOBP_Peroxiome_Fission (M22828) signatures were accessed on the MSIG database and modified slightly. The IGS is the type I IFN core signature previously reported 35 .

All human metabolism, non-hematopoietic cell, and IFN gene sets were converted to murine gene sets using the homologene R package and human2mouse function. Genes that were not converted programmatically were manually converted by GeneCards and Mouse Genome Informatics orthologs. Murine hematopoietic cell gene sets were curated by literature mining. For murine datasets with expression of two or more anergic/activated T cell markers, the aergic/activated T cell signature was combined with the T cell signature for GSVA analysis.

Although the same gene sets were input for each category in each tissue (Supplemental Data S3), the genes used in calculation of the GSVA enrichment score for each tissue differ slightly based upon the gene measurement platform and expression within that sample. All reported GSVA enrichment scores, except for the HIF1A gene signature, were calculated based upon a minimum of three genes. The HIF1A gene signature only includes HIF1A because there were no additional genes determined to be specific to hypoxia only. Although two genes were well co-expressed, the DC signature in LN TI did not meet the minimum three gene requirement for GSVA, nor did the LDG signature in LN GL and LN TI.

Hierarchical clustering. Human lupus tissue samples were hierarchically clustered by the Euclidian distance of their GSVA enrichment scores into two (k = 2) or four (k = 4) clusters using the heatmap.2 function in R.

Regression models. For all linear models, GSVA scores for cellular signatures in all tissue samples were input as independent variables, and the pathway GSVA score (metabolism signature, mitochondrial signature, or peroxisomal signature) was input as the dependent variable. As GSVA scales the expression of a signature from -1 to -1, the value for each input cellular or metabolic signature in each sample is relative to the same signature in other samples and to the other signatures in the same sample. To ensure that collinearity between immune cell signatures did not confound results for stepwise regression analyses in DLE, we combined the pDC, skin-specific DC, monocyte/MC, T cell, anergic/activated T cell, B cell, and plasma cell signatures into the "inflammatory Scientific Reports | (2021) 11:14789 | https://doi.org/10.1038/s41598-021-93034-w www.nature.com/scientificreports/ cell" signature, because the genes were highly co-expressed. The list of all genes used as the "inflammatory cell" signature can be found in Supplemental Data S3. This reduced the number of input variables for DLE stepwise analysis to ten cell signatures.

Visualization. Final figures were generated in GraphPad Prism or Adobe Illustrator.

Overlap p-values and ORs for functional enrichment of DEGs were calculated in R using twosided fisher.test with confidence level = 0.95. Because of known heterogeneity in lupus patients, the number of lupus patients in each tissue who fell above or below the control mean ± 1 standard deviation (SD) were then reported in order to determine whether individual patients exhibit an increased or decreased signature when the population did not achieve statistical significance. Calculation of mean and standard deviation (SD) for the control samples for each GSVA score in each tissue was performed in Microsoft Excel. All analyses in GraphPad Prism were carried out with version 8.2.0 (435) or later versions. Control and lupus sample populations of GSVA scores for each gene set were assessed for normality using the D' Agostino-Pearson test in GraphPad Prism and the distributions for 75% or more of the gene sets in each population were determined to be normal. Welch's t-test with Bonferroni correction for GSVA enrichment in human samples was performed in GraphPad Prism. Bonferroni correction for metabolic signatures, immune cell signatures, non-hematopoietic signatures, or mitochondrial/peroxisomal signatures were performed separately. The Mann-Whitney U test for GSVA enrichment in murine samples was performed in GraphPad Prism. Univariate (simple) linear regression and Pearson correlation analyses were carried out in GraphPad Prism. Hedges' g effect sizes were calculated in R using the cohen.d function with Hedges' correction under the "effsize" package. Stepwise regression was performed in R using the lm function followed by the stepAIC function. Variance inflation factor (VIF) < 10 was confirmed for each independent variable (Supplemental Data S4-S12). Although some data points were determined to be influential to the stepwise equation using Cook's D, no samples were removed from the models in efforts to capture the heterogeneity present in lupus. CART analysis was performed in R using the rpart function with the "anova" method. Each resulting decision tree was pruned once except for the glycolysis and PPP signatures in LN TI. GSVA scores for all individual samples (patient or control) are presented as individual data points in either dot or violin plots. The number of samples for each group can be found in Supplemental Data S1 and the figure legends. Information regarding the statistical comparisons made and level of significance is mentioned in the figure legends.

All microarray datasets in this publication are available on NCBI's GEO database (Supplemental Data S1). All bioinformatic software used in this publication is open source, and freely available for R. Example codes used in this paper (LIMMA, GSVA, stepwise regression, and CART) are available at figshare, https:// www. figsh are. com as "AMPEL BioSolutions LIMMA Differential Expression Analysis Code, " "AMPEL BioSolutions GSVA AFFY nonzeroIQR Code", and "AMPEL BioSolutions Stepwise and CART Code, " respectively.

Received: 4 January 2021; Accepted: 17 June 2021

Pathogenesis and targeted treatment of skin injury in SLE

Pathogenesis of kidney disease in systemic lupus erythematosus

Immunometabolism in systemic lupus erythematosus

Normalization of CD4+ T cell metabolism reverses lupus

Metabolic Factors that Contribute to Lupus Pathogenesis

Effects of inflammation on cholesterol metabolism: Impact on systemic lupus erythematosus

Mitochondrial dysfunction: A neglected component of skin diseases

Acute activation of oxidative pentose phosphate pathway as first-line response to oxidative stress in human skin cells

Dysregulated bioenergetics: a key regulator of joint inflammation

Global metabolic profiling of human osteoarthritic synovium

Defective fatty acid oxidation in renal tubular epithelial cells has a key role in kidney fibrosis development

Alteration of fatty acid oxidation in tubular epithelial cells: From acute kidney injury to renal fibrogenesis

Metabolic reprogramming in macrophages and dendritic cells in innate immunity

Metabolic Regulation of T Lymphocytes

Metabolic pathways and immunometabolism in rare kidney diseases

Kidney-infiltrating T cells in murine lupus nephritis are metabolically and functionally exhausted

Lipids rule: resetting lipid metabolism restores T cell function in systemic lupus erythematosus

T cell metabolism: new insights in systemic lupus erythematosus pathogenesis and therapy

Targeting metabolism for lupus therapy

Personalized immunomonitoring uncovers molecular networks that stratify lupus patients

Patient ancestry significantly contributes to molecular heterogeneity of systemic lupus erythematosus

GSVA: gene set variation analysis for microarray and RNA-Seq data

Expression of M-Ficolin in human monocytes and macrophages

Transcriptome landscape of myeloid cells in human skin reveals diversity, rare populations and putative DC progenitors

Single-cell RNA sequencing identifies candidate renal resident macrophage gene expression signatures across species

The classification of glomerulonephritis in systemic lupus erythematosus revisited

Robust stepwise regression

Classification and regression trees

Impaired B-oxidation and altered complex lipid fatty acid partitioning with advancing CKD

The immune cell landscape in kidneys of patients with lupus nephritis

Metabolic interplay between peroxisomes and other subcellular organelles including mitochondria and the endoplasmic reticulum. Front

Hypoxia and Hypoxia-Inducible Factors in Kidney Injury and Repair

Interferon-inducible gene expression signature in peripheral blood cells of patients with severe lupus

Interferon and granulopoiesis signatures in systemic lupus erythematosus blood

Gene expression analysis delineates the potential roles of multiple interferons in systemic lupus erythematosus

Transcriptional profiling suggests extensive metabolic rewiring of human and mouse macrophages during early interferon alpha responses

Direct type I IFN but not MDA5/TLR3 activation of dendritic cells is required for maturation and metabolic shift to glycolysis after poly IC Stimulation

Type 1 interferons induce changes in core metabolism that are critical for immune function

Autocrine IFN-I inhibits isocitrate dehydrogenase in the TCA cycle of LPS-stimulated macrophages

Pristane-accelerated autoimmune disease in (SWR X NZB) F1 mice leads to prominent tubulointerstitial inflammation and human lupus nephritis-like fibrosis

Identification of stage-specific genes associated with lupus nephritis and response to remission induction in (NZB × NZW)F1 and NZM2410 mice

Kidney injury marker 1 and neutrophil gelatinase-associated lipocalin in chronic kidney disease

Lipocalin 2 is essential for chronic kidney disease progression in mice and humans

The metabolic signature of macrophage responses

Initial B cell activation induces metabolic reprogramming and mitochondrial remodeling

What is damaging the kidney in lupus nephritis?

Activated renal macrophages are markers of disease onset and disease remission in lupus nephritis

Metabolic stress drives keratinocyte defenses against staphylococcus aureus infection

AMPK-mediated increase of glycolysis as an adaptive response to oxidative stress in human cells: Implication of the cell survival in mitochondrial diseases

Metabolic regulation of dermal fibroblasts contributes to skin extracellular matrix homeostasis and fibrosis

Metabolomics and gene expression analysis reveal down-regulation of the citric acid (TCA) cycle in non-diabetic CKD patients

Metabolomics reveals signature of mitochondrial dysfunction in diabetic kidney disease

Cross-species transcriptional network analysis defines shared inflammatory responses in murine and human lupus nephritis

Abnormal angiogenesis in diabetic nephropathy

Glomerular endothelial cell differentiation

Glomerular endothelial cell injury and damage precedes that of podocytes in adriamycin-induced nephropathy

Endothelial cell metabolism in normal and diseased vasculature

Quiescent endothelial cells upregulate fatty acid β-oxidation for vasculoprotection via redox homeostasis

Mitochondrial energetics in the kidney

Role of hypoxia in the pathogenesis of renal disease

Effects of interferons and viruses on metabolism

Tubulointerstitial damage predicts end stage renal disease in lupus nephritis with preserved to moderately impaired renal function: A retrospective cohort study

A subset of CD163+ macrophages displays mixed polarizations in discoid lupus skin

Btk-specific inhibition blocks pathogenic plasma cell signatures and myeloid cell-associated damage in IFNαdriven lupus nephritis

Transcriptomic analysis uncovers novel synergistic mechanisms in combination therapy for lupus nephritis

Identification of alterations in macrophage activation associated with disease activity in systemic lupus erythematosus

Comprehensive transcriptomic analysis of COVID-19 blood, lung, and airway

Evolving gene/transcript definitions significantly alter the interpretation of GeneChip data

Cytoscape: A software Environment for integrated models of biomolecular interaction networks

ClusterMaker: A multi-algorithm clustering plugin for Cytoscape

Microarray data processing techniques for genome-scale network inference from large public repositories

Analysis of gene expression from systemic lupus erythematosus synovium reveals myeloid cell-driven pathogenesis of lupus arthritis

Analysis of the human tissue-specific expression by genome-wide integration of transcriptomics and antibodybased proteomics

Tissue-based map of the human proteome

Transcriptional profiling of epidermal keratinocytes: Comparison of genes expressed in skin, cultured keratinocytes, and reconstituted epidermis, using large DNA microarrays

The kidney transcriptome and proteome defined by transcriptomics and antibody-based profiling

PanglaoDB-A Single Cell Sequencing Resource For Gene Expression Data

Molecular signatures database (MSigDB) 3.0

Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles

Transcriptomic meta-analysis of lupus affected tissues reveals shared immune, metabolic, and biochemical dysregulation

The authors would like to thank Anthony Yeo, Erika Hubbard, Adam Labonte, and Kate Owen for helpful discussion about the data and analysis techniques. Furthermore, the authors would like to thank Nicholas Geraci and Pankhuri Wanjari for help with DEG analysis version control. The authors would additionally like to acknowledge presentation of a preliminary portion of this work at the 2019 American College of Rheumatology Annual Meeting 80 .

Correspondence and requests for materials should be addressed to K.M.K.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.