Single-cell transcriptomics of clinical grade adipose-derived regenerative cells reveals consistency between donors independent of gender and BMI

Abstract

Adipose-derived regenerative cells (ADRCs) also referred to as the stromal vascular fraction, provide an ample source of stem cells with widespread regenerative therapeutic use. Being heterogenous in nature, possibly affecting the clinical outcome after stem cell treatment, the ADRC- donor, -BMI, and -gender may have a large impact on ADRC composition and quality but this remains largely unexplored. Herein, we provide a comprehensive single-cell RNA sequencing ADRC mapping across two cell trial intervention studies but found no gender- or BMI-related variations, except for a minor female increase in PI16/CD55-expressing stem cells. Indeed, ADRC heterogeneity was surprisingly minimal between donors. This provides important decision-making support on adipose stem cell donor selection for stem cell treatments, and suggest that donor, gender and BMI should be regarded as less influential.

Introduction

Subcutaneous adipose tissue (SAT) remains an inexhaustible source of stem cells that are used for various stem cell interventions [1]. Adipose-derived regenerative cells (ADRCs, a mixture of regenerative cells also referred to as stromal vascular fraction (SVF)) remain the cellular origin of SAT homeostasis and -expansion, and may be isolated for immediately use in autologous cell transplantation [2,3,4]. Alternatively, ADRCs can be expanded in culture to obtain adipose-derived stromal/stem cells (ASCs) before stem cell transplantations into patients in an autologous or allogeneic design and currently, more than 350 trials using adipose tissue-derived cells are registered as either completed, or recruiting at Clinical.trials.gov. However, since the ADRCs encompass all the cell types present in adipose tissue except mature adipocytes, the specific cell composition and heterogeneity may contribute significantly to the clinical outcome after stem cell treatments [1, 5,6,7]. For therapeutic use, several attempts have been utilized to obtain a more homogenous population of the stem cells such as adherence [5] and pre-/deselection of ADRC subpopulations based on e.g. their expression of membrane proteins (e.g. PECAM1-positive endothelial- [8] or Integrin α10β1-selected mesenchymal stem cells [9]). To define ADRC or ASC heterogeneity, Kostecka and colleagues [1], most recently reviewed ADRC single-cell RNA sequencing (scRNAseq) studies and reported that fat deposit, donor-age and health status such as obesity, Type 2 diabetes (T2D) etc., as well as ADRC isolation and expansion protocols may impact heterogeneity and also clinical stem cell intervention outcome. Yet, studies in the field of metabolomics [10] have suggested that ADRC cell composition indeed is similar between fat deposits, despite that certain CD55/PI16-progenitors may correlate with the metabolic phenotype [10]. Such discrepancies may be explained by the small scRNAseq sample- and cell size [1, 10]. However, another caveat to these previous studies is the numerous approaches used for isolating ADRCs [10], procedures that are generally not compatible with clinical stem cell production. Thus, despite some data exist, ADRC heterogeneity between donors remains relatively uncertain. Moreover, it is surprising that no previous investigations focus on gender as a component underlying ADRC heterogeneity. To stem cell therapy, donor cell diversity and its relation to gender is a very important issue when selecting the most proper ADRC donor to benefit clinical outcome.

Herein we therefore performed scRNAseq of ADRCs isolated by the same Good manufacturing practice (GMP) compatible and reproducible approach across samples from two clinical stem cell interventions. With this, we investigated overall and individual cell composition heterogeneity, with a subsequent focus on variations in cell composition between men and women to aid in donor selection.

Results and discussion

GMP compatible ADRCs for cell therapeutics reveal high quality scRNAseq data

ADRCs were isolated from SAT obtained in relation to two phase II randomized, controlled trials (RCTs), conducted at Odense University Hospital, Denmark. The RCTs were designed to test the efficacy of autologous non-cultured ADRC therapy versus placebo (randomization ratio 1:1) in relation to treatment of breast cancer-related secondary lymphedema (BCRL, n = 80, NCT03776721 [11]) and post-radical prostatectomy related erectile dysfunction (pRP-ED, n = 70, EudraCT 2015–005140-33), respectively. From each trial, samples from 7 individuals were randomly selected (BCRL, females, cell therapy n = 7 and pRP-ED, males, cell therapy n = 4; placebo: n = 3) upon availability, for further scRNAseq processing. As previously exploited by us and others [2, 11,12,13,14] in several cell therapy trials, the CE-marked, GMP-compliant device (Cytori Celution IV) cell isolation platform was used for isolation of the ADRCs. Cell numbers and -sizes were similar between donors (data not shown) and in agreement with our previous studies [2, 11,12,13,14]. Immediately after isolation, we applied the ADRCs to a recently in-house validated methanol-fixation and freezing step [8, 15, 16] based scRNAseq protocol. With this, we facilitated high RNA integrity of the ADRCs (Supplementary Table 1, Additional File 1) while preserving in situ conditions as closely as possible. High-throughput scRNA-seq using 10X Genomics combined with NovaSeq was followed by quality controlling (Supplementary Table 1, Additional File 1) and data integration. Resulting quality-filtered, single-cell data (Supplementary Table 1, Additional File 1) from the 14 patients included a total of 111,983 human ADRCs that were integrated. This number of high quality GMP compatible ADRCs exceeded previous scRNAseq datasets from trials, thus providing an opportunity to investigate clinically relevant ADRC heterogeneity in detail between donors and genders.

ADRCs enclose a mixture of regenerative vascular cell types

Prior to donor and gender comparison, we initially performed unbiased scRNAseq clustering based on ADRC expression of the canonical markers PECAM1, PTPRC, TAGLN and CFD, respectively (Supplementary Figs. 1–4, Additional File 2). As such, we identified 24 cell clusters (Fig. 1A–E and Supplementary Table 2, Additional File 1) falling into four major well-separated subgroups of endothelial-, immune-perivascular mural- and adipose stem/progenitor cells. In line with Massier and colleagues [10], as well as Gene Ontology enrichment analyses (Supplementary Fig. 2–4, Additional File 2 and Supplementary Table 3, Additional File 1 using ShinyGO [17]) and selected cell type annotation tools (Supplementary Table 4, Additional File 1: CellMarker2024, Tabula Sapiens [18], Tabula Muris [19], using Enrichr [20]) we annotated the ADRC clusters.

Fig. 1

Single-cell RNA sequencing landscape of ADRCs integrated from 14 patients. A Left; Uniform manifold approximation and projection (UMAP) 2D visualization of 111,983 ADRCs in 24 clusters falling into four major and one minor cell categories: endothelial- (EC), immune- (IC), perivascular mural- (PC), adipose stem/progenitor- (AC) and red blood (RBC) cells. Right; Dendrogram showing the hierarchical relationship between the 24 ADRC clusters. B–E DimPlot UMAP with subcluster annotation,—fractions (%) and violin plot of selected genes for each of the major cell types i.e. B endothelial- (EC1-4), C perivascular mural- (PC1-3), D immune- (IC1-9), and E adipose stem/progenitor (AC1-7) cells. MyoFCs: myofibrocytes, vSMC: vascular smooth muscle cells, Mϕ: Macrophages, Mo/DC: Monocyte/dendritic cells, LAM: lipid associated macrophages, NKT: natural killer T-cells, APC: adipose precursor cells, CPA: committed preadipocytes, FAP: fibroblast and adipogenic progenitor cells

A major subgroup of ADRCs includes endothelial cells (EC), which are considered to hold the high angiogenic potential of ADRCs [8]. The top-three genes separating ECs from the remaining major subgroups were VWF [21], CLDN5 [22] and ADGRL4 [23, 24] (Fig. 1B and Supplementary Fig. 1C, Additional File 2). ECs could however, be further branched into 4 subsets classified according to known endothelial cell types along the arteriovenous axis of the microvasculature (Fig. 1B and Supplementary Fig. 2A, Additional File 2). The largest EC clusters, EC1 and EC2 expressed arterial- (SEMA3G, HEY1, EFNB2), and to a greater extent, capillary (RGCC, CA4, CD36, FABP4, GPIHBP1) marker genes. PLVAP (endothelial cell-specific plasmalemma vesicle-associated protein) [25] was expressed by all EC-clusters although at higher levels in the minor EC clusters, EC3 and EC4. PLVAP thus marks a transition to post-capillary venules represented by these two subsets that were also characterized by their expression of additional postcapillary venule-associated genes i.e. SELE, ACKR1, ICAM1 and VCAM1. The EC progenitor markers CD34 and ENG, appeared to be uniformly expressed across all 4 EC clusters although the expression pattern of angiogenic and tip-cell EC markers suggests that the majority of potential tissue-resident endothelial progenitors [26] are found within the EC1 and EC2 subsets. The arterial/capillary/venous EC features were similar to those reported by others [10] however, as previously with CD31/PECAM1-pre-selected ADRCs from SAT [8], we failed to identify any LYVE1/PROX1 expressing ECs. Notably the gene signature of our SAT-originating ECs was very similar to ECs from visceral adipose tissue (VAT) [27] (Supplementary Fig. 4B, Additional File 2).

Besides endothelial cells, perivascular mural cells (PC) exhibit vast regenerative potential [28]. Among all ADRCs in our dataset, the PC subgroup was characterized by high expression of RGS5 [29, 30], NDUFA4L2 [31] and ACTA2 (Fig. 1C and Supplementary Fig. 1C, Additional File 2) and encompassed 3 subclusters, PC1-3 with PC1 and PC2 being closest neighbors (Fig. 1A,C). PC2 consists of pericytes defined by the genes PDGFRB, CSPG4 (NG2), KCNJ8, ABCC9, NOTCH3 and RGS5 [31] and as for cluster EC1 and EC2, FABP4, CD36 and PPARG transcripts are abundant suggesting involvement in fatty acid transport across vessels as previously described in human skin [32]. Of the mural cells, PC2 expresses the highest levels of NDUFA4L2 and COX4I2, both encoding subunits of mitochondrial proteins reported to mark human brain pericytes and vascular smooth muscle cells (vSMC) [33] and which decline along the arteriole-arterial smooth muscle cell axis [33]. Thus, the lower NDUFA4L2- and COX4I2-expression in the mural cluster PC1 along with higher levels of genes involved in vSMC contraction including ACTA2, TAGLN, MYH11, MYLK, and CNN1 [31], suggest that PC1 contains cells representing vSMCs. LMOD1- and SYNPO2 have both previously been associated with vSMC maturity [34] and while SYNPO2 is slightly more pronounced in PC3 than in PC2, PC3 can however, be further distinguished by high levels of the chemokines CCL21, CCL19 and CCL2 (top-3 genes) identifying this cluster as myofibrocytes [35] (Fig. 1C).

Altogether, these data for EC and PC subtypes underscores that ADRCs in its fresh form contain a very defined vascular potential, which may be exploited when replenishing ischemic tissues upon ADRC treatment. It will be interesting in the future to further delineate which of these EC and PC clusters that embrace the highest regenerative potential or whether they work in synergy to accomplish revascularization of ischemic tissues.

Stem- and immunomodulatory cells constitute substantial fractions of ADRCs

A large number of studies have reported an immunomodulatory mode of action underlying SVF/ADRC function in stem cell treatments [36]. This capacity is thought to be enclosed mainly within the mesenchymal cell-like population (below) but also the ADRC-derived immune cells (IC) seem involved in this function [37, 38]. In our data, the transmembrane Immune Signaling Adaptor TYROBP (TYROBP), CD45 (PTPRC) and Serglycin (SRGN) separated the IC subgroup from the remaining ADRCs by being the top-3 highest expressed genes (Supplementary Fig. 1C, Additional File 2). The ICs further branched into 9 clusters of immune cells expressing markers of various known immune cells of both myeloid and lymphoid lineages (Fig. 1D and Supplementary Fig. 3, Additional File 2). Most of the identified ICs are of the myeloid lineage (20% of all ADRC, 71% of all ICs) with further specification into monocyte/dendritic/macrophage (IC1,-2,-4) subtypes as well as a neutrophil-like cluster (IC3), mast cells (IC5) and a small cluster of contaminating red blood cells (RBC). A mixed cell cluster (IC9) appeared related to the RBC cluster, both being involved in oxygen transport (GO-biological function, Supplementary Fig. 3, Additional File 2 and Supplementary Table 3–4, Additional File 1).

Finally, three clusters belonging to the lymphoid lineage, encompassed T/NKT-cells (IC 6 + IC7) and B-cells/plasma cells (IC8). Of note, the frequency distribution and absolute cell counts of circulating immune cells will depend on initial blood infiltration and following washing steps (if any) and may, from a clinical perspective be important [39]. The differences in number- and subtypes of infiltrating immune cells can simply also reflect an actual disease state like in the case of diabetic HIV patients [40]. It is however still speculative, which IC subtypes that facilitates the immunomodulatory action, but both Tregs and monocytes have been suggested to enclose this activity [37, 38]. Besides ICs, also the mesenchymal stem cells in ADRCs possess immunomodulatory action in addition to other functions such as regeneration and wound healing [41]. Our data are in agreement with other studies [10] and showed that a rather large fraction (43%) of ADRCs were classified as adipose stem/progenitor cells (AC) distinguished herein from remaining ADRCs by their high expression of C1S, MFAP4 and NOVA1 (Fig. 1e and Supplementary Fig. 1C, Additional file 2). While Complement C1S, C3, CFD and Microfibril Associated Protein 4, MFAP4 mRNAs are known to mark adipogenic/fibroblast progenitors, the expression of NOVA1 [42, 43] has only very recently been highlighted in human adipose tissue and the stromal vascular fraction (SVF, equal to ADRCs) hereof [42]. Interestingly, Yang and colleagues recognized NOVA1 as a regulator of chromatin-remodeling during- and a prerequisite for adipogenesis in humans [42]. This function appeared to be species-specific and could not be recapitulated in mice, which can possibly add to explain the lack of NOVA1-observations in adipose tissue-derived sc/snRNAseq data that are dominated by mouse data. GEO Profile data (NCBI, GEO Profile ID: 111023595; DataSet ID: 5056 [44] and GEO Profile ID: 83214995; DataSet ID: 4276 [45]) indicate NOVA1 expression to be more pronounced in subcutaneous than visceral fat and that NOVA1 mRNA levels may be higher in human adipose tissue or -derived ASCs from lean as compared to obese individuals. From a large set of integrated human sc- and snRNAseq data [10] NOVA1 was noted in one of the minor omental- but not subcutaneous fat clusters. Herein, we reinspected their original data, and found that NOVA1 indeed was the second highest ranking gene (fold-change/adj. p-value ranking) in one of the largest subcutaneous fat clusters (sfC01 [10]). Our AC class of cells further segregated into 6 common (AC1-6) clusters while one (AC7) was specific to a single individual (male, #M2. Fig. 2C, D and Supplementary Fig. 1A, Additional File 2). Subclusters AC1 and AC3-5 appeared closely related as did subclusters AC2, AC6 and AC7 (Fig. 1A) and all are enriched for Gene Ontology (GO) terms and/or genes related to ECM/fibroblast/stromal cells (Supplementary Fig. 4A, Additional File 2 and Supplementary Table 3–4, Additional File 1). As noted by several groups [10, 46], the nomenclature for the adipose progenitor/stem cell class (termed fibroblast and adipogenic progenitors (FAPs)) is not well-established possibly due to heterogeneity between species and fat deposits. Despite of this heterogeneity, Massier et al. [10] reported the VAT and SAT clusters from their large integrated sc/snRNAseq object to be very similar. Also our SAT derived AC clusters recapitulate the expression of genes used for clustering visceral fat FAPs [27] including BMPER, recently recognized as a positive modulator of adipogenesis and highly expressed in VAT [27] (Supplementary Fig. 4B, Additional File 2). Jaccard similarity scores between our AC clusters and the described 17 SAT-derived FAPs [10] (Supplementary Fig. 4C, Additional File 2), suggests that clusters AC2, AC6 and to a less degree AC7 represent a central APC cluster, marked by the genes PI16/CD55 characteristic of non-committed APCs. AC1 and AC4-5 confine cells moving towards a committed preadipocyte (CPA) phenotype, while AC3 appears to represent terminally differentiated fibroblast/stromal cells likely arising from a different cell type other than PI16/CD55 stem cells [10]. Finally, we did a comprehensive comparison of our AC data to the AC-marker lists from other substantial human scRNAseq adipose tissue studies (Supplementary Figs, 4–5, Additional File 2). These analyses demonstrated varying degrees of overlap for the AC clusters between studies, which further underscored the design of our clustering strategy. Still, we would like to emphasize a clear importance of input data quality and differences in initial data resolution which may explain reported differences in gene enrichment profiles. Moreover, the selected markers and lack of consensus regarding nomenclature in the field complicate comparisons and cell type annotations especially in the case of the stromal cell types [46].

Fig. 2

Comparative analyses of ADRCs between donors, gender and BMI. A Baseline characteristics (age and BMI) were similar between the analyzed male and female subjects (2way ANOVA with Šídák's multiple comparisons test). B The 24 ADRC clusters (AC1-7, EC1-4, IC1-9, PC1-3 and RBC) depicted for each individual (male, black circle and female, open circle) as the fraction (%) of all ADRCs. Multiple Mann–Whitney tests revealed no differences between genders. C, D Distribution of C clusters (% of ADRCs) and D subclusters (as % of their respective major cell class) in each individual (male M1–M7 and female, F1–F7). Overall, cluster distribution was similar between donors. E Comparison of subclusters (as % of their respective major cell class) between male and female individuals showing only adipose stem/progenitor cluster AC6 to be significantly different (Multiple Mann–Whitney test, q = 0.003531). F Distribution the major cell classes AC, EC, IC and PC (% of ADRC) as a function patient BMI. No correlation was observed (Pearson). G Nebulosa UMAP showing ADRCs co-expressing PI16 and CD55. H Comparison of the percentage of CD55/PI16 co-expressing cells (% of ADRCs, left and % of ACs, right) showing significantly higher numbers in females (unpaired t-tests, p = 0.0002 in both cases). I Analysis of BMI (kg/m²) and percentage of PI16/CD55 co-expressing ADRCs or ACs revealed no correlation (Pearson, p = 0.9879 and p = 0.8544, respectively). For multiple Mann–Whitney tests the FDR was controlled using the Two-stage step-up method of Benjamini, Krieger and Yekutieli

Donor and gender dependent ADRC heterogeneity is minimal

With the established ADRC sub-clustering design described above, we next investigated donor and gender impact on ADRC composition. Our scRNAseq data were derived from 7 males and 7 females that aside from gender, were comparable regarding age and BMI (Fig. 2A). The impact of gender on the cluster distribution was initially tested using all clusters individually (Fig. 2B-C) and subsequently for each of the four major cell classes: EC, IC, PC, and AC (Supplementary Fig. 6A, Additional File 2). Overall, these analyses, did not reveal significant differences in ADRC composition between genders. Indeed, the variation in subcluster size was minimal between donors (Fig. 2C, D). Also, the distribution of subclusters within each cell class (Fig. 2D, E) appeared similar between genders except for AC6, that constituted a higher fraction of ACs in females versus males (Multiple Mann–Whitney test, q value = 0.0035). The latter seems to be explained by a significantly higher number of CD55/PI16 co-expressing stem cells in females than in males, but this difference was independent of BMI (Fig. 2G–I). This contradicts previous bulk RNAseq studies [10] where the PI16/CD55-stem cell cluster was suggested to positively correlate to BMI (and a negative metabolic status in general). In fact, we observed no BMI-related differences in the distribution of cell classes (Fig. 2F) or subclusters among the donors in general (Supplementary Fig. 6A, Additional File 2). This may possibly be due to the rather discrete variation in BMI (28.7 ± 2.8 kg/m² (mean ± SD); range: 24.7–34.7 kg/m²) or because of a more metabolically healthy status of our patient groups (Fig. 2F) as compared to those previously reported [10]. This was also true when comparing cell class distributions in patients having the significantly (t-test, p = 0.0003) lowest versus highest BMI’s (Supplementary Fig. 6C, Additional File 2). Thus, despite a small difference in CD55/PI16 co-expressing stem cells between genders, our data suggest that heterogeneity in ADRC cluster distribution is minimal between cell donors, at least within the discrete, but clinically relevant BMI and age ranges investigated herein. We can, however, not exclude that the difference in CD55/PI16 co-expressing stem cells could be due to the given disease, which differs between genders. Nonetheless, as a direct cell therapeutic product, the smaller very low-frequency populations often described and emphasized in scRNAseq data such as the CD55/PI16 co-expressing stem cells, may seem less physiologically relevant although the interplay between minor and major populations could potentially be important. Previously, Noreen-Thorsen et al. used RNAseq of unfractionated human adipose tissue and observed no obvious gender-related gene expression differences in SAT [47] underscoring the importance of scRNAseq data to identify even small differences as seen herein. However, this issue awaits to be further addressed in the context of clinical efficacy. For instance, one clinical study did report that the injected hematopoietic stem cell (HSC)/ASC ratio along with the presence of endothelial cells was important for a beneficial clinical outcome and that efficacy apparently was not linked to gender, ethnicity, age or BMI although the yield of SVF cells was higher in females than in males. Higher yields from females were also reported by Cremona et al. who analyzed the composition of uncultured ADRCs/SVF isolated in an ATMP-compliant setting from 300+ patients treated at their hospital. Using an enzyme-based method similar to ours [48], they too observed higher ADRC/SVF yields from females versus males, while the distribution of cell types (surface marker by flow cytometry) and other parameters were similar between genders as also was the case with the patient reported outcome measures. Whether the higher percentage of CD55/PI16 co-expressing stem cells observed here for female donors could impact expansion and the cell composition after ADRC culture, when aiming for allogeneic ADRC transplantations, also remains speculative, but is an interesting point to follow in the future as well. For freshly isolated ADRCs, we recently reported that, at the single-cell level, individual autologous cell therapeutic products were indistinguishable regardless of the clinical response [49]. While this supports small variations in ADRC composition to have less impact on the clinical outcome, it may be limited to the specific method of isolation, disease and patient groups among other factors [1]. Therefore a more general knowledge on how the composition of cell therapeutic products may affect clinical outcomes are highly demanded to develop the stem cell field further.

Conclusion

In conclusion, despite minor individual variations may exist, the ADRC cell therapeutic product our patients have received is very similar and with no or very limited impact of gender and BMI on the distribution of cell types. Indeed, the variation was low between donors in general. This is important for regenerative medicine because it suggests that males and females and donors overall have the same prerequisites for a successful treatment when it comes to autologous cell therapy as well as to be equally suited as cell donors for allogeneic cell therapy using cultured ASCs. However, the finding that PI16 and CD55 co-expression is significantly higher in females versus males may potentially impact therapeutic efficacy and should therefore be considered in mixed gender trials using ADRCs/ASCs.

Methods

Human samples

Lipoaspirates for ADRC isolation (n = 14) were obtained from patients enrolled in two clinical RCT studies using non-cultured ADRCs for cell therapy of 1) male participants suffering from erectile dysfunction (ED) following radical prostatectomy (Clinical trials.gov ID:NCT02240823 [50]) and female participants with breast cancer related lymphedema (BCRL). All patients gave written informed consent before participation.

Regulatory approvals

The pRP-ED-trial followed ATMP guidelines and was approved by The Danish Health and Medicines Authority (2015-005140-33), the Danish National Ethics Committee (2017070103) as well as the Danish Data Protection Agency (no 16/2816). The BCRL-trial was approved by The Regional Committees on Health Research Ethics for Southern Denmark (S-20180117) and registered with the Danish Data Protection Agency (18/51767). Both studies were registered at ClinicalTrials.gov (pRP-ED-trial: NCT02240823 and BCRL-trial: NCT03776721).and performed in accordance with the Declaration of Helsinki and ICH-GCP guidelines. All patients gave written informed consent before participation.

Isolation and preparation of ADRCs for single-cell RNA sequencing.

Adipose derived regenerative cells (ADRCs) were isolated from lipoaspirates utilizing a CE-marked, GMP-compliant device (Cytori Celution IV, Lorem Cytori) as described in details [12, 13, 51]. Excess cell material not being returned to the patient was immediately processed for futher downstream single-cell RNA sequencing. To this end and as recently reported, single cell suspensions were obtained by addition of PBS containing 1% ultrapure BSA and 1 U/µl RNAsin PLUS RNase Inhibitor (Promega, Cat.no. N2615) and passage through a 40 µm Flowmi® Cell strainer (VWR, cat.no. 734-5950) before fixation in methanol and storage at − 80 °C until further processing [8, 15]. In both trials, all cell isolations and preparation of ADRCs for scRNAseq were performed by the same staff using the same procedures.

Single-cell RNA sequencing

Prior to sequencing, fixed cells were thawed, rehydrated and immediately loaded onto a 10 × Genomics Chromium controller (10X Genomics, PN110203). Libraries were prepared following the manufacturer instructions using the 10 × Genomics Single-Cell 3’ v2/v3, Chromium Single Cell B Chip Kit, 48 runs (10X Genomics, 10X Genomics, PN-1000073) and sequenced on an Illumina NovaSeq 6000 System (10X Genomics, 20012850). Cell samples were kept for RNA purification, and RIN-values were measured using a RNA 6000 Nano kit (Agilent, cat no. 5067-1512) (Supplementary Table 1, Addional File 1).

scRNAseq data analysis

Single-cell feature count matrices were created using Cell Ranger software (v. 2.1.1 or 3.1.0), for demultiplexing the raw sequencing output files and aligning to the GRCh38 human reference genome. The output metrics generated by Cell Ranger are listed in Extended Data Table 1. Overall, an average sequence depth of 195 million total reads per sample equalling 32,796 reads per cell was obtained. All subsequent analyses were performed using R Version 4.2.3 and the Seurat R software package (Version 4.3.0) in a standard workflow as previously described [8]. Initially, data was converted to single-cell Seurat objects in which genes were expressed in at least five cells, followed by a filtering step excluding low-quality cells with mitochondrial RNA contents above 15%, or expressing less than 200 genes. Next, doublets were predicted and filtered out using DoubletFinder [52], following the authors' standard workflow with an estimated doublet formation rate of 5%. A total of 111,983 cells passed these quality filtering steps and all datasets merged, normalized, scaled and the top 3000 most variable genes identified for use in the following dimensional reduction step. Possible batch-effects were corrected by Harmony integration using the top 30 principal components, before computing UMAP embedding and clustered using Louvain clustering with the resolution parameter set to 0.5. For patients represented by two samples (resuspension in PBS and SSC following methanol fixation), the average number of cells within each cluster was calculated and used for further statistical analyses.

Differential gene expression analysis was performed for each individual cluster using the FindMarkers-command in Seurat and alldifferentially expressed genes are listed in Extended Data Table 2. Differentially expressed genes for each cluster were analyzed for significantly enriched gene ontology (GO) terms using ShinyGO [17] version 0.77 (and beta-version 0.80) (Extended Data Table 3). GO-terms were also extracted from Enrichr [20] as was cell identity predictions from the CellMarker2024-, Tabula Sapiens [18]- and Tabula Muris databases [19] (Extended Data Table 4). For comparing data sets using Jaccard index, the marker genes were selected based on a log2 fold-change > 0.5, and adjusted p value < 0.05 [10]. Data visualization and analysis of scRNAseq data included the use of scCustomize [53] and Nebulosa [54] packages for R.

Statistics

Statistical analyses were performed using the GraphPad Prism software version 10. Statistical significance was set at p < 0.05. and in each case determined using the appropriate statistical method following normality testing, as denoted in the corresponding figure legends.