Multi-organ single-cell analysis reveals an on/off switch system with potential for personalized treatment of immunological diseases

Summary Prioritization of disease mechanisms, biomarkers, and drug targets in immune-mediated inflammatory diseases (IMIDs) is complicated by altered interactions between thousands of genes. Our multi-organ single-cell RNA sequencing of a mouse IMID model, namely collagen-induced arthritis, shows highly complex and heterogeneous expression changes in all analyzed organs, even though only joints showed signs of inflammation. We organized those into a multi-organ multicellular disease model, which shows predicted molecular interactions within and between organs. That model supports that inflammation is switched on or off by altered balance between pro- and anti-inflammatory upstream regulators (URs) and downstream pathways. Meta-analyses of human IMIDs show a similar, but graded, on/off switch system. This system has the potential to prioritize, diagnose, and treat optimal combinations of URs on the levels of IMIDs, subgroups, and individual patients. That potential is supported by UR analyses in more than 600 sera from patients with systemic lupus erythematosus.


In brief
Lilja et al. report that treating immune diseases is complicated by involvement of thousands of genes. Combined analyses of mouse arthritis and human immune diseases show organome-, cellulome-, and genome-wide changes, which are switched on or off by pro-or anti-inflammatory upstream regulators (URs). Targeting URs may contribute to personalized medicine. ''I never feel completely well.'' This is a common complaint from patients with immune-mediated inflammatory diseases (IMIDs), despite state-of-the art treatment. This sentiment reflects a general health care problem: according to the US Food and Drug Administration, medication is deemed ineffective in 40%-70% of patients with common diseases. 1 Genome-wide analyses down to the single-cell level indicate that this limited responsiveness depends on both complexity and heterogeneity. Clinical studies have shown that predicting treatment response based on omics data from IMID patients is challenging. 2 Each disease can involve thousands of genes across multiple cell types, which vary between patients with the same diagnosis, and even between the same patient at different time points. 3,4 The clinical manifestations of IMIDs suggest an added layer of heterogeneity, namely, variable organ involvement in the same disease. IMIDs encompass more than 80 diseases, which include rheumatoid arthritis (RA), ulcerative colitis (UC), Crohn disease (CD), psoriasis (PSO), systemic lupus erythematosus (SLE), and many others. 5 As an example of variable organ involvement, RA can affect not only joints but also the skin and many internal organs, including the kidney, heart, and spleen. Successful pharmacological treatment of such variable organome-wide disease manifestations would ideally require answering questions such as: How many organs are affected? How complex and heterogeneous are the underlying molecular changes? Can those changes be organized into an overriding structure, which permits systematic, and increasingly detailed, analysis? Is there a hierarchy in the structure? Can that hierarchy be exploited to prioritize diagnostic and therapeutic targets?
These questions have not been investigated on combined organome-, cellulome-, and genome-wide scales. This type of analysis would involve challenges that are close to, or beyond, the limits of current understanding of the architectural principles of diseaseassociated changes in the multidimensional genome, for example: 1. Characterization of genome-, cellulome-, and organomewide changes: This can be achieved by single-cell RNA sequencing (scRNA-seq), which allows creation of atlases of all cell types in all organs in healthy mice and humans. 6,7 One reason as to why no similar effort has been made in disease states is that many organs are difficult or impossible to investigate in living human patients. Another reason is that internal organs may not give rise to specific symptoms. As an example, pathogenic mechanisms in lung have been proposed to have a primary role in RA, 8 but clinical and research foci are on joints. It is thus possible that important disease mechanisms, biomarkers, and drug targets are missed. Another problem with focusing on only one organ is that all organs may interact through the hematological, lymphatic, or nervous systems. Thus, they should ideally be studied together, rather than as individual parts. 2. Organization of organome-wide scRNA-seq data: We and others have previously described methods to organize scRNA-seq data from individual organs into multicellular disease models (MCDMs). 9,10 MCDMs are network models that show directed molecular interactions between cell types based on differentially expressed genes (DEGs) in each cell type and their predicted upstream regulators (URs) in other cell types. However, organization of MCDMs on an organome-wide scale is an unresolved challenge. 3. Prioritization of regulatory mechanisms in organome-wide scRNA-seq data: Because MCDMs have not been characterized on an organome-wide scale, their potential interactions have not been systematically investigated, nor is there any form of molecular or cellular hierarchy between organs. However, in a previous study, we found that interactions in an MCDM from one inflamed organ were multidirectional, without any evident hierarchy. 10 This complicated prioritization of diagnostic and therapeutic targets, which emphasizes the need to search for overriding structure to systematically find and prioritize regulatory mechanisms.
Here, we performed multi-organ scRNA-seq of a mouse model of collagen-induced arthritis (CIA) to develop a systemslevel strategy to define such structures, which could be validated in human IMIDs (Figure 1). Although disease mechanisms may differ between mouse models and human diseases, we reasoned that overriding structures would be comparable. In summary, we found complex and heterogeneous organomewide changes in CIA. 11 Those changes could be organized into a multi-organ MCDM (MO-MCDM) in which all organs interacted without evident hierarchy. However, despite the widespread molecular changes across all organs, only joints showed signs of inflammation. This contrast led to the question whether there could be an overriding structure in which complex mechanisms are required not only to activate but also to inhibit inflammation. If so, could that structure be systematically analyzed to prioritize such mechanisms and their URs? Combined analyses of multiorgan data from the mouse model and 10 human IMIDs supported that shared transcriptional programs were switched on or off by variable combinations of URs. Subsequent analyses of IMID patients who did or did not respond to treatment with anti-TNF (tumor necrosis factor), as well as more than 600 blood samples from SLE patients, supported that variable combinations of URs have the potential for personalized diagnostics and therapeutics in IMIDs. We propose prospective clinical studies to examine this potential and have made the data and methods freely available for such studies.

RESULTS
scRNA-seq shows highly diverse cellulome-and genome-wide expression changes in joints and multiple other organs in a mouse model of arthritis In order to search for systems-level principles to organize and prioritize disease mechanisms on organome-, cellulome-, and genome-wide scales, we performed Seq-Well-based massively parallel scRNA-seq of six DBA1/J mice with CIA and four healthy control mice. Three of the CIA mice developed mild arthritis (per limb arthritis score: 1-3) and three severe arthritis (per limb arthritis score: 4). We first analyzed 10 different organs, namely, joint, blood, draining lymph nodes, lung, thymus, skin, limb muscle, spleen, liver, and kidney, from at least one mouse with severe arthritis and one healthy mouse (Data S1). Despite only joints showing macroscopic signs of disease, we found DEGs between sick and healthy mice in multiple organs. The highest numbers of DEGs were found in muscle, joint, lung, skin, and spleen ( Figure S1A). We proceeded to analyze these five organs from all sick and control mice. After filtering and quality control, we recovered 2,230, 814, 4,565, 1,167, and 3,320 cells from joint, lung, muscle, skin, and spleen, respectively (see ''method details'' in STAR Methods; Data S1). Clustering and cell-type annotation revealed 13 cell types, namely, B cells, dendritic cells, endothelial cells, erythrocytes, fibroblasts, granulocytes, macrophages, monocytes, natural killer (NK) cells, T cells, myocytes, basal III cells, and neutrophils (see ''method details'' in STAR Methods; Figure 2A; Data S1). Cell-type proportions and DEGs varied greatly between organs (Figures 2B and 2C; Data S1). Of the total number of DEGs identified in macrophages and T cells, which were the only cell types identified in all five organs, 5% and 4%, respectively, intersected over all organs (Data S2). Pathway analysis of the cell types in the different organs resulted in a total of 501 pathways being significantly enriched in at least one cell type, although they were variably upregulated/downregulated in the different organs and cell types in which the direction could be inferred ( Figure S2A; Data S3).
The daunting complexity and heterogeneity of the molecular changes across multiple organs and cell types highlighted the overarching questions behind this study: how can disease-associated changes on organome-, cellulome-, and genome-wide scales be organized and prioritized? We reasoned that one A B Figure 1. Overview of the study (A1) Single-cell RNA sequencing (scRNA-seq) of a mouse model of collagen-induced arthritis (CIA) showed thousands of differentially expressed genes (DEGs) across all organs despite only joints showing signs of inflammation. (A2) Multicellular disease models (MCDMs) were constructed based on scRNA-seq data from the organs with the most DEGs. (A3) Transcriptional programs were identified in joints and muscle. These were turned on or off by partially shared combinations of upstream regulators (URs). (B1) Meta-analysis of multiple immune-mediated inflammatory diseases (IMIDs) showed a similar on/off switch that was (B2) regulated by different UR combinations in different IMIDs. (B3) Those URs have potential for personalized diagnostics and therapeutics, either using single-drug or combinatorial drug treatments. Histological analysis shows signs of inflammation in joints, but not in other organs To investigate the inflammation status in multiple organs and cell types, we conducted microscopic analyses of different organs from independent mice with severe CIA (clinical scores >8) and control mice. The results showed inflammation only in joints: significant infiltration of leukocytes in cartilage and synovium, together with bone destruction and synovial hyperplasia ( Figures 2D and S1B).
MCDMs show multi-directional networks in each organ without evident hierarchies The presence of macroscopic and microscopic signs of inflammation only in joints, despite cellulome-and genome-wide  (C) Heatmap presenting the similarity of DEGs. Rows and columns represent different cell types in respective organs, and the color scale corresponds to the Jaccard index. (D) Representative H&E images of the joint and muscle from control and CIA mice shown at a magnification of 1003 (scale bars, 100 mm). Red, black, and blue arrows indicate synovial hyperplasia, bone destruction, and synovial infiltration of inflammatory cells, respectively. BM, bone marrow; C, cartilage; S, synovial cavity. changes in all analyzed organs, suggested an overriding structure underlying the organization and prioritization of DEGs in this system: DEGs in joints activated inflammation, whereas DEGs in other organs suppressed inflammation. If so, this structure would hypothetically act as an on/off switch for inflammation, which could help to prioritize URs and downstream target genes for that switch. To test this hypothesis, we constructed MCDMs of each organ. The MCDMs described predicted molecular interactions between cell types in each organ. The interactions were bioinformatically inferred by linking the DEGs in each cell type with their predicted UR 12 (Data S4). DEGs linked with URs were referred to as downstream targets. 13 Because the interactions were directed, they could potentially be traced to prioritize an UR and cell type with a hierarchically superior role, as well as its downstream target genes in other cell types. We began by analyzing whether the joint MCDM had such a pro-inflammatory UR with an ''on'' role for the switch. The joint MCDM included URs of known pathogenic importance for both mouse CIA and human RA, such as Il1b and Tnf. 13 However, the MCDM showed multi-directional interactions mediated by many other URs without evident hierarchy ( Figures 3A and 3B). A similar, multi-directional organization was found in MCDMs from lung, spleen, muscle, and skin ( Figures 3C-3F). The complex, heterogeneous, and apparently non-hierarchical changes across multiple organs and cell types led us to attempt to prioritize URs on a multi-organ scale.
Ranking indicates that altered balance between pro-inflammatory and anti-inflammatory URs has an on/off function for inflammation To prioritize URs, we ranked the URs based on the size of their predicted molecular and cellular effects across all analyzed organs (see ''method details'' in STAR Methods; Figure 3G). The disease relevance of the ranking system was supported by four of the top-ranking URs in joints, Il1b, Tnf, Dsc3, and Ltf, being either therapeutically, functionally, or genetically associated with RA and other IMIDs. [13][14][15][16] In further support of the clinical relevance of those URs, the expression of their downstream targets differed significantly between mild and severe arthritis in several cell types in joints: Tnf in fibroblasts (p = 2.00 3 10 À3 ), Il1b in T cells (p = 4.14 3 10 À5 ), Ltf in T cells (p = 9.51 3 10 À5 ) and macrophages (p = 1.57 3 10 À6 ), Dsc3 in macrophages (p = 9.9 3 10 À7 ), and Mpz in macrophages (p = 1.84 3 10 À14 ) ( Figure S3; Data S4). Although the increased activity of pro-inflammatory URs in joints was expected, another finding was not: both Il1b and/or Tnf were also differentially expressed and predicted URs in organs that did not show macroscopic and microscopic signs of inflammation. However, in contrast with joints, the expressions of those URs varied and were counterbalanced by anti-inflammatory URs. For example, in muscle, Il1b increased, whereas Tnf decreased, and the anti-inflammatory UR Tgfb increased ( Figure S1C). By contrast, both Il1b and Tnf, but not Tgfb, increased in joints. Thus, the altered bal-ance between pro-and anti-inflammatory URs could act as an on/off switch for inflammation. Such a switch would explain why only joints showed signs of inflammation, despite organome-wide expression changes. Another unexpected finding was that although Tnf was downregulated in muscle, its predicted downstream targets in the same organ were activated in monocytes (p = 6.9 3 10 À4 ; Z score = 0.69), and in T cells (p = 1.02 3 10 À3 ; Z score = 1.182). A potential explanation was that the downstream targets were regulated by TNF derived from other organs and transported via the blood. If so, inflammatory mechanisms in different organs are interconnected. We hypothesized that this concept could be developed by searching for molecular interactions between MCDMs, such that a MO-MCDM would be formed.

MO-MCDMs connect inflammatory mechanisms in different organs into a multi-directional network
To investigate systematically if molecular interactions between MCDMs in different organs could be organized into an MO-MCDM, we used the same methods as for individual MCDMs. However, we included only URs predicted to be released into the blood (based on the Human Protein Atlas). We identified 1,966 of such inter-organ interactions, which were mediated by 48 URs (Figure 4; see ''method details'' in STAR Methods). The resulting MO-MCDM formed a multi-directional network in which all MCDMs were interconnected.
To validate that URs could mediate interconnectivity in the MO-MCDM, we performed protein analyses of high-ranking URs and interacting cytokines in sera from independent CIA mice ( Figure S4). Because CIA may variably involve different organs, these analyses were performed at different time points during disease progression. 17 In support of interconnectivity, all analyzed URs and interacting cytokines were found in sera ( Figure S4). TNF increased at early time points, whereas it decreased to normal expression at later stages of the disease. For interleukin (IL)-1b, no systemic changes in protein expression level could be seen at different stages of disease progression. However, for IL-1a, a significant drop in protein expression level was seen at later stages of the disease. IL-6 and interferon g (IFN-g), which are known to interact with IL-1b and TNF, 18-21 also showed variable changes in expression level at different time points of disease. Such variations could be consistent with dynamic changes in organ inflammation in CIA. 22 We next analyzed whether the altered balance between pro-and anti-inflammatory URs would be associated with an altered balance between downstream proand anti-inflammatory pathways.
Connective pathway analysis supported a graded switch system in CIA To systematically test whether the altered balance between downstream pro-and anti-inflammatory mechanisms explained why joints, and not muscle, showed signs of inflammation, we performed pathway analysis of all DEGs in different cell types from these two organs. In total, we identified 428 significantly enriched pathways in at least one cell type (Data S3). The large number of pathways complicated systematic testing of our hypothesis. A potential solution was suggested by the fact that 64% of all genes in the 428 pathways were shared by more than one pathway. This led us to hypothesize that a higher-order structure than pathways could be identified, namely, groups of pathways with partially shared genes (henceforth referred to as programs). If such programs were relevant for pathogenesis, they should (1) be enriched in genome-wide association study (GWAS) genes from human RA and (2) differ in activation profiles between joints and muscle. To find such programs, we developed a method called connective pathway analysis. This approach used the 1-Jaccard index as a distance metric for clustering of pathways (pathways that share many genes would then be more proximate to each other than those that do not and, therefore, be closer in the dendrogram; Figure 5A; see ''method details'' in STAR Methods).
We reasoned that cutting the dendrogram at different levels would provide a systematic approach to prioritizing the programs that differed most in activation profiles between joints and muscle and therefore would be most relevant for an on/off switch. At the highest level of the dendrogram, we found two main CIA-associated programs (CIA_P), CIA_P1 and CIA_P2. Both CIA_P1 and CIA_P2 were enriched for GWAS genes of human RA (p < 0.006; Data S5). We also tested whether similar programs would be found using Kyoto Encyclopedia of Genes   and Genomes (KEGG) instead of Ingenuity Pathway Analysis (IPA), and we found significant overlap (p < 0.0001; Figure S2B). Analyses of these programs did not support the existence of a discrete on/off switch: both programs included activated pathways in the non-inflamed muscle, compared with healthy control mice. This finding suggested a graded switch system in which the non-inflamed state was an intermediate in a continuous spectrum, in which healthy and inflamed organs were extremes. We next focused on CIA_P1 because pathways were mainly activated in joints and inhibited in muscle ( Figure 5B). To facilitate the identification of pathways that had the most opposing activation directions in CIA_P1 (i.e., being activated in one organ and inhibited in the other or being significantly enriched in one organ and not significant in the other; see ''method details'' in STAR Methods), we cut the dendrogram into 10 subprograms (CIA_SPs) ( Figure 5C). Of these, CIA_SP1.3, CIA_SP_1.1, and CIA_SP1.6 showed the highest percentages of pathways with opposing activation directions (79%, 61%, and 60%, respectively; Data S3). The highest GWAS enrichment was found in CIA_SP1.6 (Data S5). Further analysis of CIA_SP1.6 showed that 75% of its pathways were related to human RA (Data S3). CIA_SP1.6 contained pro-inflammatory pathways such as ''Acute Phase Response Signaling'' and ''IL-6 signaling,'' as well as anti-inflammatory pathways such as ''PPAR signaling'' ( Figure 5D). The pro-inflammatory pathways were mainly activated in joint and inhibited in muscle, whereas the anti-inflammatory pathways showed the opposite pattern ( Figure 5E). Il1b and Tnf were predicted URs of CIA_SP1.6 in both joint and muscle (Data S6). However, in muscle, the downregulation of the pro-inflammatory URs Tnf and Apoe 23 and upregulation of the anti-inflammatory Tgfb1 ( Figure S1C) could explain why no signs of inflammation were found. By contrast, increased expression of Il1b in muscle was consistent with partial activation of some pro-inflammatory pathways. For example, the ''Senescence Pathway'' showed mixed activation or inhibition in different cell types in both joint and muscle ( Figure 5D). This pathway has previously been implicated in IL-1-and TNF-induced tissue damage in RA. 24,25 The mixed activation pattern of the ''Senescence Pathway'' in joints and muscle is thus consistent with mixed activation of IL-1 and TNF in these two organs. Taken together, these findings support a graded, rather than discrete, on/off system. We next examined whether such a system could be translated to human IMIDs.

Meta-analysis of human IMIDs supported a graded switch system
To test the disease relevance of the graded switch system, we performed meta-analysis of 10 different IMIDs (Data S1): RA, UC, CD, PSO, Sjö gren's syndrome (SS), systemic sclerosis (SSc), atopic dermatitis (AD), juvenile myositis (JM), ''at risk for'' type 1 diabetes (T1D), as well as SLE. The SLE datasets included discoid lupus erythematosus (DLE), subacute cutaneous lupus erythematosus (SCLE), and lupus nephritis (LN). The meta-analysis was based on 32 bulk profiling datasets from human organ biopsies. The IMID biopsies were taken from inflamed and/or non-inflamed sites and compared with biopsies from healthy controls. Meta-analysis of DEGs from each IMID showed highly complex changes in inflamed and non-inflamed sites with 647 pathways that differed significantly compared with controls (Data S3). Similar to CIA, connective pathway analysis revealed two IMID-associated programs (IMID_P): IMID_P1 and IMID_P2 ( Figure 6A).
The disease relevance of both IMID_P1 and IMID_P2 was supported by significant enrichment for GWAS genes (Data S5). However, when we compared the pathway overlap between IMID_P1 and IMID_P2 with the corresponding programs from each individual IMID, we mainly found significant overlaps within IMID_P1 ( Figure 7A; see ''method details'' in STAR Methods). This observation indicated that IMID_P1 contained pathways that were shared across IMIDs, whereas IMID_P2 had more disease-specific pathways. In IMID_P1, 32% of pathways were activated in inflamed organ sites and inhibited, or not significant, in non-inflamed sites. The corresponding figure for IMID_P2 was 21%. Both IMID_P1 and IMID_P2 also included many activated pathways in non-inflamed sites. In agreement with a graded switch system across a spectrum between health and an inflamed phenotype, this more graduated response could predict an increased risk of a gradual shift in the balance between proand anti-inflammatory pathways toward an inflamed phenotype.
We next cut the dendrograms of both programs into subprograms (IMID_SPs) in order to prioritize the subprogram that had the most pronounced on/off pattern and GWAS enrichment. This analysis led to prioritization of IMID_SP1.6, which had 49% of pathways with opposing patterns between inflamed and noninflamed organ sites, and GWAS enrichment in 72% of the IMID datasets, with median (range) odds ratio (OR) among the significant = 5.34 ($3.54-8.94) (Figures 6, 7A, and 7B; Data S3 and Figure 5. Connective pathway analysis to systematically define and prioritize transcriptional programs in joint and muscle from the CIA mouse model (A) Outline of connective pathway analysis: (1) identification of genes that belong to a pathway; (2) mapping DEGs on that pathway; (3) pairwise comparison of pathway-associated DEGs; (4) assessment of pathway-associated DEGs overlap (Jaccard Index); (5) examples of two extreme situations, pathways X1 (green) and X2 (orange) have a high overlap of DEGs but not X1 (green) and X3 (purple); (6) hierarchical clustering using the 1-Jaccard index as distance matrix; and (7) dendrogram transformation into a tree-like structure. (B) Connective pathway analysis identified two main programs, CIA_P1 (purple) and CIA_P2 (green). Each pathway was labeled as ''activated'' (red), ''inhibited'' (blue), ''unknown activation'' (gray), or ''not significant'' (black). (C) Tree-like representation of CIA_P1 with subprograms indicated with different node colors. Color bars indicate main program CIA_P1 (purple) and subprogram CIA_SP1.6 (ochre). Node size represents the total number of cell types in which the pathway is significantly enriched. (D) Detail of CIA_SP1.6. Each node represents a pathway, and pie charts within nodes represent ratios of cell types in muscle and in joint (left and right part of the pie chart, respectively) for which the pathway was significantly enriched. Colors represent pathway activation profile. (E) Detail of selected pathways from CIA_SP1.6. Left (right) part of polar charts presents all cells in muscle (joint). Each pie sector represents one cell type. The degree to which the sector is filled with color represents enrichment Àlog10(p value), whereas the color shows pathway activation (blue for inhibition, red for activation, and gray for unknown activation).
The top-ranking pathways in IMID_SP1.6, in inflamed sites, were ''Acute phase response signaling,'' ''B cell receptor signaling,'' ''Chemokine signaling,'' and ''IL-6 signaling.'' By contrast, anti-inflammatory pathways, such as ''PPAR signaling'' and ''PPARa/RXRa Activation,'' 26-28 were inhibited ( Figure 6C). However, IMID_SP1.6 also included activated pro-inflammatory pathways in non-inflamed organ sites: ''Leukocyte extravasation,'' ''Natural Killer Signaling,'' and ''MS-RON Signaling,'' all of which can contribute to chronic inflammation, and thereby a switch from off to on. 29 Analysis of the URs of IMID_SP1.6 agreed with the graded switch system being regulated by variable combinations of pro-and anti-inflammatory URs. This could have important basic and clinical implications, namely, that the stepwise characterization of programs and subprograms, as described above, could help to prioritize, diagnose, and treat optimal combinations of URs on the levels of IMIDs, subgroups, and individual patients.
Combinatorial regulation of the graded switch system has diagnostic and therapeutic implications For all IMIDs combined, we found a total of 389 predicted URs (Data S6). Specifically, for each disease we found a median    Fas, IFN-g, IL-1a, IL-1b,  TLR3, and TNF. In agreement with the graded switch system depending on altered balance between pro-and anti-inflammatory URs, Fas, IFN-g, IL-1a, IL-1b, TLR-3, and TNF are mainly pro-inflammatory, whereas AR and ES-b are anti-inflammatory. [30][31][32][33][34][35] Unexpectedly, however, the predicted effects of these URs, based on Z score, contrasted with their measured fold changes (FCs; Figures S5A and 5B 5B; Data S6). These predictions agree with the clinical experience that anti-TNF treatment is effective in the two former diseases, but not in the latter ( Figure 7C). This led us to examine the effects of anti-TNF treatment on subprograms in UC and CD. The pathways analyses of DEGs after treatment between CD patients who responded to anti-TNF (GEO: GSE52746, 10 treated anti-TNF responders versus 7 untreated patients) showed significant enrichment of anti-TNF-targeted pathways among IMID_P1 and subprograms IMID_SP1.6, IMID_SP1.9, and IMID_SP2.5 (false discovery rate [FDR], <4.58 3 10 À2 ) (Data S6). The corresponding analyses of DEGs from patients with UC (GEO: GSE92415; 29 treated anti-TNF responders versus 32 untreated patients) showed significant enrichment of anti-TNF-targeted pathways in the subprograms IMID_SP1.2 and IMID_SP2.7 (FDR, <1.99 3 10 À2 ). The enriched subprograms are henceforth referred to as affected subprograms and the others as non-affected. The relevance of the affected subprograms was supported by TNF being predicted to regulate all of them in inflamed organs (Data S6). For subprogram IMID_SP1.6, which we in this study identified as highly relevant for IMIDs, six of seven pathways (where a |Z score| > 0 could be inferred) showed the opposite direction of activation as a result of treatment response compared with how they were affected by the disease (Figure S5C), indicating effective treatment response. We next tested the hypothesis that the non-affected subprograms could be explained by URs whose downstream targets overlapped with TNF, having redundant, ''backup,'' functions. In UC, we found 14 URs predicted to co-regulate non-affected subprograms. For example, NR4A2 was predicted to co-regulate 8 out of 16 non-affected subprograms ( Figure S6A). In support of this prediction, NR4A2 was significantly upregulated in inflamed samples of UC (FDR, 2.07 3 10 À4 ; logFC = 1.31). By contrast, affected subprograms were not predicted to be regulated by NR4A2. 31 We further tested the potential of other URs to take over the effect of TNF in UC patients who did not respond to anti-TNF treatment. The pathways enriched among the DEGs of re-sponders, respectively non-responders versus controls, showed similar main and subprogram associations ( Figure S6B). We identified 92 predicted URs for the DEGs of non-responders versus controls. By predicting the potential for each alternative UR predicted for non-responders to take over the downstream effect of TNF (see ''method details'' in STAR Methods), we identified TLR6 as a potential UR to take over the effect of TNF among the non-responders after treatment.
To verify that the downstream genes of TLR6 and TNF followed the expected FCs before and after treatment, we additionally analyzed the DEGs for responders and non-responders after treatment versus controls. As expected, we found a tendency for a higher |FC| of the TLR6 downstream genes among the non-responders compared with the responders, before and after treatment ( Figure S6C). However, only CXCL8 showed a significant difference between responders and non-responders after treatment (p = 3.62 3 10 À2 ) ( Figure S6C). Furthermore, the TNF downstream genes showed smaller differences compared with control for the treated responders, compared with any of the untreated groups or treated non-responders ( Figure S6D).
Because UC and CD primarily affect the intestine, we next analyzed whether different UR combinations were associated with variable organ involvement in another IMID that often shows multi-organ involvement, namely, SLE. The clinical relevance lies in that this could indicate the need for diagnostic and therapeutic targeting of different URs in patients with different forms of organ involvement.
Different combinations of UR proteins in sera were associated with different subtypes of SLE, as well as with disease severity We analyzed 18 predicted UR proteins of IMID_SP1.6 in sera from two clinical visits of 304 patients who had been prospectively seen by the same rheumatologist (C.S.) at the tertiary referral unit, Linkö ping University Hospital, according to standardized criteria. 36 The American College of Rheumatology (ACR) criteria were used to define disease phenotype. 37 We constructed regression models to estimate relationships of UR proteins with ACR criteria, as well as measures of disease activity and organ damage (SLE disease activity index-2000 [SLEDAI] and Systemic Lupus International Collaborating Clinics/ACR index [SDI], respectively) ( Figure S7). Clinical variables, including results from physical examinations, laboratory values, age, gender, and duration of the disease, were collected and used in all our models. We also included the treatment information in regression models that predicted SLEDAI and SDI. We did not use treatment information in models predicting ACR criteria, because the treatment is dependent on the patient phenotype, e.g., patients with LN are widely treated by mycophenolate mofetil (MMF). We showed that different UR proteins were associated with different ACR phenotypes (Data S1). For example, cases classified with ACR-1 (malar rash) were positively associated with CD40-L and negatively associated with Fas. Cases classified with ACR-7 (renal disorder/LN) were positively associated with FAS but negatively associated with hepatocyte growth factor (HGF). TNF was positively associated with ACR-9 (hematological disorder) and ACR-10 (immunological disorder). By contrast with ACR-7 and ACR-9 (hematological disorder), TNF and HGF were significantly associated with ACR-6 (serositis), TNF negatively and HGF positively.
Next, we examined whether different combinations of potential URs were associated with disease activity in patients with (n = 80) and without (n = 224) LN. The regression models showed that different proteins were associated with the disease activity in the two groups (Data S1). IL-1a, IL-4, Fas, and oncostatin M (OSM) were associated with SLEDAI in patients with LN, whereas TNF, HGF, and CD40 were associated with SLEDAI in patients without LN. In agreement with our hypothesis, the URs that were negatively correlated with SLEDAI may have anti-inflammatory roles. However, testing whether altered expression of these URs inhibit or activate inflammation is complicated by their context-dependent pleiotropic roles. For example, OSM has an anti-inflammatory role on synovial cells: it reduces IL-1 and TNF expression, but it has a pro-inflammatory role on endothelial cells by inducing leukocyte recruitment and IL-6 production from endothelial cells. 38,39 A similar pleiotropy has been described for HGF and IL-4. Focusing on organ damage (Data S1), we found that TNF, IL-27, OSM, and TGF-b1 were associated with SDI in patients with LN. In contrast, IL-6, IL-1a, IL-2, HGF, and CD40-L were associated with SDI in patients without LN.

DISCUSSION
The main problem behind this study is that many patients with IMIDs do not respond adequately to treatment. 1 An important reason for this inadequate response is the daunting complexity and heterogeneity of the molecular changes in these diseases. scRNA-seq studies of IMIDs have shown altered expression of thousands of genes across multiple cell types in individual organs with phenotypic signs of disease. [40][41][42] Despite detailed information about those changes, translation of the data to personalized medicine has proven difficult. Thus, there is a wide gap between the complexity of disease-associated changes and health care today. Bridging that gap will likely involve great challenges for researchers and clinicians.
Variable multi-organ involvement in each IMID adds to the complexity and heterogeneity. This indicates the need for characterization, organization, and prioritization of molecular changes on organome-, cellulome-, and genome-wide scales. Our multi-organ scRNA-seq analyses of a mouse model of an IMID, namely, CIA, showed extensive changes on all those scales. Although those changes could be organized into an MO-MCDM, this analytical approach showed no evident molecular or cellular hierarchy that would allow prioritization of molecular changes. An unexpected finding led us to a potential solution: despite the organome-wide changes, only joints showed signs of disease. That contrast led us to hypothesize that the expression changes were organized into an overarching structure designed to switch inflammation on or off. We developed an analytical strategy that supported such an on/off switch and showed that it depended on altered balance between proand anti-inflammatory URs. Such a switch has been previously suggested in inflammatory responses and validated by functional studies of individual genes and cell types. [40][41][42] Given the complexity and heterogeneity of the organome-, cellulome-, and genome-wide changes, ranking of URs is crucial for understanding and prioritization of disease mechanisms. Our strategy provided a solution for systematic characterization and prioritization of URs, as well as their downstream target genes. We found that URs could be ranked based on the size of their effects on the downstream genes. In support of clinical relevance, the top-ranking URs included known therapeutic targets in IMIDs, including IL-1 and TNF. The downstream target genes could be organized into two main programs, and their subprograms, which permitted increasingly detailed analyses of pathways. In support of disease relevance, both programs were significantly enriched for genes identified by GWAS of human RA. However, large molecular changes in non-inflamed organs, including partially activated pathways, did not support a discrete on/off switch, in which non-inflamed organs corresponded to an ''off'' state. Instead, the activated pathways could increase the risk of a switch to an ''on'' state. This observation is consistent with a graded on/off switch, in which non-inflamed organs are intermediates on a spectrum, where healthy and inflamed organs represent extremes. Such an intermediate ''risk'' state could explain an important characteristic of both CIA and human IMIDs, namely, variable organ involvement during disease progression. A graded on/off switch has been previously described in model organisms and proposed to be generally applicable to biological systems. The relevance of molecular gradients in disease is supported by previous findings that variable expression or dysregulation of interconnected genes will define whether an organ is affected by disease. 43 The translational relevance of the strategy was supported by meta-analyses of human IMIDs. These showed a similar organization as in CIA, with URs, programs, and subprograms that agreed with a graded on/off switch. High-ranking URs included known drug targets such as IL-1 and TNF. However, except for a core group of URs, these varied across diseases and organs. For example, TNF was a predicted UR in IBD, but not in SS, which agrees with the clinical experience that anti-TNF treatment is effective in the former, but not in the latter. Clinical implications may be that characterization and ranking of URs will be needed for successful treatment, on the levels of IMIDs, subgroups, or even individuals. Those implications were supported by our analyses of the effects of anti-TNF treatment in patients with IBD who did or did not respond to that treatment. We found that lack of response could be explained by overlapping downstream effects of the URs, such that the effect of inhibiting one UR could be diminished by one or more other, functionally redundant, ''backup'' URs. Prospective clinical studies are warranted to test whether predictive classifiers for treatment response can be developed based on high-ranking URs. To examine whether UR combinations would vary in patients with different forms of organ involvement, we focused on an IMID with highly variable multi-organ manifestations, SLE. Analyses of predicted UR proteins in more than 600 sera from SLE patients did show that different combinations of UR proteins correlated with different forms of organ involvement. Furthermore, different combinations of URs were associated with disease severity in SLE patients with and without renal involvement. A clinical implication may be that different URs should be targeted in these two subgroups. Moreover, some shared URs showed opposing associations. Interestingly, we found that HGF was positively associated with the damage index but negatively associated with organ involvement in SLE patients without renal involvement. This supports our message that a systems-level strategy for prioritization of URs on an organome-wide scale is important. Similar to CIA, several subprograms in non-inflamed organs included activated pro-inflammatory pathways. The pathogenic relevance of those subprograms was supported by enrichment of genes identified by GWAS. In agreement with a graded switch system, these subprograms could enhance risk of altered balance between proand anti-inflammatory pathways, resulting in an inflammatory phenotype or amplification of that phenotype. A potential clinical implication is the development of combinatorial diagnostic and therapeutic targeting of such URs during disease remission to prevent a gradual switch to active disease. 3,43 We propose that our strategy to organize and prioritize disease-associated changes on organome-, cellulome-, and genome-wide scales has significant potential for future studies aimed at personalized combinatorial diagnostics and therapeutics. We have made the methods and data freely available for such studies.

Limitations of the study
The analytical strategy was derived from multi-organ analyses of a mouse model of CIA, which may not be representative of human disease. Moreover, the CIA model includes use of an adjuvant to enhance the inflammatory response, which could induce pro-inflammatory pathways in tissues other than joints. The construction of MCDMs and the predicted effects of combinations of URs on downstream genes were based on previously described or predicted protein interactions, which may be confounded by knowledge bias. Although the analyses of the CIA mouse model were based on scRNA-seq, the analyses of human IMIDs were performed on bulk RNA-seq. From a translational perspective, future experimental and clinical studies are warranted to examine the diagnostic and therapeutic implications of the study, in particular the effects of URs on downstream genes.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:  Human data Samples were obtained from 304 patients (263 women, 41 men) classified with SLE according to the 1982 American College of Rheumatology and/or the Fries's diagnostic principle 83 (Data S1). All subjects had provided oral and written informed consent. The study protocol was approved by the Regional Ethics Review Board in Linkö ping (M75-08/2008). All subjects were included in the prospective and observational research program Clinical Lupus Register in North-Eastern Gothia at the Rheumatology Unit, Linkö ping University Hospital48. Patients were not involved in the design, conduct, reporting or dissemination plans of our research. Serum was available from each patient at two different time-points from which disease activity had been assessed by the clinical SLE disease activity index (SLEDAI) and damage accrual by the Systemic Lupus International Collaborating Clinics/ACR damage index (SDI). 84,85 The recent treatment of the patients prescribed at the previous visit was included in the clinical information (Data S1).

Study design
Our aims were to characterize, organize and prioritize disease-associated organome-, cellulome-and genome-wide changes in IMIDs (Figure 1). We combined single cell and bulk multi-organ profiling of mouse and human IMIDs. We found complex and heterogeneous organome-wide changes in a mouse model of CIA, which could be organized into a multi-organ multicellular disease model (MO-MCDM). In this MO-MCDM all organs interacted, without evident hierarchy. Despite the organome-wide molecular changes only joints showed signs of inflammation. This contrast led to the identification of an overriding structure in which shared transcriptional programs were switched on or off by variable combinations of URs. Analyses of IMID patients who did or did not respond to treatment with anti-TNF, as well as more than 600 blood samples from SLE patients, supported a graded on/off switch regulated by variable combinations of URs, which have the potential for personalized diagnostics and therapeutics.

CIA mouse model generation
For scRNA-seq analysis, CIA was established following a previously described method. 11 Six mice were immunized with 100 mg (50 mL) bovine collagen II (BC-II, Chondrex, USA) emulsified with 50 mL Complete Freund's Adjuvant (CFA) (Sigma-Aldrich, USA) in 1:1 ratio, via intradermal injection near the base of the tail. A booster immunization was administered on day 20 with 100 mg of BC-II emulsion (prepared with 1:1 incomplete Freund's adjuvant (IFA)) and injected at the base of the tail. 100 mL of Phosphate-Buffered Saline (PBS) was injected similarly to control mice. The severity of arthritic limbs was scored on a 0-4 scale, 0: normal; 1: swelling and redness in one digit; 2: swelling and redness in more than one digit or swelling and redness in one digit, wrist and ankle; 3: Swelling and redness presenting in paw and digits; 4: maximum inflammation of limb involving all joints and digits as described in the protocol by Brand et al. 11 The arthritic score for each mouse was the sum of the scores of arthritic limbs. The mice were sacrificed when they achieved scores of 8-12 or after they had been immunized for 60 days under isoflurane anesthesia via cervical dislocation. The joint, blood, draining lymph nodes, lung, thymus, skin, limb muscle, spleen, liver, and kidney were collected for further analysis. For the independent histology analysis, four 8-week male DBA1/J mice were immunized intradermally in the proximal tail with 100 mL of emulsified chicken type II collagen (2 mg/mL, Chondrex, USA)/CFA (1 mg/mL). The clinical arthritis score was evaluated for each limb from 0 to 4 with a maximal score of 16 for each mouse. 86 The healthy control and severe CIA mice (clinical score >8) were sacrificed, and the knee joints, lungs, livers, kidneys, skin, and hindlimb muscles were collected after heart perfusion. 17,87 Histological analysis Whole knee joints, lungs, livers, kidneys, skin, and hindlimb muscles from healthy control and CIA mice were fixed in 4% formaldehyde. Joints were further decalcified with Decalcification Solution (ServiceBio, G1107, China) for 7 days. The specimens were then embedded in paraffin and sagittal sections (4 mm) were cut. The sections were stained with hematoxylin and eosin (H&E, Sigma-Aldrich, USA) for the histology analysis. Histological sections were assessed for infiltration of cells into the synovial cavity resulting in inflammation, proliferation of cells in the synovial layer, and bone erosion.

Sample cryopreservation
All dissected organ samples were placed into suitable tubes with freezing solution (10% DMSO and 90% FBS), placed into a CoolCell LX box (Corning, USA), and frozen with gradually decreasing temperature (1 C/min) to À80 C. The samples were then stored at À175 C until further analysis.

Sample thawing
Before digesting all the organs and harvesting the single cells, the cryopreserved samples were thawed following 88 with slight differences. Briefly, cryopreserved samples were quickly thawed in a 37 C water bath with continuous agitation, then transferred into 15 mL centrifuge tubes with 1 mL pre-warmed thaw solution (90% Hibernate-A and 10% FBS) and incubated at room temperature for 1 minute. Next, 2 mL, 5 mL, and 5 mL thaw solutions were added into the centrifuge tube, separated by 1-minute incubation. The samples were then centrifuged at 350 3 g for 5 minutes. Lastly, the samples were resuspended with 1 mL Hibernate-A, after which the supernatant was removed and incubated until the next step.

Single cell suspension
The thymus, spleen, and lymph node were thawed as described above, and were passed through a 70 mm strainer to collect cells in a 50 mL centrifuge tube. After being centrifuged at 350 3 g for 5 minutes, cells were resuspended with 5 mL red blood cell lysis buffer for 5 minutes. The lysis reaction was quenched by adding medium (90% RPMI-1640 with 10% FBS). Cells were centrifuged at 350 3 g for 5 minutes and washed thrice to remove the lysis buffer. Single cell suspensions were prepared with RPMI-1640 at a density of 1 3 10cells/mL. Samples from different organs (whole knee joint, muscle, lung, skin, liver, and kidney) were quickly transferred into 75 mm dishes with 1 mL DMEM after thawing, and then minced into $1mm pieces with scissors. Next, pieces of organ samples were transferred into 15 mL centrifuge tubes containing 5 mL DMEM. Different organ samples were treated with different enzymes for different durations (Data S1). After dissociation, another 5 mL DMEM with 10% FBS was added to the 15 mL centrifuge tubes. Dissociated cells were centrifuged at 350 3 g for 5 minutes after passing through a 70 mm strainer, after which the cells were washed thrice with PBS. Single cells were resuspended into RPMI-1640 at a density of 1 3 10 5 cells/mL for cell loading. Peripheral blood mononuclear cells (PBMCs) were isolated as previously described. 89 Briefly, 0.5 mL peripheral blood was diluted with an equal volume of PBS (calcium free), which was further loaded on the top of 1 mL Lymphoprep followed by centrifugation at 800 3 g for 30 minutes at room temperature in a swinging bucket rotor with the brake off. PBMCs were retrieved and washed with PBS. Erythroid cells of different organs and peripheral blood of each mouse were removed using RBC Lysis Buffer (Bio Legend, USA) (Data S1). Single cell suspensions were prepared by resuspension of PBMCs with RPMI-1640 at a density of 1 3 10 5 cells/mL. scRNA-seq wet-lab protocol All scRNA-seq experiments were performed using the Seq-Well technique. 90 Briefly, prepared single cell suspensions were coloaded with barcoded and functionalized oligo-dT beads (Chemgenes, USA; cat. no. MACOSKO-2011-10) on microwell arrays synthesized as described. 90 For each sample, 20,000 live cells were loaded onto an array to bind with oligo-dT beads. The arrays, covered with plasma-treated polycarbonate membranes, were placed in a 37 C incubator for 30 minutes. Next, beads were collected to perform cell lysis, hybridization, reverse transcription, and whole transcriptome amplification. Libraries were then prepared for each sample using the Nextera XT DNA Library Preparation Kit (Illumina, USA; cat. no. FC-131-1096) according to the manufacturer's instructions. Libraries from three samples were pooled together and sequenced using the NextSeq 500/550 system, and sequencing results were analyzed as described below.

Cytokine analyses in peripheral blood
Approximately 100 mL of blood were collected from ten healthy control DBA1/J mice, as well as CIA mice at week three (before symptom onset, n = 12), week eight (early stage after symptom onset, n = 12) and week 15 (later stage after symptom onset, n = 9) after CIA induction by retro-orbital bleeding. Twenty-five mL of serum were used for assaying inflammatory cytokines using LEGENDplex TM Mouse Inflammation Panel (13-plex) (CAT: 740446, BioLegend, USA), including IL-1a, IL-1b, IL-6, IL-10, IL-12p70, IL-17A, IL-23, IL-27, MCP-1, IFN-b, IFN-g, TNF, and GM-CSF. The assay was performed according to the manufacturer's protocol and the data were collected on a BD FACSAria III flow cytometer and analyzed by Flowjo. The mean fluorescence intensity of each cytokine of the standards was used for calculating the standard curve for each cytokine using a log-log curve fit. The difference in concentration of each cytokine between the different time points was calculated using the Wilcoxon rank sum test, as described above.

The analysis of the UR protein expressions in SLE sera
We used the clinical variables, patient information, drug treatment and protein levels (Data S1) to estimate the association between the protein levels and the patient phenotype, organ, damage, and disease activity. To preprocess the data, we log-transformed the protein levels, and used random forest imputation 79 to impute the missing values. We removed the observations where the response variable was not available, and we scaled the input variables to zero mean and unit variance. In case of treatment information, we only included drugs that were used to treat at least ten patients. As the distribution of the response variables differed, we fitted different regression models for each respective response variable. We used a logistic regression, 91 negative binomial regression 92 and zero inflated negative binomial regression 93 to predict the patient phenotype, organ damage, and disease activity respectively. We also  94 to perform variable selection. These models returned the coefficient estimates as well as a p value, which tested a hypothesis that the estimate is zero. A positive coefficient estimate reflects that the patients with high response variables had high protein levels, and we further referred to it as a positive association. In contrast, a negative coefficient estimate reflects that the patients with high response variables had low protein levels, and we further referred to it as a negative association. To assess whether the tested models follow their assumptions, we standardized the residuals between 0 and 1, then compared them against the assumed distribution and performed outlier and dispersion checks. We also plotted the standardized residuals against the rank transformed predicted variable where we expected a uniform distribution. 95 R functions missForest, 79 glm, 79 glm.nb, 96 were used to perform the analysis, and function simulateResiduals 95 was used to investigate the residuals.

Measurement of the UR proteins in SLE sera
We measured the protein levels of the 18 predicted UR proteins of most disease relevant IMID_SP in the serum of 304 patients at two separate phlebotomies using human magnetic multiplex beads assay (R&D, Bio-Techne, USA) according to the manufacturer's instructions. The five-parameter logistic curve was used to generate the standard curve. The Limit of Detection (LOD) was defined as the standard point with the lowest concentration of an analyte that can reliably distinguish signal from background noise. The Lower Limit of Quantification (LLOQ) is defined as the standard with the lowest concentration. The Upper Limit of Quantification (ULOQ) was defined as the standard point with the highest concentration. If protein values lay outside of the interval between LLOD and ULOQ but higher than LOD, we used extrapolated data for further analysis. We excluded IL-17 and GM-CSF from the analysis as 86% and 75% of measurements respectively were below LOD.

Organ prioritization
The relevance of the organs for CIA development was tested in a pilot study, in which at least one sample from each organ was sequenced (Data S1) and processed as described. 10 In short, the data from each sample were extracted, and poor-quality cells were sorted out as described below. Sequencing data from two of the organs, liver and kidney, did not meet quality criteria (%25 cells with 10,000 reads per cell) and were excluded from further analysis. The data from the remaining samples was knn-smoothed (k = 12), whereafter DEGs were identified between sick and healthy individuals for each organ separately using Monocle 70,97 as described in. 10 For comparative analysis between organs, 40 cells were bootstrapped from each group of healthy and sick individuals, for ten sampling rounds, for the differential expression analysis. The number of DEGs (Benjamini-Hochberg adjusted p value (FDR) < 0.05) for each sampling was compared and the five organs with the highest number of DEGs (joint, lung, muscle, skin, and spleen) were selected for downstream analysis (Data S1 and Figure S1A).

scRNA-seq data processing
The single cell data from the different mice samples were processed into digital gene expression matrices following James Nemesh, McCarrol's lab Drop-seq Core Computational Protocol (version 1.0.1, Drop-Seq tools v1.12) (http://mccarrolllab.com) using bcl2fastq (v2.19.1) conversion and Picard software (v2.9.0). The indexed reference for alignment of the reads was generated from GRCm38 (June 2017, Ensembl) using STAR software (v2.5.3). 72 Only primary alignments towards the reference genome were considered during downstream analyses, according to the mapping quality using STAR software. The quality of cells was assessed by having a minimum of 10,000 reads, 400 transcripts, 200 genes and less than 20% of mitochondrial genes per cell. The five organs with the highest number of DEGs based on the organ prioritization described above, namely, joint, lung, muscle, skin and spleen, were then analyzed together. Outliers were removed based on an overestimation of transcripts count (i.e., cells with more than 6,000 transcripts) due to the risk of duplicates in the library. For a gene to be included in the data, it needed to be identified in at least 10% of the cells.

Clustering and cell type identification
We used a reference-based approach to identify cell types. As a reference dataset, we used mouse bulk expression data of sorted cell populations available in the R package SingleR (v1.0.6). 73 To preprocess the data for cell type identification, we only retained 6,395 shared genes for both the bulk expression data and the single cell data. The resulting single cell data was denoised by the deep count autoencoder (DCA, v0.2.3) 74 with the default settings. This method has an in-built data normalization and outputs 1) a denoised dataset corrected for dropouts and varying library sizes, where each data value represents the expected (denoised) gene expression, and 2) a latent representation of the denoised data in a 32-dimensional latent space. Existing reference-based single cell type identification such as 73 or 98 uses correlation measures, either Spearman or Pearson, to match the single cell observations to the reference data. While Pearson correlation measures linear relationships, Spearman correlation is more general in the sense that it accounts for monotonic (i.e., non-linear) relationships. Accordingly, for each reference point, a single cell observation having the highest Spearman's correlation to that reference was found, and then a Monotonic Regression (MR) 99 was computed with the reference vector as the input variable and the natural logarithm of the denoised single cell expression of the selected cell as the output variable in the regression. The exponents of the predicted values from the MR were then treated as reference data expressions in the scale of the single cell data. This means that the rounded exponents of the fitted MR values for all reference data points were used as inputs in the DCA Cell Reports Medicine 4, 100956, March 21, 2023 e5 Article ll OPEN ACCESS that was estimated in the data validation step, and the latent representations outputted by the autoencoder were used in the following step. The latent space observations obtained from both the single cell data and from the bulk data were clustered together by the Leiden's algorithm. 100 In the Leiden's algorithm, we set the size of local neighborhood to 30 and the resolution parameter (which determines the number of clusters) was set in such a manner that 70% of cells in each resulting cluster matched the same cell type. Briefly, we started with a resolution parameter of 0.5, resulting in a low number of clusters. For each single cell observation in each computed cluster, a match was computed by finding a reference point that belongs to that cluster and has the highest Spearman correlation to the denoised gene expression of the single cell. The purity of a cluster was determined by computing the proportion of single cell observations within that cluster that had the same reference type. The resolution was increased until each cluster had a satisfactory purity, which resulted in resolution = 1. Finally, all single cell data points within a cluster were labeled with the reference label corresponding to the dominant reference match within that cluster (Data S1). All clusters recognized as the same cell type were merged for further analysis, resulting in 13 groups of cells (Figure 2A). The cell types identified, as well as those unidentified due to a lack of cell types in the reference, were further validated, or identified, using marker genes as described below.

Data normalization and DEG analysis
To calculate DEGs in the dataset, we used the single cell variational inference (scVI, v0.7.1) framework. 76 First, the variational inference model was set-up based on the UMI-count data, reducing batch effect based on the input samples, i.e., individual mice and organs, after which the model was trained using default parameters. The DEGs were then identified between CIA and healthy mice for each cell type in each organ separately, using the 'change'-mode. The significant DEGs were identified as those with 'is_de_fdr_0.05 = True' in the scVI differential expression output. To infer the direction of change, we used the mean log(fold change) ('lfc_mean') values produced, where a positive fold change (FC) indicated upregulation while a negative FC indicated downregulation in cells from CIA compared to healthy mice (Data S2). The data used for differential expression analysis were normalized by the scVI autoencoder, correcting for variation in sequencing depth. The normalized expression matrix and the latent space representation data were used for single cell downstream analyses as further described. To compare how much the lists of DEGs differed between cell types and organs, the Jaccard index was calculated for each pair of gene lists.
Cell type identification using marker genes Marker genes were used to validate the cell types identified by the Leiden algorithm, and to obtain the identity of clusters that were not represented by any reference data point. For all organs combined, the marker genes were calculated based on DEG analysis as described above. The marker genes were defined as those being significantly differentially expressed between each cell type or cluster of unknown identity (based on the previous identification) and all other cell types in the dataset (Data S1). We then searched for the known marker genes (Table S1) within the sets of cell type-specific marker genes. If marker genes of a certain cell type were enriched in a cluster, i.e., a positive FC, the cluster was identified accordingly.
For each cell type identified in any of the organs, the variation in cell type proportion over organs was calculated using ANOVA, adjusting the p values by Holm correction.

Cell-cell interaction analysis
The interactions between cells were identified by analyzing the data using NicheNet (v1.0.0). 12 NicheNet is an R package developed for identification of inter-cellular interactions based on lists of potential interacting genes and a database of known upstream regulators (URs) to target interactions. For these analyses, we used the lists of DEGs between CIA and healthy mice to find the interactions which change due to arthritis. As NicheNet requires human gene symbols as input, the mouse genes were translated into their human orthologs based on the Ensembl genes (http://www.ensembl.org/biomart/martview, GRCh38.p13, Downloaded June 6, 2020). These human orthologs were used for all downstream analyses based on these cell-cell interactions.
The cut-off to define the expressed genes in the data was set according to the author's recommendation, to give 5,000 to 10,000 expressed genes for the sender and receiver cell population independently. First, the previous transformation of the scVI-normalized expression data inverse logarithm was applied (10-1). Next, genes with a mean expression level R1 3 10 À5 in the population of cells were defined as expressed. The interactions were then identified, based on the lists of DEGs, between each pair of cell types using the default analysis set-up. For each interaction identified, all of the potential target genes in the source cell type were identified using the get_weighted_ligand_target_links() function and its default settings. To focus our analyses on the strongest interactions of the networks, we included only those with a Pearson Correlation Score (PCC) > 0, meaning that the target genes of the interactions are enriched among the differentially expressed genes.
When analyzing inter-organ interactions, the UR-target interactions identified by NicheNet were curated only to include those that are biologically feasible between organs. For this aim, we used Ingenuity Pathway Analysis (IPA, Q4 2020, Qiagen, Germany) to specify the cellular location known for each potential UR (Data S6). The list of inter-organ interactions was then curated to include only URs located in extracellular space, as they have the potential to be transported through the blood.

MCDM construction and UR prioritization
For each organ separately, a MCDM was constructed based on predicted molecular interactions between all cell types in the organ. The predictions were inferred using NicheNet 12 as described above. The MCDM thus consists of cell types as nodes and unique Article ll OPEN ACCESS interactions as lines. A unique interaction represents one cell type-UR-cell type combination, thus enabling multiple edges, based on different UR, between each pair of nodes. The MO-MCDM was created in the same way, but only based on inter-organ interactions as described above. Thus, the nodes represent cell types in each organ, and the lines represented interactions between cell types in different organ.
The URs were then prioritized based on their downstream effect. First, they were ranked based on the total number of predicted downstream target genes that they were predicted to regulate, in all cell types and organs combined. Secondly, they were ranked based on the number of downstream cell types and organs which they were predicted to target.
Pathway enrichment analysis Identification of pathways was performed with IPA. We used the core analysis in the IPA software to identify canonical pathways based on a list of DEGs (e.g., all the potential targets of all the URs in the different organs). When this analysis was completed for each cell type and organ separately, the Bayes factors from the differential expression analysis were included to define the direction of change due to arthritis.
IPA consists of a global network that is based on manual curation of a vast body of medical literature and biomedical databases, which is continuously updated. 101 The core analysis in IPA (parameters: species = mouse) was used to identify pathways that were significantly enriched among the list of genes. Statistical analysis was performed using Fisher's exact test, right tailed, within the IPA software (Q1 2021 and Q4 2020 version). 77 All pathways with p < 0.05 were considered significantly enriched. Pathway activation direction was indicated by IPA activation z-score as activated (z-score >0) or inhibited (Z score <0).
The association between UR-target expression and arthritis score The association between UR-target expression and arthritis score was calculated by pairwise comparison of the target mean expression levels between cells from mouse joints with mild arthritis (score 1-3), or severe arthritis (score 4) (Data S1), as well as healthy control mice, using a Wilcoxon rank sum test. The gene names of the normalized data from scVI analysis were first translated to their human orthologs as previously described. Thereafter, the data were standardized, producing a mean expression of zero and standard deviation of one, over all genes within each cell. The mean expression level of all UR target genes was calculated for each cell. The differences in target expression were then computed, in a pair-wise manner between the groups using the wilcox_test() and add_significance() functions in the R package rstatix (v0.7.0), producing p values adjusted for multiple testing by the Holm correction.

Connective pathway analysis
Pathway analysis was first conducted using IPA, as described above, on all DEGs for each dataset/cell type. To systematically assess similarities and differences between pathways in the different groups, (i.e. joint and muscle for CIA or inflamed and non-inflamed organ sites for IMIDs), we first clustered the pathways based on gene list similarities. To do so, the Jaccard-index was calculated, and values of 1-Jaccard-index used as distances for hierarchical clustering. Clustering was performed using the hclust() function in R with Ward's method (i.e. parameter ward.D2). Only genes differentially expressed in at least one cell type (CIA)/dataset (IMIDs) within each group (active and inactive disease organ) were considered for clustering. The dendrogram from the hierarchical clustering was next transformed into a tree-like structure, in which each node represents one pathway. To prioritize clusters for downstream analyses further, each pathway was labelled as ''activated,'' ''inhibited,'' ''no_direction,'' or ''not_significant'' in each group separately. Labelling was based on the ratio of datasets showing a specific direction of change (indicated by the IPA activation z-score). For example, in the CIA data, the ''Acute Phase Response Signaling'' pathway was found significantly enriched (overlap p < 0.05) in six cell types out of 12 in muscle: activated (z-score >0) in two cell types and inhibited (z-score <0) in six cell types. Since the pathway was inhibited in a higher number of datasets, ''Acute Phase Response Signaling'' was labelled as ''inactive'' in muscle. In cases in which there was an equal number of cell types with a pathway that was activated and inhibited or if the direction could not be predicted (z-score = NA), the pathway was labelled ''no_direction'', and if none of the cell types showed significant enrichment of the pathway, it was labelled ''not_significant''. This oversimplification helps to achieve an overview of general behavior of pathways in the distinct groups, and to prioritize clusters of pathways for further analyses.
For further validation, KEGG pathway analysis was also conducted on all DEGs for each cell type in muscle and joint using ''en-richKEGG'' function from R package ''clusterProfiler''. Then, connective pathway analysis was performed on the KEGG pathways with p < 0.05. First, similarities and differences were assessed between pathways using the Jaccard-index. Clustering was performed as previously described and divided into two programs. To assess whether connective pathway analysis using KEGG pathways has comparable results as using IPA, similarities between programs from connective pathway analysis of KEGG and IPA were compared using enrichment analysis (two-sided Fisher's exact test). Since all pathways in KEGG have different pathway names, pathways were not directly compared between KEGG and IPA; instead, the enrichment analyses were performed on gene content of each program using Fisher's exact test. All genes enriched in KEGG pathways were used as background for the analyses.
Prediction of URs regulating specific programs/subprograms in connective pathway analysis To identify potential URs of each program/subprogram from the connective pathway analysis, we performed enrichment analyses. For each potential UR-program/subprogram pair we computed the enrichment of the URs downstream targets among the pathwayassociated genes (right tail Fisher's exact test). All genes in connective pathway analysis were used as background. Enrichment