Multiplexed digital spatial protein pro�ling reveals a unique protein signature for advanced liver �brosis

Background and Aims: Intrahepatic mononuclear phagocytes (MPs) are critical for the initiation and progression of liver �brosis. In this study, using multiplexed digital spatial protein pro�ling, we aimed to derive a unique protein signature predicting advanced liver �brosis. Snap-frozen liver tissues from various chronic liver diseases were subjected to spatially de�ned protein-based multiplexed pro�ling (Nanostring GeoMX™). Single-cell RNA sequencing analysis was performed using Gene Expression Omnibus (GEO) datasets from normal and cirrhotic livers. Results Sixty-four portal regions of interest (ROIs) were selected for the spatial pro�ling. Combined analysis of single-cell RNA sequencing data from GEO datasets (GSE136103) and spatially-de�ned, protein-based multiplexed pro�ling revealed that most proteins upregulated in F0–F2 livers in portal CD68 + cells were speci�cally marked in tissue monocytes whereas proteins upregulated in F3 and F4 livers were marked in SAMacs and tissue monocytes. Internal validation using mRNA expression data with the same cohort tissues demonstrated that mRNA levels TREM2, CD9, and CD68 are signi�cantly higher in livers with advanced �brosis. Using the results from the CD68 + area, a highly sensitive and speci�c immune-related protein signature (CD68, HLA-DR, OX40L, phospho-c-RAF, STING, and TIM3) was developed to predict advanced (F3 and F4) �brosis.


Introduction
Liver brosis can eventually progresses to cirrhosis, which is a major concern to cause approximately 2% of deaths globally. [1] Non-resolving liver injury or liver brosis can cause cirrhosis, a condition that can progress to hepatocellular carcinoma. Furthermore, decompensated liver cirrhosis is accompanied by serious complications resulting from portal hypertension and/or opportunistic infections and leads to shorter life expectancy and deterioration in quality of life. [2] In an effort to alleviate the burden of chronic liver disease, many studies have sought to nd novel biomarkers of advanced liver brosis. [3] A very recent study from our group also recently presented a promising immune-related gene signature for predicting advanced brosis. [4] Although several approaches to reduce brosis are still under investigation, no drug has yet been demonstrated to be effective in reversing brosis.
Liver brosis involves many non-parenchymal cells (NPCs) including immune cells, endothelial cells, and mesenchymal cells. [5] These cells interplay with each other in scarring tissue, also known as the brotic niche. Among immune cells, intrahepatic mononuclear phagocytes (MPs) are thought to play a crucial role in exacerbating hepatic in ammation and brosis. The mononuclear phagocytic system (MPS) is comprised of circulating monocytes, dendritic cells, and tissue-resident macrophages, also known as Kupffer cells (KC) in the liver.
[6] KCs, which reside in hepatic sinusoids, are ontogenically different from other circulating monocytes. Originating from the embryonic yolk sac, they are capable of self-renewing independently of circulating monocytes. [7] KCs play an important role in innate immunity as gatekeepers and drive the in ammatory response to liver injury. By contrast, circulating monocytes are likely dispensable for replenishing intrahepatic macrophages in homeostasis. However, in the setting of hepatic in ammation, massive in ltration of monocyte-derived macrophages (MoMFs) to the liver is triggered.
[8] MoMFs, which originate from bone marrow, are recruited from the portal triad to the injured liver by cytokines or chemokines secreted by activated KCs including IL-1β, tumor necrosis factor α, CCL2, and CCL5. [9] As in ammation prolongs, MoMFs then may differentiate into pro-brogenic TREM + CD9 + scar-associated macrophages (SAMacs) and participate in the process of scarring, which leads to liver brosis.
[10] SAMacs is known to play a critical role in the initiation and progression of liver brosis and continue to expand in the brotic liver while brosis proceeds.
Overall, these ndings suggest that MPs play a critical role in the process of liver brosis, thus it is feasible to nd out biomarkers based on the MPsassociated transcriptomes. accuracy to predict liver brosis are still lacking. Here, using multiplexed digital spatial pro ling, we aimed to nd out the novel monocyte/macrophage-based protein biomarkers for advanced brosis.

Patients
This study included benign liver tissues from 83 patients who underwent surgical resection for hepatocellular carcinoma. All tissues were obtained between 1996 and 2015 at the Ajou Medical Center (AjouMC). The inclusion criteria were as follows: (1) benign liver tissues with chronic liver diseases (2) tissues collected during the surgical procedures such as liver resection. The retrospective study protocol was approved by the Institutional Review Boards of Ajou Medical Center (AJIRB-BMR-KSP-18-444) and The Catholic University of Korea (XC20EEDI0034). The brosis stages of every tissue sample enrolled in this study were determined by one pathologist (E.S.J), using the METAVIR scoring system.

RNA extraction and gene expression assay
Total RNA was extracted from the liver tissues using a RNeasy Mini Kit (QIAGEN, Hilden, Germany) with DNase I treatment (QIAGEN). Gene expression pro les were analyzed using nCounter MAX (NanoString Technologies, Seattle, WA, USA). The nCounter PanCancer Immune Pro ling Panel (NanoString Technologies) was used for gene set pro ling as previously described. [4] Tissue microarray (TMA) construction TMAs were constructed using 2-mm diameter cores punched from formalin-xed, para n-embedded (FFPE) blocks. The TMA blocks were sectioned 5-mmthick. Six TMA slides were constructed with 15 cores placed in a 5×3 arrangement on each slide, and its representative ROIs are presented in gure 2a.
Digital spatial pro ling (DSP) assay Protein expression pro les were analyzed using GeoMx DSP (NanoString Technologies, Seattle, WA, USA). The TMA slides were stained with a mixture of detection and morphological markers. Morphological markers included Syto13 for nuclei, CD68 for macrophages, CD3 for T cells, and alpha-SMA for smooth muscle. The detection antibodies comprised one core panel and six modules of the GeoMx assay (GeoMx immune cell pro ling panel, GeoMx io drug target module, GeoMx immune activation status module, GeoMx immune cell typing module, GeoMx pan-tumor module, GeoMx cell death module, and GeoMx MAPK signaling module). In total, 88 ROIs were selected around the portal tract. Each ROI was divided into CD68+, CD3+, and SMA+ areas. Probes attached to the detection antibodies were collected sequentially from the CD68+, CD3+, and SMA+ areas.

Analysis of DSP data and validation of the protein signatures
The left part of the gure 1 shows a owchart of protein signature development. Differentially expressed proteins were analyzed by comparing brosis stages 0, 1, 2, 3, and 4 in CD68+ areas of samples. Of the differentially expressed proteins, those that were signi cantly different in the logistic regression analysis were included in the protein combination. For the derivation of protein signatures, logistic regression coe cients for each protein were identi ed and weighted according to protein expression values. Candidate protein signatures (AUC > 0.85, accuracy > 90%, P < 0.05) were validated by k-fold cross-validation to identify the optimal protein combination. The patients were randomly separated by 2-folds (training and test sets) 300 times.

Single-cell RNA sequencing analysis
The right part of gure 1 depicts a ow diagram of the scRNA-seq analysis. The analysis was performed using the GSE136103 dataset and the Seurat package version 4.0.5 (https://satijalab.org/seurat/index.html). Pre-processing followed the GSE136103 method. For shared nearest neighbor clustering, variable features were determined using the variance stabilizing transformation (VST) of the selection method and 2,000 variable counts. After scaling the data, principal component analysis (PCA) was performed using the identi ed variable features. With the analyzed principal components, the optimal dimensions were analyzed using an elbow plot. According to the determined dimensions and principal components, single cells were clustered and visualized as a t-distributed stochastic neighbor embedding (t-SNE) graph. Clusters were identi ed using markers used to divide the clusters. Primary clustering was performed using data from ve healthy and ve cirrhotic patients. Secondary clustering was performed using cells presumed to be MPs. To identify the distribution in which cells express genes that code for differentially expressed proteins detected by DSP, the corresponding genes were marked in the secondary cluster using SCINA (Semi-supervised Category Identi cation and Assignment).

Statistical analysis
Continuous data are presented in means with standardized deviation and categorial variables are expressed as number and percentage. The categorial variables between brosis stage and clinicopathological variables were assessed using the chi-squared test or Fisher's exact test and continuous variables between groups were evaluated using Wilcoxon rank-sum test. The predictive accuracy of the threshold values for classifying brosis stages 0-2 and stages Of the 88 ROIs, 64 with low CD3 expression were presumed to be "non-in ammatory" regions whereas 24 ROIs with high CD3 were classi ed as "in ammatory" regions ( gure 1). We excluded in ammatory ROIs because the phenotypes of the immune cells may not actually re ect the brogenesis process but the liver injury process. The baseline characteristics of the enrolled patients are shown in Table 1. A total of 31 ROIs from the enrolled patients were classi ed as having early stage brosis (seven as F0, 17 as F1, and 7 as F2), and 33 ROIs were classi ed as having advanced brosis (15 as F3 and 18 as F4). The enrolled patients were predominantly male, and the mean ages were 54.23 and 50.91 years for patients with brosis stages 0-2 or stages 3-4, respectively. The most common etiology of chronic liver disease was hepatitis B infection, which accounted for 67.74% and 75.76% of patients with brosis stages 0-2 or stage 3-4, respectively. Because nucleotide/side analogs for HBV therapy are reimbursed for patients with HBV DNA > 2000 IU/mL during cirrhosis in Korea, a larger proportion of patients with advanced brosis received antiviral therapy compared to those with early brosis (9.68% vs. 24.24%, p=0.2255). There were no signi cant differences between the early and advanced brosis groups regarding albumin, AST, ALT, and platelet levels.

Multiplexed DSP of protein expression level according to the brosis stage
Representative ROIs stained with CD3, CD68, and SMA according to the brosis stage are depicted in Figure 2A. There were liver tissues with various brosis grades, namely F0 (n = 7), F1 (n = 17), F2 (n = 7), F3 (n = 15), and F4 (n = 18). A volcano plot for protein marker expression is depicted in Figure 2B. The gure compares protein expression levels in early (F0, F1, and F2) and advanced (F3 and F4) brosis. The results are summarized in Table 2. Compared to early stage brosis, advanced brosis showed upregulation of CD68 and HLA-DR, which were highly expressed in the livers with advanced brosis with 1.50-and 1.32-fold changes, respectively. By contrast, protein markers other than CD68 and HLA-DR, such as phospho-c-RAF, stimulator of interferon genes (STING), OX40 ligand (OX40L), V-domain IgG suppressor of T cell activation (VISTA), pan RAS, and T cell immunoglobulin and mucin domain-containing protein 3 (TIM3) were downregulated in brosis stages 3 and 4, with fold changes up to 2.65 ( Table 2).
Validation of DSP protein analysis using mRNA expression data from snap-frozen livers using NanoString nCounter MAX system We also applied NanoString nCounter MAX mRNA expression analysis using the same liver tissues (snap-frozen) used in DSP analysis to validate our protein data. Figure 3D delineates the differentially expressed genes between brosis stage 0-2 and 3-4. A Previous report demonstrated that TREM2 and CD9 are selectively upregulated in SAMacs and can be used as protein markers for this MP subset. [15] In our NanoString nCounter MAX mRNA expression analysis, TREM2, and CD9 were also shown to be upregulated in F3-4 than F0-2 (TREM2, fold change=1.70, p=0.16; CD9, fold change=1.36, p=5.54×10 -4 ). Moreover, CD68, a protein that was newly identi ed to be upregulated in MPs in advanced brosis by DSP analysis, was also signi cantly upregulated by mRNA expression analysis (fold change=1.29, p=1.84×10 -4 ). CD74, also known as HLA-DR antigens-associated invariant chain, showed a tendency of higher expression in advanced brosis compared to early brosis although statistical signi cance was not met.
Predicting related immune cells in different brosis stage using single cell RNA sequencing database Next, we tried to visualize our marker proteins from the DSP results in a t-SNE map derived from a publicly available scRNA-seq dataset (GSE136103). By utilizing the non-parenchymal liver cell atlas established by Ramachandran et al. [15], we have marked our DSP driven proteins into the atlas to specify the upregulated subpopulation of MPs. We re-analyzed the dataset and classi ed MPs into ten clusters, as depicted in Figure 3A. Each cluster showed high similarities to the t-SNE map previously presented by Ramachandran et al. [15] Most protein markers upregulated in MPs of F0-2 by DSP were marked in clusters of tissue monocytes, as shown in left panel of gure 3B. Most protein markers upregulated in MPs of F3-4 by DSP were marked in clusters of SAMacs, KCs, and tissue monocytes ( gure 3B, right panel). The scaled gene expression in each cluster of MPs is shown in Figure 3C. CD68 and CD74, which were two genes that showed higher protein expression levels in advanced brosis in our DSP analysis, were also highly expressed in SAMacs1 and SAMacs2 rather than in tissue monocytes.
Protein signatures for the advanced brosis derived from the DSP analysis Using the DSP results from the CD68 + area, we identi ed unique immune-related protein signatures that re ect the advanced brosis. Table 3 shows the candidate protein signatures derived from the DSP. Of the ve candidate protein signatures, one that was composed of the genes CD68, HLA-DR, OX40L, phospho-c-RAF, STING, and TIM3 showed the highest positive predictive value for the advanced brosis, with an AUC in the ROC curve of 0.873 (0.791-0.955).
The p value of the protein signature by logistic regression analysis was 2.53×10 −4 . The positive predictive value of the developed protein signatures was 93.10, and the negative predictive value was 82.86. The sensitivity and speci city of this protein signature were 81.82% and 93.55%, respectively (Table 3).
We also evaluated the factors associated with high-grade brosis using logistic regression analysis (Table 4). Variables such as clinical parameters, pathological states, and protein signatures were included in the analysis. In the univariate analysis, the protein signature was the only factor associated with predicting advanced brosis (p=1.16×10 -6 ). Two variables with acceptable p-values, the protein signature and BMI, were included in the multivariate analysis.

Discussion
In this study, using multiplexed DSP protein pro ling, we have identi ed a novel protein signature that predicted advanced brosis with a high reliability.
Moreover, using scRNA-seq database, we demonstrated the phenotypical heterogeneity of portal MPs according to the brosis stages.
Liver brosis is a common pathological consequence of most chronic liver diseases. Fibrosis is associated with many NPCs, including in ammatory, endothelial, and mesenchymal cells. Numerous reports have elucidated the different phenotypes of NPCs depending on the presence or absence of liver brosis. [19,20] In addition, different roles of the MPs according to brosis stage have been suggested in several previous studies. [15] KCs, which dominate the hepatic macrophage pool, are central to intrahepatic immunological tolerance and provide an anti-in ammatory micromilieu to the liver during homeostasis. [9,21] However, in acute or chronic liver injury, KCs secrete CCL2 and thereby recruit circulating monocytes to the liver, which then differentiate into MoMFs. These MoMFs prevail during liver injury and stimulate stellate cells. [22] This results in excessive deposition of extracellular matrix, which leads to hepatic brosis. MacParland et al. mapped the cellular landscape of the human liver via scRNA-seq and reported that CD68 + macrophages have two distinct phenotypes that are classi ed as having pro-in ammatory or immune-regulatory roles. [12] Moreover, recent studies using scRNA-seq demonstrated that TREM + CD9 + SAMacs were derived from circulating monocytes and demonstrated a pro-brogenic phenotype. [15,23] Collectively, these studies proposed distinct phenotypes of intrahepatic cell populations through scRNA-seq and suggested that the pro-in ammatory phenotype of intrahepatic macrophage switch to anti-in ammatory or pro-brogenic phenotype during the process of liver brosis.
To the best of our knowledge, this is the rst study to show the different phenotypes of MPs between the early and late stages of liver brosis using spatially de ned protein-based multiplexed pro ling. The DSP transcriptome in our study was matched and analyzed using the publicly available RNA-seq dataset (GSE136103) from Ramachandran et al. [15] Tissue monocytes appeared to be highly abundant in the portal area of livers with brosis stage 0-2 whereas KCs, SAMacs, and tissue monocytes were highly abundant in the portal area of livers with brosis stage 3 and 4. In addition, mRNA expression analysis using nCounter gene expression assay showed higher expression level of representative markers of SAMacs, namely TREM2, CD68, and CD9 in advanced brosis, supporting the DSP results. These results are consistent with previous reports demonstrating the pro-brogenic phenotype of SAMacs, which are thought to be derived from circulating monocytes.
In this study, we also proposed a novel protein signature derived from DSP data for advanced liver brosis. The proposed protein signature was composed of six different proteins: CD68, HLA-DR, OX40L, phospho-cRAF, STING, and TIM3. CD68 and HLA-DR were upregulated in brosis stages 3 and 4 whereas the other four proteins were downregulated in the advanced stages compared to stages 1 and 2. CD68, a type 1 transmembrane glycoprotein of 110 kDa, is a representative macrophage marker. A recent study using mice livers demonstrated increase in CD68 + macrophages in advanced liver brosis or cirrhosis compared to normal liver. It also found that CD68 + macrophages were predominantly concentrated in scars during advanced brosis, suggesting its probrogenic role in the process of liver brosis. [4,24] HLA-DR is routinely used to identify macrophage lineages such as Kupffer cells and circulating macrophages. [25] HLA-DR is also a typical marker for macrophages and is widely present in antigen-presenting cells. It is known to be upregulated upon immune stimulation and is proposed to be a monocyte activation marker.
[26] Our very recent study demonstrated higher level of HLA-DR in intrahepatic monocytes in human NASH livers with advanced brosis than that with early brosis. [22] OX40L, phospho-cRAF, STING and Tim3, which are shown to be downregulated in the advanced stage of brosis, are known to be related to the in ammatory change of intrahepatic monocytes. OX40L, a member of the tumor necrosis factor superfamily, interacts with OX40 and is associated with the secretion of pro-in ammatory cytokines in the setting of non-alcoholic steatohepatitis in mice. [27] Raf kinase is thought to promote cell growth through the direct phosphorylation of mitogen-activated protein kinase (MAPK) and activation of its downstream signaling. [28,29] In terms of STING, it is an important innate immune protein that coordinates with multiple immune responses, including the induction of interferons. [30,31] In addition, many studies on the role of STING in the activation of Kupffer cell and non-alcoholic fatty liver disease (NAFLD) have been proposed. [32][33][34] Luo et al. [32] and Yu et al. [33] have demonstrated STING upregulation in mice NASH livers and have suggested its possible role as a therapeutic target on the liver brosis. Lastly, TIM3 is a surface marker for terminally differentiated T-cells, is also expressed in monocytes, and is thought to have a regulatory role in liver brosis. [35] Several studies have shown that expression level of TIM3 in monocytes are decreased in cirrhosis and also suggested that high level of TIM3 blocks monocyte activation. [36,37] The protein signatures are tissue driven biomarkers and for that reason, are di cult to apply in the everyday clinical practice. However, they can be considered as bridges in establishing liquid biopsy biomarkers to predict advanced brosis. Blood-based biomarkers may be identi ed by analyzing the miRNAs or gene signatures associated with our tissue-driven protein signatures.
This study developed the novel protein signature predicting advanced brosis with high reliability. Using DSP, we have speci ed the region of protein analysis to the periportal area where MPs are abundant. Furthermore, this is the rst study that used spatially de ned protein-based multiplexed pro ling to show the critical difference in the phenotypes of portal MPs between livers with early-or late-stage brosis, and the results were validated using internal gene expression data and publicly available scRNA-seq database. Further studies are essential to validate the proposed protein signature and to recreate the differences in the phenotype of MPs between early-and late-stage brosis.

Con icts of interest
The authors declare that they have no con icts of interest. The funders had no role in the study design; collection, analyses, or interpretation of data; writing of the manuscript; or decision to publish the results.
Author Contributions: S.H.B., P.S.S. and E.S.J were responsible for conceptuallization, methodology and supervision of the overall study. J.L. and C.M.K were responsible for visualization and writing the original manuscript. H.J.W. was responsible for resources. J.H.C., J.Y.P., and Y.S.Y. were responsisble for formal analysis, investigation and data curation. All authors have read and agreed to the nal version of the manuscript.