Spatial transcriptomic characterization of COVID-19 pneumonitis identifies immune circuits related to tissue injury

Severe lung damage in COVID-19 involves complex interactions between diverse populations of immune and stromal cells. In this study, we used a spatial transcriptomics approach to delineate the cells, pathways and genes present across the spectrum of histopathological damage in COVID-19 lung tissue. We applied correlation network-based approaches to deconvolve gene expression data from areas of interest within well preserved post-mortem lung samples from three patients. Despite substantial inter-patient heterogeneity we discovered evidence for a common immune cell signaling circuit in areas of severe tissue that involves crosstalk between cytotoxic lymphocytes and pro-inflammatory macrophages. Expression of IFNG by cytotoxic lymphocytes was associated with induction of chemokines including CXCL9, CXCL10 and CXCL11 which are known to promote the recruitment of CXCR3+ immune cells. The tumour necrosis factor (TNF) superfamily members BAFF (TNFSF13B) and TRAIL (TNFSF10) were found to be consistently upregulated in the areas with severe tissue damage. We used published spatial and single cell SARS-CoV-2 datasets to confirm our findings in the lung tissue from additional cohorts of COVID-19 patients. The resulting model of severe COVID-19 immune-mediated tissue pathology may inform future therapeutic strategies. One Sentence Summary Spatial analysis identifies IFNγ response signatures as focal to severe alveolar damage in COVID-19 pneumonitis.

The histological and immune cell landscape within COVID-19 lung tissue from three patients with fatal disease was investigated to establish the extent of intra-tissue variation in cellular pathology (mean specimen area 1.78cm 2 ; Supplementary Table 1 for patient information). Each sample featured a spectrum of DAD, from mild to severe, comprising pneumocyte hyperplasia, hyaline membrane formation, and interstitial expansion. There was a non-uniform distribution of immune infiltrates ( Figure 1A Immunofluorescent staining for CD3, CD68 and pan-cytokeratin (panCK) enabled the identification of lymphocytes, macrophages, epithelial cells, and general tissue architecture ( Figure   1B). Based on this staining we selected 46 areas of interest (AOIs) with a broad range of inflammatory features (Figure 1B and 1C). DAD histopathological severity was assessed by two independent pathologists and each AOI categorised as having either mild/moderate injury with some conservation of alveolar architecture or severe injury with a loss of alveolar structure and significant inflammation. In addition, for each AOI we quantified the number of CD3+ lymphocytes, CD68+ macrophages and cells (by nuclear stain). The SARS-CoV-2 nucleocapsid (N) protein positivity of adjacent tissue from sequential sections was also recorded. Areas of interest with severe DAD were associated with higher frequencies of T cells, although CD68+ macrophage numbers were not consistently increased ( Figure 1C Figure 1D). Genes upregulated in the severe areas 7 showed a significant over-representation (BH adjusted p<0.05) of gene ontology (GO) biological process terms related to T cell activation and differentiation, antigen presentation, cytokine production, cytotoxicity, and response to interferon gamma ( Figure 1E). By contrast, genes enriched in areas of mild damage (n=40, BH adjusted p<0.05, |fold change| >1.5) had a significant over-representation (BH adjusted p<0.05) of pathways associated with wound healing, epithelial migration and extracellular matrix organisation ( Figure 1D).

Subhead 2: Network analysis implicates CD8 + T cells, mononuclear phagocytes and active TLR, interferon and IL-1 signaling in COVID-19 lung inflammation
To perform an unbiased exploration of the cellular and phenotypic variation present in the set of the profiled AOIs, we employed weighted gene correlation network analysis (WGCNA) (13). This  Table 4). We began a systematic characterization of the identified modules by investigating their association with specific cell types.
To do so we correlated expression of the module eigengenes (the representative module expression patterns) with separate cell type abundance estimates that were determined for each AOI by automatic cell type deconvolution analyses (14) (Figure 2A  To confirm the involvement of specific cell types we examined the correspondence of the cell type abundance predictions with the expression of a curated set of known cell type marker genes. While we generally saw a good agreement between the marker genes and the automatic predictions of cell type abundance, correlated expression of MPO, ELANE and CTSG suggested that neutrophils may be present in some of the AOIs (Supplementary Figure 5A-B). To help resolve this discrepancy we applied mass-cytometry imaging to sequential sections and examined the areas aligned with selected severe AOIs. This analysis revealed the presence of a substantial number of CD15 + neutrophils in an area of severe damage from patient B, whose tissues also showed strong staining 8 for the SARS-CoV-2 N protein. Additionally, this analysis confirmed the expected presence of CD68 + macrophages and CD8 + T cells within severe areas, with only a small number of CD4 + cells being detected (Supplementary Figure 6). Next, we identified biological processes and pathways over-represented amongst the gene members of each WGCNA module (Figure 2B and   Supplementary Table 5) and explored the correlation of the module eigengenes with genes involved in immune signaling, the complement system, viral infection and fibrosis ( Figure 2C).
Based on these analyses we named each module according to its cell type associations and biological pathway enrichments.
As expected for lung tissue, we found a set of stromal modules representing (i) "Epithelial cells" (containing EPCAM), (ii) "Type 2 pneumocytes" (containing the surfactant encoding genes SFTPB, SFTPC and SFTPD), (iii) fibroblasts ("Fibroblast phenotype", containing COL1A1, COL3A2, COL5A1 and THY1) and (iv) "Vasculature" (containing CDH5, THBD, ENG), which all showed corresponding pathway and cell type associations (Figure 2A Figure 3C). These comprised (i) an "Alveolar macrophage" module that also displayed high expression correlation with the macrophage receptor MARCO and enrichment for the "phagocytosis, engulfment" GO biological process, (ii) a "Macrophage identity" module that also showed high expression correlation with the mannose receptor MRC1 (a marker of alternative "M2" macrophages) and SIRPA, and (iii) an "Antigen presentation" module that showed a high correlation with predicted inflammatory monocytederived macrophage (monoMac) cell abundance (Figure 2A, Supplementary Figure 4D) along with high expression correlation with SIGLEC1 and C1QA/B (Supplementary Figure 3B). We also noted a "IFITM2/HSP/ECM" module that contained MERTK and PDGFRA and a module with pathway enrichments for "Apelin/mTOR signalling".   Critical COVID-19 is associated with massive lung immune cell infiltration, and in keeping with this we discovered modules of genes with clear associations with lymphocytes and mononuclear phagocytes. The "Cytotoxicity and T cells" module was associated with the GO "Interferon-gamma production" biological process, as well as the KEGG "T cell receptor signalling" and "Natural Killer mediated cytotoxicity" pathways ( Figure 2B). The expression of this module was positively associated with predicted presence of CD8 + cytotoxic T cells, NK cells and activated dendritic cells (DCs) (Figure 2A) as well as the presence of CD3 + cells (Supplementary Figure 3C). This module contained known T cell markers including CD3D/E, CD2 and CD8A as well as genes associated with cytotoxicity such as PRF1 and GNLY (Supplementary Figure 3B). It also showed a high correlation with the expression of chemokines including CXCL9, CXCL10 and CXCL11, and had the highest correlation with the expression of IFNG ( Figure 2C) as well as the "Interferon responses" module (Supplementary Figure 3D). The "Antigen presentation" module was 11 associated with areas of high CD68 expression by immunofluorescence microscopy (Supplementary Figure 3C) and showed positive associations with the "Macrophage identity", "Cytotoxicity and T cells", "Interferon response" and TLR signalling related gene modules (Supplementary Figure 3D). The "TLR signaling and monocytes" module contained the gene encoding the classical monocyte marker CD14 along with CD163 (Supplementary Figure 3B). In addition, we discovered several modules that contained genes associated with innate inflammation including an "Interferon responses" module, a "TLR and IL-1 signaling" module and an "IL-1 response: IL-6/IL-8" module that were named according to their pathway enrichments ( Figure 2B) and correlations with chemokine/cytokine expression ( Figure 2C). The "Vasculature" and "IL-1 response: IL-6/IL-8" modules contained genes associated with mature neutrophils (ELANE, MME, Figure 3B). The "Vasculature" module was also associated with the GO.BP platelet degranulation pathway ( Figure 2B). Finally, we found three modules that were associated with more general cellular processes. The presence of these modules, which we termed "Cell cycling", "Chromatin remodelling" and "Hypoxic response" was indicative of immune cell proliferation and oxygen stress.  (15), as well as IL-32 which is known to stimulate TNFα and IL-6 secretion from macrophages (16), and CCL19 which acts via CCR7 to promote DC and central memory T cell migration (17). Next, we sought to assess the extent to which of the discovered features of severe DAD might represent common features of critical COVID-19. Within AOIs from the each of three patients, we found that severe DAD was consistently associated with the expression of the "Interferon responses", "Cytotoxicity and T cells", "Chromatin remodelling", "Antigen presentation", and TLR signaling ("TLR and IL-1 signalling" or "TLR signaling and monocytes") module eigengenes ( Figure 2E). The "IL-1 response: IL-6/IL-8" module did not show an obvious correlation with DAD severity, while higher expression of gene modules associated with "Epithelial cells" and "Vasculature" were associated with the lowest DAD scores ( Figure 2E). To investigate commonalities and differences of the tissue pathology of the three patients in more detail, we repeated the 'cellular phenotype network analysis' separately for the AOIs of each patient. This 13 analysis demonstrated that a core association between mono-Mac, CD8 + T cells, the expression of IFNγ target chemokines CXCL9, CXCL10 and CXCL11, the expression of IL32, the expression of the regulator of extrinsic apoptosis TRAIL (TNFSF10) and the expression of B cell activating factor, BAFF (TNFS13B) was present in the lung tissue from all three patients (Figure 3A-C).

MPO, CTSG) (Supplementary
To help prioritise immune signalling factors involved in severe DAD we ranked secreted immune The elevated expression of TRAIL in regions of severe DAD is consistent with a previous reports of apoptosis pathway upregulation in alveolar areas in late stage COVID-19 (6). The correlation of antigen presentation related genes with DAD severity (Figure 2E and Supplementary Figure 8) suggested a directed cytotoxic response, but we did not find an obvious association between the severity and the levels of SARS-CoV-2 N protein RNA (Supplementary Figure 1E). This suggested that CD8T/NK cell activation might be triggered by APC presentation or PRR (e.g. relative to healthy controls along with increased expression of CXCL9, CXCL10, CXCL11 and BAFF in myeloid cells and neutrophils ( Figure 4A). We also noted a broad increase in STAT1 and TRAIL expression along with induction of CCL3 and CCL4 expression in macrophages and neutrophils. The cytotoxic molecule encoding genes GNLY and PRF1 were upregulated in NK and T cells and neutrophils, which were not captured in the healthy controls, constituted an additional source of S100A8/9. The cellular sources and expression levels of these genes were consistent with those observed in other lung tissue single-nuclei and BALF single cell datasets (Supplementary Figure 11A, B). Inspection of a previously published spatial analysis of alveolar tissue (which did not discriminate between mild and severe DAD) confirmed up-regulation of STAT1, CXCL10, BAFF and TRAIL in COVID-19 (Supplementary Figure 11C). Increased expression of these genes in COVID-19 was also observed, albeit weakly, in a bulk analysis of macrophages purified from BALF (Supplementary Figure 11D). Together, these previously published observations and our targeted spatial analysis suggest the existence of a cellular circuit in which IFNγ production by T/NK cells drives CXCL9/10/11 and BAFF production by myeloid cells in COVID-19. If such a 17 circuit exists, we reasoned that the expression levels of these cytokines should co-vary in the relevant cell types in lung tissue across COVID-19 individuals. To investigate this possibility, we analysed single-nuclei data from the lung tissue of n=16 SARS-CoV-2 infected COVID-19 autopsy donors (7) (Supplementary Figure 12). In keeping with the existence of the proposed circuit, levels of IFNG in per-donor T/NK nuclei pseudobulks showed significant correlations (p <0.05) with the levels of TRAIL, BAFF and CXCL10 in the paired myeloid nuclei pseudobulks ( Figure   4B).  (19)(20)(21)(22)(23)(24). While much has been learnt from analysis of blood and BALF (10,18,(25)(26)(27)(28)(29), study of tissue is needed to understand the cellular causes of COVID-19 associated lung damage. Further insight has been gained from single-cell and single nuclei approaches (7, 30) but these approaches do not retain spatial information, which is vital for deciphering the interplay between different cell types. Initial applications of spatial proteomics and transcriptomics have started to reveal the spatial landscape of lung damage in COVID-19 (6, 7) but molecular details of the signalling circuits that perpetuate pathology remain to be fully elucidated. In this study, we sought to generate new insights by applying network-based analysis approaches to the analysis of rich spatial transcriptomics data generated with Nanostring GeoMx DSP platform. To do so, we used correlation networks to integrate WGCNA module eigengenes, cytokine gene expression levels and computationally predicted cell type abundances. Use of this flexible and extensible 'cellular phenotype network analysis' approach uncovered new links between the cell types, biological pathways and cytokines that are associated with lung tissue damage in severe COVID-19.
Our investigations found substantial differences in the cellular and molecular pathology of the lung tissue sampled from the three patients we studied. Tissue from patient A displayed a stronger 20 signature of type 2 pneumocyte and alveolar macrophages. This patient also had the highest expression of a gene module associated with hypoxic response, an observation which was not unexpected given the detection of low oxygen response and p53 stress pathways in the BALF of critical COVID-19 patients (31). In contrast with that from the other two cases, the tissue from patient B showed the strongest interferon response and TLR and IL-1 signalling signatures which corresponded with immunohistochemical evidence of high levels of viral infection and the presence of neutrophils. This is consistent with the idea that neutrophil extracellular trap formation (NETosis) may contribute to ongoing inflammation in some COVID-19 patients (32). Finally, the samples from patient C showed a lower expression of interferon and hypoxic response signatures and were distinguished by elevated expression of a vasculature associated gene module. Despite these broad differences, our initial within patient analysis revealed a shared association of severe DAD in COVID-19 with interferon signalling, cytotoxicity and T cells, cell proliferation and antigen presentation related genes.
To further elucidate common features of severe COVID-19 we performed within-patient 'cellular phenotype network analysis' and explored the reproducibility of our findings using data from published single cell, single nuclei and spatial transcriptomic studies of COVID-19 tissues (6,7,18,25). The results provide the basis of a model of severe lung tissue damage in COVID-19 in which IFNγ production by CD8+T and NK cells (i) activates macrophages and other innate immune cells and (ii) induces expression of CXCL9/10/11 (15,33) promoting further recruitment of CXCR3 + immune cells (including NK and cytotoxic T cells) in into lymphoid-rich areas ( Figure   4C). The presence of cytotoxic lymphocytes and elevated expression of TRAIL (TNFSF10) suggest that severe lung damage in COVID-19 may involve cytolysis and extrinsically regulated apoptosis.
In keeping with this model, myeloid cell dysregulation is a hallmark of severe or progressive COVID-19 infection (34)(35)(36)(37) and there is strong evidence for increased numbers of macrophages in 21 COVID-19 lung tissue (6,7,30). An increased ratio of CD14 + HLA-DR lo inflammatory monocytes to tissue-resident alveolar macrophages has also been noted (18) and macrophage hyperactivation by persistent IFNγ production has been previously suggested to be a possible mechanism in COVID-19 (38). Consistent with the proposed model, despite the well-characterised peripheral blood T and NK cell lymphopenia (26), numbers of CD8+ T cells in COVID-19 lung tissue are comparable to those found in healthy individuals, and higher than those found in pneumonia (6).

NK cells are less abundant than CD8+ T cells in the lung tissue of COVID-19 patients but appear
to show an increase in mild disease that is reduced to, or below, healthy levels in severe cases (6,7,18). Our data extends these observations by showing that, within the lung tissue, cytotoxic CD8+ T cells can localise to interstitial immune cell infiltrates with inflammatory, likely monocyte-derived macrophages and neutrophils within areas of severe COVID-19 associated damage. Based on our transcriptomic analysis we consider that such regions are also likely to contain cytotoxic NK cells but did not find transcriptional or immunohistochemical evidence for CD4 T cell involvement. Our model is likely to be incomplete because important limitations of our work include the number of patients studied, the targeted panel used for the transcriptomic analysis and the limited spatial resolution of the GeoMx DSP platform. Further large-scale studies of COVID-19 lung tissue using higher-resolution proteomic and whole-transcriptome spatial platforms to complement the throughput and targeted sampling that the GeoMx DSP platform affords will be essential for fully deciphering the fine cellular and molecular details of inflammation and severe tissue damage in this disease (39).
A key question that arises from the proposed model is the nature of the upstream mechanism(s) by which CD8+ T and NK cells are stimulated to release IFNγ in areas of severe damage. In a similar circuit proposed by Grant et al. based on the analysis of BALF, activation of SARS-CoV-2 reactive T cells in lung alveoli was proposed to be sustained by continued SARS-CoV-2 infection of 22 recruited monocyte derived macrophages (25). However, the fact that only one of the three patients studied here showed convincing evidence of viral infection in the tissue suggests that viral antigens may not be the only trigger for cytotoxic lymphocyte activation in severe cases. In support of this hypothesis, a similar observation was made in a study of lung tissue from fatal COVID-19 cases where "virus-independent immunopathology, rather than direct viral cytotoxicity" was proposed to be a primary pathogenic mechanism (8). Additionally, it has also been reported that the CD8+ T cell repertoire is more diverse in BALF from severe COVID-19 patients than it is in moderate cases (18). Together, these observations suggest that in severe COVID-19, activation of CD8+ T cells in lung tissue may be spatially uncoupled from, and inappropriately sustained after, clearance of viral infection. Candidate mechanisms for this process include exposure to endogenous DAMPs and/or pro-inflammatory cytokines. Possible endogenous DAMPs include the release of S100A8/9 from macrophages (and/or neutrophils) killed by cytotoxic lymphocytes or TRAIL mediated extrinsic apoptosis if not adequately cleared by phagocytosis (40). With respect to the possible involvement of inflammatory cytokines, we noted prominent expression of BAFF (TNFSF13B) in regions of severe DAD, where, based on inspection of single cell data (18), the most likely source was the myeloid cells. It is known that BAFF can promote CD4 + T cell IFNγ production and CD8 + T cell cytotoxicity in chronic obstructive pulmonary disease (COPD) (41), supporting the concept that it may be part of a positive feedback loop that sustains non-specific cytotoxic lymphocyte activation in COVID-19. It is also likely that elevated levels of IL-1, IL-6 and TNF may contribute to dysregulation of CD8 T cells, NK cells and macrophages as part of the so-called 'cytokine storm' (26,42).
Currently deployed therapeutics such as dexamethasone and the anti-IL-6 therapy tocilizumab are effective but not curative in all patients (43). While the cases that we examined provided some evidence of a role for IL-6 in areas of milder damage, we did not find an obvious link between expression of this cytokine and severe tissue pathology. Our data, suggests that therapeutic 23 targeting of immune cell recruitment via the CXCL9, CXCL10, CXCL11/CXCR3 axis may be a valuable therapeutic strategy for resolution of inflammation in severe COVID-19 where it has become uncoupled from viral clearance. Systemic and lung tissue upregulation of these IFNγ induced cytokines in COVID-19 has been noted in many studies, and this axis has already been proposed a therapeutic target by others for COVID-19 and other serious diseases including cancer (6,10,15,18,26,38,44). Less attention has been paid to BAFF (TNFSF13B) which has been shown to be significantly up-regulated in COVID-19 patient plasma, and for which monoclonal blocking antibodies have been developed (45,46). While further study of the role of this interesting cytokine in COVID-19 is required, our data suggest that it may be a valid therapeutic target in severe cases. In two of the three patients studied we also noted a robust up-regulation of IL-32 in areas of severe damage, and the role of this cytokine, which has important functions in anti-viral responses and can induce expression of pro-inflammatory cytokines may warrant further study in . Overall, the heterogeneity of the tissue pathology that we observed both within and between the patients studied underscores a need for personalised approaches where the choice of therapy (or therapies) is guided by careful assessment of the stage of disease and viral clearance.
For example, while innate immune activation via treatments such as intranasal IFNB1α/β (48) are likely to be important for patients unable to clear the virus, our model predicts that they may exacerbate non-specific lymphocyte and myeloid cell activation and tissue damage in virus-free lung tissue and in patients who remain critical after successful viral clearance. As a more general strategy, our findings strongly support the use of combinatorial regimes that pair therapeutics directly targeting the virus such as molnupiravir, remdesivir, ritonavir and recombinant soluble ACE2 (49,50) with those that can temper the dysregulated host immune responses that drive tissue damage.

Design
To delineate the tissue-specific immune pathology in severe COVID-19, we assessed the transcriptomic profile across a spectrum of DAD in well-preserved tissue samples obtained at the point of death from three COVID-19 patients. Medical records of patients were retrospectively reviewed (51) and three COVID-19 patients selected for in-depth analysis based on their clinical manifestation of ARDS (Supplementary Table 1), typical COVID-19 histology (Supplementary Table 2) with a 4-5 score on the Brescia-COVID Respiratory Severity Scale, and a lung-restricted presence of SARS-COV-2 (absence in heart, liver and kidney biopsies). None of the patients were vaccinated against SARS-CoV-2. At least 14 AOIs from each COVID-19 tissue were selected for analysis, spanning on average 0.2mm 2 (range: 0.05 to 0.33mm 2 ) with exclusion of empty space.
Areas were selected to represent the spectrum of alveolar injury within each tissue covering regions of a) mild/moderate injury with some conservation of alveolar architecture and b) severe injury with a loss of alveolar structure and significant inflammation. Selected AOIs occupy alveolar and interstitial spaces, except from A_16 and A_17 containing bronchiolar epithelium. The severity grade of each AOI was confirmed post-hoc by 2 pathologists.

Ethics statement
This study was approved by the ethics committee of the University of Navarra, Spain (Approval Patients and tissue processing 25 Post-mortem lung tissues were obtained through open biopsy at the point of death and processed as described (51). In brief, tissues were immediately fixed in neutral buffered formalin for <24 hours and then paraffin embedded. Sections (5μm each) were cut for H&E staining, ISH and DSP analysis. Six pathologists reviewed the histology in and agreed on the gross histological characteristics (Supplementary Table 2). RNA was extracted from 4-8 x 5μm sections for quantification of SARS-CoV-2 nucleocapsid (N) and envelope protein (E) transcripts (Supplementary Table 2).

In situ hybridization
In-situ hybridization was conducted using the RNAscope2.5 LS Reagent Red Kit, according to manufacturer's instructions and using the Leica BOND-RXm system. Deparaffinisation and heatinduced epitope retrieval were performed with BOND Epitope Retrieval Solution 2 (ER2, pH 9.0) for 25 minutes at 95°C. Hybridization for SARS-CoV-2 RNA was carried out using the

Analysis of immunofluorescent images for cell counts
The number of nuclei, CD3 + and CD68 + cell counts were determined using CellProfiler software.
A pipeline was designed to quantify circular objects within RBG files for each AOI (53). Global manual intensity thresholds were set for object identification of nuclei and CD3 + cells, whilst 27 CD68 + cells were identified by adapted intensity thresholds. The efficacy of object identification for each AOI was visually confirmed.

Quality control and pre-processing of GeoMx transcript expression data
Quality control and initial data exploration was conducted using the GeoMx DSP Analysis Suite.
Sequencing quality per AOI was examined and an under-sequenced area (B_04) with zero deduplicated reads was excluded. Expression of each transcript was measured by 5 or more  WGCNA modules were further characterized by assessing the correlation of the module eigengenes with a) histological severity (Pearson correlation) b) predicted cell abundances (as estimated using Spatial Decon; Spearman correlation) and c) that of selected genes (Spearman correlation, R library Hmisc).

Cell deconvolution of GeoMx transcript expression data
Deconvolution of cell types from the gene expression data was performed using the SpatialDecon R library (14). For this analysis the full matrix of quantile normalized gene counts (without correction), mean negative probe counts and the cell profile matrix 'Lung_plus_neut' were applied as input. The cell profile matrix retrieved from the package was generated from the Human Cell Atlas lung scRNAseq dataset and appended with neutrophil profiles as described in Desai et al.
(2020) (5). For correlation analyses, the abundance of cell types is reported for cell types present with > 2 relative abundance in > 2 AOIs. Complete cell type output is shown in Supplementary   Figures 5, 7, 9, 10.

Construction of AOI group correlation networks
We constructed correlation networks to investigate the relationship between the WGCNA modules, the estimated cell abundances (from SpatialDecon) and the expression levels of immune signalling genes for each of the n=5 AOI spatial groups (R igraph library; layout = layout_with_dh). Nodes were scaled in size and/or colour according to cell abundance, normalized gene expression or WGCNA module eigengene expression. Drawn edges represent significant correlations (p < 0.05) and were weighted according to the value of the correlation coefficient (Spearman's Rho).
WGCNA modules were included in the network if they had a median module eigengene expression >0 within the relevant group of AOIs. Immune signalling genes from the KEGG 'cytokinecytokine receptor interaction' pathway (hsa04060) (human) were included if showed a median expression above the expression detection threshold (as defined in the pre-processing section) in a given AOI group. Cell type abundance estimates were included if they had a median abundance of >2 the given AOI group. For inclusion in the network, nodes (modules, cell types and genes) had higher median values than the thresholds stated above and at least 1 positive correlation to another node type.