Exploring inflammatory signatures in arthritic joint biopsies with Spatial Transcriptomics

Lately it has become possible to analyze transcriptomic profiles in tissue sections with retained cellular context. We aimed to explore synovial biopsies from rheumatoid arthritis (RA) and spondyloarthritis (SpA) patients, using Spatial Transcriptomics (ST) as a proof of principle approach for unbiased mRNA studies at the site of inflammation in these chronic inflammatory diseases. Synovial tissue biopsies from affected joints were studied with ST. The transcriptome data was subjected to differential gene expression analysis (DEA), pathway analysis, immune cell type identification using Xcell analysis and validation with immunohistochemistry (IHC). The ST technology allows selective analyses on areas of interest, thus we analyzed morphologically distinct areas of mononuclear cell infiltrates. The top differentially expressed genes revealed an adaptive immune response profile and T-B cell interactions in RA, while in SpA, the profiles implicate functions associated with tissue repair. With spatially resolved gene expression data, overlaid on high-resolution histological images, we digitally portrayed pre-selected cell types in silico. The RA displayed an overrepresentation of central memory T cells, while in SpA effector memory T cells were most prominent. Consequently, ST allows for deeper understanding of cellular mechanisms and diversity in tissues from chronic inflammatory diseases.


Method optimization for RA and SpA tissue for Spatial Transcriptomics
The optimal conditions for tissue permeabilization and tissue removal for the Spatial Transcriptomics protocol were determined. This process uses cDNA synthesis with Cy3-labeled dCTPs to create a fluorescent footprint of cDNA, as previously described. 1 For RA and SpA tissues, the following parameters were selected: The tissue biopsies were sectioned at 7μm, fixed in 2% NBF (neutral buffered formalin) for 10min and stained with Hematoxylin for 7min. The prepermeabilization step was performed with a 20min incubation at 37°C in 14U of collagenase type I (Life Technologies, Paisley, UK) diluted in 1xHBSS buffer (Thermo Fisher Scientific, Life Technologies, Paisley, UK) supplemented with 1% BSA (NEB, Ipswich, MA, USA). For permeabilization, the sections were incubated in 0.1% pepsin-HCl (pH 1) for 8min at 37°C. For the tissue removal, the tissue was first incubated for 1h at 56°C in beta-mercaptoethanol which was diluted in RLT buffer (x4). After that, the tissue was treated with proteinase K diluted in PKD buffer at a ratio of 1:7 for 1h at 56°C, as previously described. The cDNA footprint was revealed by scanning the slide in an InnoScan 910 Microarray Scanner with 5μm resolution and gain at 10% or with a Metafer VSlide system at 20x resolution. An image showing the cDNA-footprint from one of the tissue optimizations can be seen in Fig. S2A-C. Figure S1. Overview of the patient samples and the spatial transcriptomics assay (A) Depiction of the joint biopsies (knee and hip) and brief patient characteristics. (B) Age, gender and disease duration at the timepoint of collection. Treatments at time of sampling include 1; DMARDs, 2; TNF inhibitors, 3; anti-IL-6R blockers, 4; glucocorticoids. (C) Depiction of the main steps of the spatial transcriptomics method starting with cryosectioning (three sections per biopsy), cDNA synthesis, tissue removal and release on surface cDNA, library preparation, DNA-sequencing and analysis. Figure S2. Spatial Transcriptomics optimization assay (A) An overview of the cDNA synthesis step for the ST optimization assay. The surface probes are distributed all over the surface area, and cDNA synthesis is conducted with Cy3-labeled dCTPs in the nucleotide mixture and thereafter scanned with a 532 nm light source. (B) Shows the same section scanned with the Metafer VSlide system (H&E image to the left) and scanned on 532 nm after the tissue section was removed (to the right). The cDNA footprint can be seen with no signal diffusion and also less signal in areas with less cells (shown with an arrow). The small white circles 1-4 is areas in which the signal intensity was measured seen in (C) which shows the average signal intensities measured with GenePix Pro 5.0 from the circular areas with error bars.

Figure S3. Cluster selection for RA-A 3
Black lines are regions annotated as lymphocyte infiltrates. Each spot represent areas in which mRNA was captured from. The numbers and colours represent the different clusters that were assigned. For the bulk analysis all spots were selected. For the infiltrate analysis cluster 4 was selected. For the outside-infiltrate analysis clusters 1 and 3 were selected.

Figure S4. Larger set of genes from the DE analysis on the bulk data.
The dendrograms were ordered so that all RA-samples are to the left and all SpAsamples to the right. For each selection, all genes that had a log2fold change >2 and p-value lower than 0.01 from the DE-analysis are shown. N.B. the reference genome GRCh38 only maps to the haplotype HLA-DR1, which contains HLA-DRB1 and HLA-DRB5 which is a homolog of HLA-DRB4. It is possible that faulty mapping can occur for specific variable genes like the ones building the MHC structures. The genes known to have this issue have been marked with *.

Figure S5. Larger set of genes from the DE analysis of selected clusters correlating to infiltrate regions.
The dendrograms were ordered so that all RA-samples are to the left and all SpAsamples to the right. For each selection, all genes that had a log2fold change >2 and p-value lower than 0.01 from the DE-analysis are shown.N.B. the reference genome GRCh38 only maps to the haplotype HLA-DR1, which contains HLA-DRB1 and HLA-DRB5 which is a homolog of HLA-DRB4. It is possible that faulty mapping can occur for specific variable genes like the ones building the MHC structures. The genes known to have this issue have been marked with *. The dendrograms were ordered so that all RA-samples are to the left and all SpAsamples to the right. For each selection, all genes that had a log2fold change >2 and p-value lower than 0.01 from the DE-analysis are shown. N.B. the reference genome GRCh38 only maps to the haplotype HLA-DR1, which contains HLA-DRB1 and HLA-DRB5 which is a homolog of HLA-DRB4. It is possible that faulty mapping can occur for specific variable genes like the ones building the MHC structures. The genes known to have this issue have been marked with *.

Figure S7. Top upstream regulators
Top upstream regulators for RA and SpA in the bulk RNA (A) and infiltrate areas (B).

Figure S8. Circos plot and heatmap of enriched terms outside infiltrates (A)
The Circos plot shows how genes (hits -RA, 304 genes; hits 2 -SPA, 113 genes) from the input gene list overlap (please see description in Figure Table S1. Top canonical pathways from IPA on bulk data. In the tables, similarities are marked in blue.