Characterisation of HBV and co-infection with HDV and HIV through spatial transcriptomics

Patient samples and ethics

Archived formalin-fixed paraffin-embedded (FFPE) liver biopsies from three patients with chronic HBV (PT1), one with HDV/HBV co-infection (PT2) and one patient with HIV/HBV (PT3) co-infection were analysed in this study. Participants were recruited at The Royal London Hospital, Barts Health NHS Trust and samples were obtained during routine diagnostic procedures following written informed consent. Table 1 includes the demographic and clinical parameters of each participant. Each biopsy was evaluated by a histopathologist and scored for the degree of fibrosis (Ishak) and Hepatitis Activity Index.

NanoString GeoMx digital spatial profiling

Slides were submitted for DSP and whole transcriptome sequencing, according to the manufacturer’s recommendations for GeoMx-NGS RNA BOND RX slide preparation (manual no. MAN-10 131-02). The panel of morphology markers was designed to include HBsAg ((Polyclonal) 488 1:200 from Bioss #bs-1557G-A488), CD45 (D9M8I) 594 1:100 from Cell Signaling Technology #13 917BF, CD3 (clone UMAB54; Origene #UM500048CF) and a DNA dye (Syto83dye; Invitrogen). Deparaffinisation, rehydration, heat-induced epitope retrieval (for 20 min at 100°C) and enzymatic digestion (1 μg/mL proteinase K for 15 min at 37°C) were carried on the Leica BOND-RX. Tissues were incubated with 10% neutral buffered formalin for 5 min and for 5 min with NBF stop buffer. The tissue sections were hybridised with the oligonucleotide probe mix (Whole Transcriptome Atlas) overnight, then blocked and incubated with the four fluorescently labelled markers for 1 hour. Tissue sections were loaded into the GeoMx platform and scanned for an immunofluorescent signal. A geometric strategy was employed for selection of ROIs, which were segmented on the basis of HBsAg staining and CD45 staining to capture the presence of immune cell infiltrate (online supplemental figure 1A–C). ROI metadata are included in online supplemental table 1. The DSP barcodes were UV-cleaved and collected for each selected ROI and counted. During the library preparation, the DSP barcodes were tagged with their specific ROI location and RNA target identification sequence. Sequenced oligonucleotides were then processed and imported back into the GeoMx DSP platform for integration with the slide images and ROI selections for spatially resolved RNA expression. FASTQ sequencing files were processed into digital count conversion digital files using Nanostring GeoMx NGS Pipeline software.

Quality control and normalisation

Quality control checks and data analysis were performed in the GeoMx DSP Data Analysis suite. All ROIs had acceptable sequencing quality metrics after examination of the raw, trimmed, stitched, aligned and deduplicated reads, and the sequencing saturation (>60%; online supplemental figure 2B, C). An expression value per gene target per ROI was generated in GeoMx DSP Analysis Suite through the compilation of multiple probes, within the NanoString GeoMx Human Whole Transcriptome Atlas probe mix, for each transcript.

The manufacturer’s recommended threshold for the detection of genes is referred to as the limit of quantification (LOQ); the LOQ is 2 geometric SD above the geometric mean. The ROIs from PT2, PT1 and PT3 tissues were filtered to include only the genes detected over the LOQ in 3 or more ROIs: 9672 genes. The data were then quantile normalised and the performance of the normalisation was assessed by the standardisation of housekeeping genes and distribution of total gene expression between ROI (online supplemental figure 2D, E).

Principal component analysis

The variance in the dataset was assessed by principal component analysis (PCA) of log₂(n+1) transformed of the quantile-normalised and gene-filtered expression values from the selected tissue samples and variance was associated with biological factors, such as tissue identity, estimated immune infiltrate and estimated HBsAg expression. The contributing genes to the PCA analysis are included in online supplemental table 2.

Over-representation analyses

Differential gene expression was calculated using linear mixed models with fixed effects of immune presence or HBsAg presence in R package Limma V.3.56.2. P values were unadjusted in the figures and both unadjusted and adjusted using the Benjamini-Hochberg method in online supplemental table 3. Over-representation analysis for Gene Ontology Biological Processes (GO.BP) terms was performed on genes showing >1.5-fold change and p value <0.05 using R package clusterProfiler (V.4.8.2). Analysis parameters included a p value cut-off of 0.05 and a q value cut-off of 0.1. Redundant and duplicated pathways were excluded from figures. For the complete list of pathways with both unadjusted and adjusted values, refer to online supplemental table 4.

Cell deconvolution

Cell deconvolution analyses were generated in the GeoMx DSP control centre, using the spatialdecon geoscript (V.1.1, updated April 2021) available at Nanostring’s Geoscript Hub (https://nanostring.com/products/geomx-digital-spatial-profiler/geoscript-hub/). The analyses were run using the adult liver landscape 10x matrix (https://github.com/Nanostring-Biostats/CellProfileLibrary). Cell abundances were estimated using the SpatialDecon R library, which performs mixture deconvolution using constrained log-normal regression. The 0.75 quantile-scaled data were used as input.23

Weighted gene correlation network analysis

The weighted gene correlation network analysis (WGCNA) was applied to the quantile normalised and filtered expression values, which were then transformed by log₂(n+1). WCGNA was applied to 37 ROIs. The parameter values were set as follows: minimum fraction of non-missing samples for a gene to be considered good=0.5 (default), minimum number of non-missing samples for a gene to be considered good=4 (default), minimum number of good genes=4 (default), cut height for removing outlying samples=80 (no samples removed), minimum number of objects on branch to be considered a cluster=2 (cutreeStatic), network type=signed-hybrid, soft power=5, adjacency correlation function=bicor, maxPOutliers=0.05 adjacency distance function=dist (network adjacency default), TOM type=signed, minimum module size (number of genes)=30, dissimilarity threshold used for merging modules: 0.25. The analysis identified 8 distinct modules labelled with colours, plus a module of 2780 unassigned genes (grey module). The WGCNA analysis was performed using R packages: WGCNA V.1.72-1, tidyverse V.2.0.0 and openxlsx V.4.2.5.2.

This analysis identified 8 modules of co-expressed genes ranging from 157 to 2280 genes per module. These modules were characterised by enrichment in GO.BP terms using R package clusterProfiler V.4.8.2, by correlation with cell type abundance determined through cell deconvolution against phenotypes identified by human liver single-cell RNA sequencing analysis and by correlation with estimates of CD45 and HBsAg immunofluorescent staining during GeoMx protocols (online supplemental table 5).

Study participants have discrete spatial signatures associated with viral and immune features

We studied FFPE liver tissue from three CHB treatment-naïve patients: HBV mono-infection (PT1); HDV/HBV co-infection (PT2) and HIV/HBV co-infection (PT3) (table 1). PT1 in the HBeAg− chronic hepatitis phase had fluctuating ALT levels (20–53 IU/mL) and HBV DNA (6850–190 000 IU/mL) over 2 years. Liver biopsy showed normal architecture with minimal inflammation. Mild steatosis was also present. PT2 with untreated HBV/HDV co-infection for 18 years, and HBeAg+ chronic hepatitis, had moderate levels of HBV DNA and HDV RNA and evidence of advanced fibrosis suggestive of cirrhotic transformation and moderate inflammation. PT3 with HIV/HBV co-infection diagnosed 1 year prior, and HBeAg+ chronic hepatitis had fluctuating ALT levels (151–348 IU/mL) and evidence of marked inflammation with retained liver architecture. Immunofluorescent staining for HBsAg, PTRPC (CD45) and DNA (nuclear stain) guided selection of ROIs for transcriptional profiling, revealing distinct/heterogeneous HBsAg patterns (figure 1A–C and online supplemental figure 1A–C) as previously described.22 24 PT1 displayed a ground glass appearance of staining, with areas of marginal/submembranous HBsAg staining (figure 1A). PT2 exhibited a clustered distribution of cytoplasmic HBsAg expression, resembling clonally expanded populations25 (figure 1B). PT3 had a diffuse membranous pattern of HBsAg expression involving the majority of hepatocytes indicative of high replicative HBV status (HBV DNA 235×10⁶ IU/mL) (figure 1C).

Request permission

We assessed immune infiltration and HBsAg expression in selected ROIs in each patient biopsy (figure 1A–C, online supplemental table 1). ROIs averaging 96 871 μm² (±44 284 μm² SD) with an estimated mean of 412 nuclei (±238 nuclei SD; online supplemental figure 2A) demonstrated excellent sequencing quality (online supplemental figure 2B, C) and minimal variation across the whole genome atlas/housekeeping transcripts (online supplemental figure 2D, E). 37 ROIs were categorised according to immune cell (CD45+) and HBsAg+ cell composition and gene expression enriched in each ROI was quantified by the GeoMx Human Whole Transcriptome Atlas assay. PCA indicated variability attributed to donor-specific factors (figure 2A), with 10.4% of the variance associated with immune cell proportions, while 30.7% was associated with HBsAg expression, Ishak score and/or unknown sample-specific differences (figure 2B,C). The first dimension was influenced by genes involved in fatty acid biosynthesis (SCD, FASN, THRSP, GPAM). The second dimension was associated with immunoglobulin genes (IGHG1/2/3/4, IGKC), human leucocyte antigen genes (HLA-A/B, HLA-DRA/DQA1/DQB1/DPA1) and T cell transcripts (TRBC1, CD8A) (figure 2D). The third and fourth dimensions only explained 5.4% and 4.4% of the global variance, respectively, and suggested ROI heterogeneity related to hepatocyte zonation (online supplemental figure 3A–C). The third dimension distinguished sample-specific phenotypes between PT1 and PT3 (online supplemental figure 3B). This observation did not associate with the quantity of immune cells or HBsAg+ cells. Instead, the top contributing genes included Kupffer/macrophage-associated transcripts (C1QC, C1QB, CD5L, HLA-DQA1) and the zonally expressed hepatocyte genes (CYP2A6, G0S2 and FASN) (online supplemental figure 3C). The fourth dimension differentiated a single ROI on the upregulation of canonical hepatocyte zone 3 markers (such as GLUL and CYP3A4), compared with the relative reduction in hepatocyte zone 1 transcripts such as SDS, A2M, LEPR, GLS2, TAT and HAL (online supplemental figure 3C).

Request permission

Differential gene expression identifies immune active profiles within patients

To profile the pathways in spatially heterogenous HBsAg expressing ROIs, in the presence or absence of co-infection, we performed differential expression analysis (figure 3A,B and online supplemental table 3). No differentially expressed genes (DEGs) were observed when comparing HBsAg+ areas with HBsAg− areas, suggesting that subtle viral-induced changes in the liver are not captured with the current number of samples. We therefore assessed DEGs without adjustment for multiple comparisons to identify potential gene candidates impacted by HBV infection. In PT1 (HBeAg−), high HBsAg+/CD45− staining areas showed an upregulation of ENKD1, AKT1S1 and ADIRF, involved in microtubule organisation, regulation of cell growth/survival in response to nutrient/hormonal signals and adipogenesis, respectively (figure 3A). PRDKC, a DNA damage sensor and effective activator of HBV cccDNA transcription,26 was also upregulated in these areas (figure 3A). The HBsAg+ areas in PT2 with HDV/HBV co-infection, were associated with immune genes IGHG1, IGHG2, IGHG3, IGHG4 participating in the ‘B cell receptor signalling pathway’. There were very few CD45+ cells visible by immunofluorescence, suggesting an underestimate of immune cell infiltration (figure 3B). No overlapping or conserved upregulated genes were found between PT1 and PT2. PT3 (HIV/HBV co-infection) displayed extensive membranous HBsAg staining with no HBsAg low regions for intratissue comparisons.

Request permission

Exploring gene enrichment associated with immune cell infiltrate, in PT1 the CD45+ regions showed an upregulation of genes associated with the constant region of immunoglobulin heavy chains (IGHG1, IGHG2, IGHG3, IGHG4); a common feature in the immune-rich areas of all samples tested (figure 4A). PT2 was characterised by additional enrichment of genes linked to viral infection and apoptotic process regulation (USP17L3, USP17L15, USP17L17, CASP1), along with cellular processing/major histocompatibility complex (MHC) class II presentation (CD74, HLA-DRB1) and chemokine genes (CXCL9, CCL19) involved in inflammatory responses/homing (figure 4B). In PT3 (HIV/HBV co-infection), immune-enriched areas exhibited increased expression of MHC class II genes (HLA-DRB1, HLA-DQB1), including classical and non-classical MHC class I genes (HLA-B, HLA-A, HLA-E, HLA-F) (figure 4C). Similar to PT2, the immune-high areas in PT3 showed increased expression of CXCL9 and CCL5, a major HIV suppressive factor, in addition to genes involved in CD8 T cell signalling (TRBC1, TRAC), T cell and natural killer (NK) cell cytotoxicity (NKG7) and interferon (IFN) transcriptional regulation and response (IRF1, GBP1). Notably, ‘cell cytotoxicity’ and ‘B cell receptor signalling pathway’ were common across all three patient tissues despite quantifiable differences in leucocyte presence (figure 4D–F and online supplemental table 4), indicating both shared and unique profiles.

Request permission

Distinct distribution of cellular phenotypes in the profiled ROIs

To assess the composition of the host’s immune response, we performed an unbiased analysis of cellular phenotypes present in the profiled ROIs, deconvoluting gene expression patterns based on a training matrix of single-cell sequencing data from the human liver cell atlas (figure 5A–C). Estimated proportions of different cell types showed contribution from hepatocytes and non-parenchymal cells with distinct immune cell distributions localising in high immune areas/HBsAg+ areas. PT1 (HBeAg−) showed lower adaptive and innate immune cell proportions consistent with an inactive disease phase (figure 5A). Differences in immune cell composition were more striking for PT2 with HDV/HBV co-infection (HBeAg+; high fibrosis) showing higher abundance of γδT cells, αβT cells, mature B cells, NK-like and inflammatory macrophages in addition to higher proportions of cholangiocytes and stellate cells relative to the low immune/low HBsAg+ areas (figure 5B). The immune-rich HBsAg+ areas in PT3 (HIV/HBV co-infection) were distinguished by higher abundance of αβT cells and inflammatory macrophages relative to PT2 (HDV/HBV co-infection), as well as the low immune areas within the same donor (figure 5C). Interestingly, the presence of antibody secreting cells was identified in PT2 and PT3 within the immune-rich/high HBsAg areas, although with low abundance.

Request permission

For unbiased assessment of diverse transcriptomic patterns across the tissues, we used a WGCNA that identified 8 modules of co-expressed genes (range 157–2280 genes per module). Modules were characterised by (i) enrichment in GO.BP terms, (ii) correlation with cell type abundance calculated through cell deconvolution and (iii) estimates of CD45+ and HBsAg+ cells from immunofluorescent staining. This analysis allowed us to establish module identity in all modules with the exception of the Black module. Each tissue displayed distinct transcriptional signatures and module-trait relationships (figure 6A,B, online supplemental figure 4A, B). PT1 (HBeAg−) was characterised by tissue-wide expression of the brown module relating to hepatocyte function and metabolism including fatty acid, cholesterol metabolic process, lipid catabolism and aerobic respiration. PT2 (HDV/HBV co-infection, HBeAg+; high fibrosis) showed tissue-wide expression of the turquoise module enriched for ribosomal transcripts. PT3 (HIV/HBV co-infection; HBeAg+) was distinguished by the red module associated with B cell mediated immunity and phagocytosis (figure 6A,B, online supplemental figure 4A, B). Notably, the red module was upregulated in areas with the highest CD45+ cell counts, distinguishing it from comparable areas in PT2, suggesting that B cell immunity and phagocytosis were more active processes in PT3 with HIV/HBV co-infection (figure 6A and online supplemental figure 4A). Transcripts associating with amino acid metabolism were enriched in the green module, shared by the HDV/HBV and HIV/HBV co-infected samples (PT2 and PT3) in HBsAg high areas with low immune infiltration (figure 6B and online supplemental figure 4A, B). The yellow module, which is associated with lymphocyte activation, B cell signalling and cytotoxicity, was also shared by PT2 and PT3. The presence of CD45+ within a given ROI was strongly associated with the yellow module, containing genes implicated in adhesion, cytokine production, αβ T cell activation, lymphocyte proliferation, cytotoxicity and responses to tumour necrosis factors and type II IFNs reflecting the presence of plasma cells, B cells, αβ and γδT cells, inflammatory macrophages and NK-like cells (figure 6B and online supplemental figure 4A, B).

Request permission

An increased hepatocyte abundance was associated with the pink, blue and brown modules (online supplemental figure 4B). The pink module that involved genes related to responses to metal, such as zinc or copper, was associated with PT1. The blue module, containing genes characteristic of hepatocyte metabolic and catabolic processes involving fatty acids, alcohol, carboxylic acids, cholesterol and amino acid catabolism, was associated with HBsAg expression and was over-represented in the HBV mono-infected and HIV/HBV samples (PT1 and PT3); however, the correlation is skewed by the high expression and lack of HBsAg heterogeneity within these samples, respectively. Interestingly, PT1 (HBeAg−, low fibrosis tissue) was associated with both the brown and blue modules, displaying comparable metabolic and catabolic pathways, whereas a negative module-trait relationship was observed with PT2 (HDV/HBV co-infection, HBeAg+; high fibrosis), suggesting transcriptional disparities in hepatocytes between samples (online supplemental figure 4A, B).