Identification and Characterisation of Infiltrating Immune Cells in Malignant Pleural Mesothelioma Using Spatial Transcriptomics

Increasing evidence strongly supports the key role of the tumour microenvironment in response to systemic therapy, particularly immune checkpoint inhibitors (ICIs). The tumour microenvironment is a complex tapestry of immune cells, some of which can suppress T-cell immunity to negatively impact ICI therapy. The immune component of the tumour microenvironment, although poorly understood, has the potential to reveal novel insights that can impact the efficacy and safety of ICI therapy. Successful identification and validation of these factors using cutting-edge spatial and single-cell technologies may enable the development of broad acting adjunct therapies as well as personalised cancer immunotherapies in the near future. In this paper we describe a protocol built upon Visium (10x Genomics) spatial transcriptomics to map and characterise the tumour-infiltrating immune microenvironment in malignant pleural mesothelioma. Using ImSig tumour-specific immune cell gene signatures and BayesSpace Bayesian statistical methodology, we were able to significantly improve immune cell identification and spatial resolution, respectively, improving our ability to analyse immune cell interactions within the tumour microenvironment.


Background
Malignant pleural mesothelioma (MPM) is a rare and aggressive tumour with less than 12 months median overall survival (OS) in advanced stages [1]. The incidence of MPM is 0.7 and 0.3 per 100,000 males and females, respectively, with a higher incidence in countries with greater exposure to asbestos such as Australia [2]. Most patients are diagnosed at an advanced stage, making them ineligible for curative surgery [1]. Consequently, systemic therapy has been a foundation for unresectable MPM [3] with the recently approved immune checkpoint inhibitors (ICIs) drastically increasing the OS by up to 18.4 months in MPM patients [4].
Increasing evidence strongly supports the key role of the tumour microenvironment on response to systemic therapy. The tumour microenvironment is a complex tapestry of immune cells [13], some of which possess immunosuppressive roles by impairing T-cell infiltration and antigen presentation or by enhancing immunosuppressive metabolism [14]. This heterogeneity may impact ICI therapy outcomes, but remains poorly understood in this context, possessing the opportunity to reveal specific risk factors impacting ICI therapy. In this paper, we describe a protocol for using spatial transcriptomics to map the tumour-infiltrating immune microenvironment. We describe step-by-step guidance on how to combine Visium spatial technology with unique bioinformatic approaches to overcome the limitations of current Visium spatial gene expression and to enable the identification of immune cell populations at significantly enhanced resolution.

Tissue Samples
Research sample collection and analysis was conducted as part of the NCT04631731 study (ICEMELT), which has been approved by the Western Sydney Local Health District (WSLHD) Human Research Ethics Committee (Ethics reference number: 2020/ETH02285). Malignant pleural mesothelioma patients were recruited at Blacktown and Westmead Public Hospitals in the WSLHD of Sydney, Australia. All patients provided written consent to participate in the study. A formalin-fixed paraffin-embedded (FFPE) tissue block was obtained from New South Wales Health Pathology (NSWHP), each containing tissue from the original tumour diagnostic biopsy prior to ICI therapy (n = 4). For the purpose of this pilot study, n = 4 FFPE blocks from metastatic MPM patients were selected. All patients were treated with IPI + NIVO and coded as A1, B1, C1 and D1. The quality of the FFPE blocks was first assessed followed by pathological annotation. Biopsies of patients A1, B1 and C1 were obtained from the original tumour site, and in D1, the biopsy was obtained from a metastatic site in temporal lobe. One section of 5 µm thickness was cut from each block and stained with Hematoxylin & Eosin (H&E) for pathological review. The clinical pathologist examined the tissue to identify an appropriate high yield as the size of the tissue which can be mounted on a Visium slide is only 6.5 × 6.5 mm. The pathologist microscopically selected a 6.5 × 6.5 mm area of interest for this study, identifying the type of tissue (tumour, non-tumour, adipose tissue, immune cell infiltrates). For the purpose of our study, the area with a significant accumulation of infiltrating lymphocytes was chosen. Tumour core and periphery represent unique immune cell niches with altered cell densities that must be considered if performing comparative analysis among patient subgroups. We will elaborate on this concept in the discussion. The pathologist's annotation was performed in one day and adjacent/tumour tissues within the examined area were characterised. Another 10 µm thick section from each block was cut for RNA extraction and quality control. Only samples with sufficient RNA quality, DV200 (percentage of RNA fragments of more than 200 nucleotides) > 50%, were used for this study.

RNA Quality Check
RNA quality was measured from 10 µm sections extracted using the RNeasy kit from Qiagen as recommended by Visium guidelines. Briefly, FFPE tissue sections were deparaffinised, treated with Proteinase K (60 µg/mL) followed by heat treatment (56 • C for 15 min followed by 80 • C for 15 min), DNase (16 µL of DNase booster buffer and 10 µL of DNase I stock solution [2700 U/mL]), and transferred to Rneasy MinElute column. Following 2 wash steps, the quality of the eluted RNA was assessed using the Agilent TapeStation to obtain a DV200 score for each FFPE tissue. In our experience, tissue source and length of storage can significantly impact RNA quality for the purposes of spatial transcriptomics. Gastrointestinal tissue FFPE blocks, for example, are subject to significantly faster RNA degradation than MPM blocks stored for a similar time frame. All four MPM FFPE blocks obtained a sufficient DV200 score that was necessary to proceed.

Visium Tissue Section Test Slide
Before placing sections on Visium slides, we first practiced mounting of the tissue sections onto 6.5 × 6.5 mm capture areas on regular Superforst Plus glass slides followed by staining with H&E to confirm adequate mounting and adhesion. We next followed the protocol provided by 10x Genomics using the Visium Tissue Section test slide to ensure adequate adhesion to the Visium slides. We observed very poor adhesion to these slides, with sections coming off the slides during H&E staining, even if processed delicately. We suspected that the absence of spatial barcode oligonucleotides impacted tissue adhesion, a complication that was confirmed by 10x Genomics representatives. We therefore decided to proceed with mounting our 4 sections onto the Visium Spatial Gene Expression slide. All sections adhered well to the Visium slides, suggesting that the bound oligonucleotides likely improve tissue adherence. These preliminary data impacted subsequent Visium guidelines, suggesting that users skip the tissue adhesion test for future experiments.

Histology and Imaging
The experimental protocol was scheduled to take 4 days (Table 1) beginning with a sectioning of FFPE blocks at a thickness of 5µm (Leica Biosystems Microtome). We strictly followed the protocol provided by 10x Genomics and no deviations were made. Each tissue section was placed within one of four 6.5 × 6.5 mm capture regions of the Visium slide. Each of the capture area contains an array of~5000 spots (55 µm in diameter) which have oligonucleotide sequences required for the capture of probe-ligated mRNA, a spatial barcode, unique molecular identifier (UMI) and Illumina sequencing primer binding site ( Figure 1). Sections were left overnight to dry and incubated the following day in a 60-degree fan-forced oven for 2 h and stained with H&E. The protocol for H&E staining was provided by 10x Genomics and was strictly followed. All reagents were newly purchased for the purpose of this analysis and were RNase-free. H&E-stained slides were coated with 85% glycerol and a coverslip, and imaged with the Olympus VS120 Slide Scanner at 40× lens magnification. After imaging, glycerol was washed off by immersing the Visium slide into the beaker containing 800 mL of Milli-Q water. The slide was immersed horizontally and held until the coverslip fully separated. Sections were dried at 37 • C and stored at 4 • C before decrosslinking. Table 1. A step-by-step plan of the Visium workflow. D-day; FFPE-formalin-fixed paraffinembedded; H&E-hematoxylin and eosin. Step

Sub-Steps Location Estimated Timing
Case selection and quality control  Overview of Visium spatial analysis workflow. Sections from FFPE tissues chosen for the spatial transcriptomic workflow are first quality checked, then freshly cut sections are adhered to the Visium slide. Following deparaffinisation and staining, they are imaged using a slide scanner followed by library construction, PCR amplification and sequencing. Spatial data are finally analysed using the Loupe Browser or other bioinformatic techniques as we describe herein. H&E-Hematoxylin & Eosin; UMI-unique molecular identifier.

Library Preparation and Sequencing
Samples were decrosslinked using the recommended TE buffer (pH 9.0) before human transcriptome probes covering 18,000 genes were ligated to target mRNA for capture and library construction according to the Visium Spatial Gene Expression for the FFPE protocol with no deviation. For PCR indexing, the optimal PCR cycle number determination was performed for each sample with the following cycles used in the final indexing PCR reaction: A1 (19 cycles), B1 (19 cycles), C1 (21 cycles) and D1 (18 cycles). DNA libraries were assessed using High Sensitivity D5000 ScreenTape on the TapeStation platform (Agilent Technologies), followed by sequencing using on a Novaseq SP flow cell with 50 bp paired end read length. The minimum sequencing depth required for each sample was calculated using the manual alignment tool in Loupe Browser to estimate the number of tissue-covered spots in each capture area and multiplying the spot number by the recommended 25,000 reads. On this basis, the minimum read pairs required for each sample were A1 (58.9 M reads), B1, (31.1 M reads), C1 (52.2 M reads) and D1 (90.4 M reads). Approximately 100 M paired reads were acquired for each sample and Fastq files were generated using the bcl2fastq 2.20.0.422 pipeline. Overview of Visium spatial analysis workflow. Sections from FFPE tissues chosen for the spatial transcriptomic workflow are first quality checked, then freshly cut sections are adhered to the Visium slide. Following deparaffinisation and staining, they are imaged using a slide scanner followed by library construction, PCR amplification and sequencing. Spatial data are finally analysed using the Loupe Browser or other bioinformatic techniques as we describe herein. H&E-Hematoxylin & Eosin; UMI-unique molecular identifier.

Library Preparation and Sequencing
Samples were decrosslinked using the recommended TE buffer (pH 9.0) before human transcriptome probes covering 18,000 genes were ligated to target mRNA for capture and library construction according to the Visium Spatial Gene Expression for the FFPE protocol with no deviation. For PCR indexing, the optimal PCR cycle number determination was performed for each sample with the following cycles used in the final indexing PCR reaction: A1 (19 cycles), B1 (19 cycles), C1 (21 cycles) and D1 (18 cycles). DNA libraries were assessed using High Sensitivity D5000 ScreenTape on the TapeStation platform (Agilent Technologies), followed by sequencing using on a Novaseq SP flow cell with 50 bp paired end read length. The minimum sequencing depth required for each sample was calculated using the manual alignment tool in Loupe Browser to estimate the number of tissuecovered spots in each capture area and multiplying the spot number by the recommended 25,000 reads. On this basis, the minimum read pairs required for each sample were A1 (58.9 M reads), B1, (31.1 M reads), C1 (52.2 M reads) and D1 (90.4 M reads). Approximately 100 M paired reads were acquired for each sample and Fastq files were generated using the bcl2fastq 2.20.0.422 pipeline.

Bioinformatic Analysis
Gene expression counts were generated from fastq files using space ranger V1.3 with refdata-gex-GRCh38-2020-A and Visium Human Transcriptome Probe Set version 1.0 GRCh38-2020-A.csv annotations. Loupe files (generated from raw data using 10x Genomics Cloud Analysis software) were analysed using loupe browser and count matrices were imported into R for analysis with Seurat [15]. Clustering and spatial enhancement was performed using BayesSpace according to the author's instruction [16]. Briefly, optimised clustering resolution was determined using the qtune function (between 5-7 clusters per section), and the optimised enhancement parameter jitter_scale was determined using the mcmcChain function on low replicate (n = 300) spatial-enhanced data; for our sections, the optimum scale was 0.35. Spatial enhancement was performed using these optimised parameters for 10,000 repetitions.
To define the tumour-infiltrating immune cells within the spatial framework, gene expression signature lists were used from the R package ImSig [17] for immune cells (B cells, T cells, NK cells, neutrophils, macrophages and monocytes). To define MPM tumour tissue, a gene list containing 51 overexpressed genes identified in MPM was used (e.g., MSLN, ALDOA, CDC2, NMU, PCNA, PDGFRB) from a systematic review of 10 independent transcriptome studies [18]. Enhanced feature expression was generated for all nominated genes (n = 401) using BayesSpace, and for each feature, a score was generated using the sum of scaled expression normalised to the size of the gene list and overall gene expression for that feature. Cell types were then assigned to features by scores greater than the background and verified by manual image analysis. A manual cut-off of 250 was used on the scores enhanced by the expression scores that were scaled, using the sum of the signature gene scaled expression values divided by the number of genes in each individual signature [19]. Neighbour counting was carried out using a custom script in R. Briefly, for each feature, all features sharing at least one vertex (up to 12 per feature) were evaluated for cell type assignments and tallied for each. In addition, Loupe browser version 6.0.0 (10x Genomics) was used for analysis of the spatial data. ConsensusPath-DB online software was used to functionally annotate the set of significant (p < 0.05) DEGs. Graphical summaries of immune cell compositions were completed using Prism version 9 (GraphPad Software, LLC). Statistics were not performed on this data due to insufficient power (n = 4 sections).

Results
Tissue sections (6.5 mm × 6.5 mm) from each block were processed to quantify spatial gene expression (VISIUM, 10x Genomics) offering 9041 barcoded arrays across four samples with an average of 3125 genes detected per spot. The summary reports from Space Ranger (10x Genomics) are presented in Supplementary Materials ( Figure S1). Figure 2A,B,E,F demonstrate H&E staining as well as the pathologist's annotation of tissue sections A1 and B1, demonstrating the heterogeneity of tissue as well as immune cell aggregation in tumour tissue. The goal of this study was to identify and transcriptionally characterise immune cells within the tumour landscape. Unfortunately, this was not immediately feasible for two reasons: low spatial resolution and depth of sequencing. The 55 µm spot diameter present on the Visium slides contained transcriptional data from multiple cells, resulting in mixedcell transcriptomes. This issue was particularly problematic when identifying immune cell clusters, representing densely packed areas of multiple immune cell populations. This was evident in tissues A1 and B1, where immune cells aggregated at the tumour periphery. Secondly, the depth of sequencing was not sufficient to identify immune cells based on the expression of a single gene (e.g., CD3+ T cells, CD19+ B cells or CD56+ NK cells). Moreover, different immune cell-specific transcripts would generate significantly different spatial maps of immune cell location (Figure 2C To overcome the limitation of immune cell identification, we used ImSig tumour-specific immune cell gene signatures [17] that enabled us to confidently identify tumour infiltrating immune cell populations (T cells, B cells, monocytes, macrophages, neutrophils, natural killer [NK] cells) (Figure 3). In addition to immune cells, ImSig enabled the identification of other cell signatures including interferon, proliferation and translation that are intricately linked to oncogenic signalling and anti-cancer responses [17]. Applied to the spatially resolved data, ImSig, immune cell transcriptomes and locations could be examined and compared within sections and between sections based on a variety of clinical parameters including treatment response and/or ICI toxicity. Due to limitations of the Visium slide, some spots were assigned to multiple immune cell types. We manually labelled each of the spots sequentially, beginning with T cells, followed by B cells, NK cells, macrophages, monocytes and neutrophils ( Figure 3G,H). To overcome the limitation of immune cell identification, we used ImSig tumourspecific immune cell gene signatures [17] that enabled us to confidently identify tumour infiltrating immune cell populations (T cells, B cells, monocytes, macrophages, neutrophils, natural killer [NK] cells) (Figure 3). In addition to immune cells, ImSig enabled the identification of other cell signatures including interferon, proliferation and translation that are intricately linked to oncogenic signalling and anti-cancer responses [17]. Applied to the spatially resolved data, ImSig, immune cell transcriptomes and locations could be examined and compared within sections and between sections based on a variety of clinical parameters including treatment response and/or ICI toxicity. Due to limitations of the Visium slide, some spots were assigned to multiple immune cell types. We manually labelled each of the spots sequentially, beginning with T cells, followed by B cells, NK cells, macrophages, monocytes and neutrophils ( Figure 3G,H). To address the spatial limitation of the Visium method, we used the BayesSpace package that utilises the Bayesian statistical method to subdivide each spatial spot and obtaining subspot resolution. BayesSpace uses transcriptional data from neighbouring spots to increase the resolution of each 55 µm spot 6-fold. As immune cell sizes are in the range of~5-20 µm and parenchymal cells are even larger, this method increased the resolution of our sections to a nearly single-cell level. Enhancement of data allowed a more precise identification of immune cell clusters within the tumour (Figure 4). BayesSpace resolved the tissue structure that is not detectable at original resolution and applicable both to infiltrating immune cells and malignant cells ( Figure S2), thus overcoming a major limitation of the current Visium spatial technology. To address the spatial limitation of the Visium method, we used the BayesSpace package that utilises the Bayesian statistical method to subdivide each spatial spot and obtaining sub-spot resolution. BayesSpace uses transcriptional data from neighbouring spots to increase the resolution of each 55 µm spot 6-fold. As immune cell sizes are in the range of ~5-20 µm and parenchymal cells are even larger, this method increased the resolution of our sections to a nearly single-cell level. Enhancement of data allowed a more precise identification of immune cell clusters within the tumour (Figure 4). BayesSpace resolved the tissue structure that is not detectable at original resolution and applicable both to infiltrating immune cells and malignant cells ( Figure S2), thus overcoming a major limitation of the current Visium spatial technology. We next sought to analyse immune cell makeup, phenotype and interactions in normal and tumour tissue using BayesSpace-enhanced images, with which we have generated coordinate-based information for all ImSig-based annotation. Figure 5A demonstrates the immune cell makeup of tumour and non-tumour tissue based on the ImSig identification of immune cells. Both A1 and B1 sections demonstrate drastically different immune cell compositions, dominated by myeloid cells in the tumour and lymphocytes in adjacent tissue. Conversely, the immune cell composition within C1 and D1 are more similar within tumour and adjacent tissue. Tumour and adjacent immune cell populations were next characterised based on whether they shared proliferation or interferon signatures defined using ImSig ( Figure 5B). The proportion of immune cell spots sharing the proliferation signature was largely increased in the tumour tissue, whereas the interferon signature was generally enhanced in the adjacent non-tumour tissue immune cells. Lastly, immune cell interactions were quantified by calculating the average number of NK and T-cell neighbouring spots held by individual immune cell populations, with a maximum of 13 neighbours (including co-localisation) following BayesSpace enhancement ( Figure 5C). Both T-and NK-cell spots generally demonstrated more interactions with myeloid cell populations (monocytes, macrophages and neutrophils) in tumour tissue and fewer interactions with other lymphocyte populations (B, T and NK cells). These data provide crucial information regarding immune cell interactions with potential clinical implications; however, the limited sample size and significant variability among sections used in this study have limited any further statistical considerations.  . (A,B) Enhanced resolution allowed us to more precisely locate and quantify tumour-infiltrating immune cells, particularly T cells, B cells, NK cells and macrophages using the ImSig algorithm. (C,D) Scales and histograms represent the cut-off values assigned to elucidate the signature scores for each section and/or cell type using an outlier analysis (hampel filter) [20].
We next sought to analyse immune cell makeup, phenotype and interactions in normal and tumour tissue using BayesSpace-enhanced images, with which we have generated coordinate-based information for all ImSig-based annotation. Figure 5A demonstrates the immune cell makeup of tumour and non-tumour tissue based on the ImSig identification of immune cells. Both A1 and B1 sections demonstrate drastically different immune cell compositions, dominated by myeloid cells in the tumour and lymphocytes in adjacent tissue. Conversely, the immune cell composition within C1 and D1 are more similar within tumour and adjacent tissue. Tumour and adjacent immune cell populations were next characterised based on whether they shared proliferation or interferon signatures defined using ImSig ( Figure 5B). The proportion of immune cell spots sharing the proliferation signature was largely increased in the tumour tissue, whereas the interferon signature was generally enhanced in the adjacent non-tumour tissue immune cells. Lastly, immune cell interactions were quantified by calculating the average number of NK and T-cell neighbouring spots held by individual immune cell populations, with a maximum of 13 neighbours (including co-localisation) following BayesSpace enhancement ( Figure  5C). Both T-and NK-cell spots generally demonstrated more interactions with myeloid cell populations (monocytes, macrophages and neutrophils) in tumour tissue and fewer interactions with other lymphocyte populations (B, T and NK cells). These data provide crucial information regarding immune cell interactions with potential clinical implications; however, the limited sample size and significant variability among sections used in this study have limited any further statistical considerations.

Discussion
Spatial transcriptomics has recently been declared the Method of Year by the Journal Nature Methods [21]. Spatially resolved transcriptomics enables the characterisation of transcriptional patterns within sections of tissue while preserving the original tissue architecture [22]. It overcomes the limitations of bulk and single-cell sequencing methods which can provide the transcriptomic data without tissue alignment. A large study (516 MPM samples) conducted by Aday et al., for example, established that T H2 and cytotoxic T cells are abundant in tumours of patients with greater overall survival [23]. Their study was nonetheless unable to identify cellular identity and phenotype, its relation to neighbouring cells and non-cellular structures which provides crucial information regarding tumour composition and the role of infiltrating lymphocytes in treatment response [24]. There are a few commercialized spatial techniques that can be used on both FFPE and fresh-frozen tissues [22]. Herein, we describe a methodology of using a Visium spatial transcriptomics method (10x Genomics) on FFPE tissue samples collected from a real-world cohort of cancer (MPM) patients. In this pilot study we were able to successfully implement a Visium-based framework to identify and characterise tumour-infiltrating immune cells in FFPE tissues using various packages in R.
While Visium spatial transcriptomics is a powerful technology, there are limitations of the current slides that must be noted. Most importantly, the diameter of each spatial spot under the tissue section is 55 µm, which will likely capture transcripts from multiple cells, thus diluting important information regarding cellular identity and phenotype. BayesSpace was able to increase resolution 6-fold; however the spatial enhancement comes at the cost of greatly increasing the amount of data produced and a greater reliance on orthogonal validation. Most importantly, the resolved images could not be exported to the Loupe browser and aligned with gene expression data. Consequently, differential gene expression could not be performed. Nonetheless, the next generation of spatial kits offered by 10x Genomics provide better resolution (much smaller spot size) with a larger tissue capture area [25]. This will surely simplify downstream analysis; however, cell identification tools such as ImSig will likely still be necessary to identify cell populations.
The second major limitation of spatial transcriptomics technology is the absence of open-source and well-documented software for performing data analysis. A creation of a unique bioinformatic tool for data collection and analysis will markedly simplify the workflow and allow researchers to focus more on understanding the biological problem rather than spending time learning new syntax or bringing in seasoned bioinformaticians to implement a new algorithm/code using R-based frameworks. In the near future, more advanced tools and techniques will likely appear to provide more in-depth analysis of spatial data, thus dramatically changing the field of clinical and biomedical research. Moreover, spatially resolved transcriptomics have a potential to be introduced into clinical diagnostic guidelines as a high-throughput tool for guiding physicians' decision making.
Tissue acquisition and selection for spatial profiling represents an additional hurdle that can significantly impair analyses as performed herein. Extensive planning and communication with medical staff is necessary to obtain tumour tissue that is well characterised histologically with regard to location within the tumour (central versus peripheral). The retrospective acquisition of FFPE biopsy tissue represents a significant limitation of this study. Biopsy tissue was obtained from different locations within the tumour, and hence possesses different immune cell profiles. Sections A1 and B1 are clearly isolated from a peripheral tumour location, as demonstrated by adipose tissue surrounding the tumour. Indeed, these sections demonstrate significantly different immune cell profiles and interactions as compared to sections C1 and D1 that are more homogenous with regard to tissue composition. Future studies sampling homogenous tumour locations will better inform how clinical parameters (e.g., treatment response, toxicity) are affected by the tumour microenvironment.
In summary, the preservation of tissue architecture and alignment with gene expression data allowed us to determine number and location of tumour-infiltrating immune cells and the upregulation of genes potentially responsible for better/lower therapeutic outcomes in examined patients. Further exploration of cell-cell interactions by implementing spatial and single-cell methods is necessary, and in the near future, it may identify the unique cellular markers correlated with immunotherapy outcomes in patients with rare cancers.

Institutional Review Board Statement:
Research sample collection and analysis was conducted as part of the NCT04631731 study (ICEMELT), which has been approved by the Western Sydney Local Health District (WSLHD) Human Research Ethics Committee (ethics reference number: 2020/ETH02285).

Informed Consent Statement:
Each recruited patient has provided written informed consent for participating in ICEMELT (NCT04631731) study.

Data Availability Statement:
The data are available upon a request to a corresponding author.