Single cell spatial transcriptomics integration deciphers the morphological heterogeneity of atherosclerotic carotid arteries

Domain: neuroscience
Published: 2025-12-18
Editorial outcome: minor_revision

Abstract
The process of arterial atherosclerosis is characterised by accumulation of lipids and fibrous material with accompanying inflammation. As plaques progress, they restrict blood flow and cause rupture, which results in life threatening organ ischemia and dysfunction. Although extensively studied, a clear understanding of plaque heterogeneity and mechanisms that trigger their destabilization remains elusive. Our study reveals the molecular micro-architecture of human carotid artery plaques, using bulk and single-cell RNA sequencing combined with single-cell spatial transcriptomics, for which we present optimized cell segmentation algorithms. We identified distinct plaque morphologies linked to different cell type compositions, impacting early and advanced lesion formation, as well as destabilization. Spatial transcriptomics enabled us to further determine an inflammatory smooth muscle cell subtype, localize regions of neovascularization, and assign hotspots for macrophage activity within distinct cellular neighbourhoods across lesions. For different macrophage substates, we propose gradual and locally contained transdifferentiation of subluminal inflammatory HMOX1+ macrophages into a lipid-handling TREM2+ phenotype within border zones of the fibrous cap and necrotic core. Our findings provide insight into the complex heterogeneity of human atherosclerosis by unravelling location and proximity of different mural and immune cell substates involved in plaque progression and vulnerability.

Introduction
Atherosclerosis is a leading cause of cardiovascular disease and contributes to myocardial infarction and stroke1. Plaques are complex structures that consist of endothelial cells (ECs), vascular smooth muscle cells (VSMCs), fibroblasts, and immune cells, including macrophages and B and T lymphocytes, with each cell type playing a distinct role during lesion formation and progression2.

Atherosclerotic plaque development involves complex cellular and molecular interactions initiated by the focal subendothelial retention of apolipoprotein B-containing lipoproteins3. Focal EC dysfunction increases permeability and inflammation, and blood-born monocytes enter the nascent lesions, differentiate into macrophages, and, along with adaptive immune cells, secrete inflammatory cytokines and promote an inflammatory environment4. Plaque growth is further stimulated by VSMCs that migrate from the media while adapting their phenotype. Both intimal macrophages and VSMCs can become lipid laden (foam cells), presumably due to ingestion of retained lipoproteins in the intima5 and other yet to be determined mechanisms. ECM components, including collagen and elastin, provide structural support, with ECM (dys)balance influencing plaque stability6.

Atherosclerotic plaques from patients undergoing carotid endarterectomy (CEA), the surgical removal of plaque buildup in carotid arteries, consist of three major subregions with distinct cellular and molecular features. The necrotic core (NC), largely acellular, contains dead cell debris and cholesterol crystals. A major mechanism of NC development is secondary macrophage necrosis due to failure of proper clearance, or efferocytosis, of apoptotic macrophages by neighboring lesional macrophages7,8. The fibrous cap (FC), mainly made of VSMCs and collagen, separates the NC from the arterial lumen. The intima-media layer, rich in VSMCs and fibroblasts, represents the pre-atherosclerotic arterial wall. The adventitia is not excised during the surgical procedure and is retained in the patient. All three substructures with their different cellular and molecular components interact in a dynamic and complex manner. However, this interaction continues to be incompletely understood, which is why research on human lesions remains highly relevant9.

Sequencing technologies, including bulk RNA-Seq, single cell (sc) RNA-Seq, and spatial transcriptomics, have helped tremendously to better determine molecular alterations and cellular composition in human diseases10–15. scRNA-Seq biases toward immune cells, therefore offering limited interpretability of fibroblasts, VSMCs, and apoptotic cells, which limits our understanding of plaque progression16,17. Integrating all three methodologies enables high-resolution gene expression analysis, revealing cellular dysfunction and interactions. Spatial transcriptomics preserves spatial context, thus providing key insights into disease mechanisms not captured by other methods.

This study aims to characterise early and late-stage human atherosclerotic carotid plaques, using an integrated multi-transcriptomic approach, combining bulk, single cell, and spatial transcriptomics. Analyzing paired early and advanced plaques from the same patient undergoing CEA, we assessed lesion progression, destabilization, and molecular alterations in plaque development, clarifying dynamic interactions and identifying drivers of plaque progression.

Results
Human carotid plaques from CEA patients were selected based on conclusive morphology, including remnant media and macroscopically assessed and dissected into advanced atherosclerotic (plaque) and adjacent early lesion (control) tissue from the same patient (matched paired diseased-control samples). Plaques were characterized for stability based on the Oxford Plaque Studies18–20. Both pieces of this unique setup (statistically powerful paired analysis possible) underwent bulk- and scRNA-Seq. In addition, formalin-fixed paraffin (FFPE)-embedded samples of neighbouring pieces, maintaining the same paired plaque vs. control configuration, were processed for histomorphological and spatial transcriptomics analyses (Fig. 1).

Transcriptomic landscape of human atherosclerotic carotid plaques
The initial bulk RNA-Seq analysis compared advanced plaque and early lesion control tissues in unpaired (n=232, n=65) and paired (n=42) designs (Fig. 2A, B). Plaque-control comparisons identified 903 dysregulated genes, with 70% upregulated, increasing to 1775 in paired analyses, with 76% upregulated. Key DEGs that were increased in advanced plaque vs. early lesion controls included immune regulators like FLT421 or B-cell markers like CD79A22 and, mitochondrial stress marker MZB123, pro-osteogenic factor STMN224, and extracellular matrix degradation enzymes MMP9 and MMP1325 (Fig. 2B). The advanced plaques were further classified into stable (n=97) and unstable (n=126) lesions based on FC thickness above or below 200µm18. Nine advanced plaque specimens (of n=232) were rendered non-classifiable, and not included into the stability-based comparison.

While a previous study11 identified 107 DEGs between stable and unstable coronary plaques (n=10 in total), our bulk transcriptomes of 223 carotid plaques identified only four DEGs, including CXCL5 and PI3K (Fig. 2C). The lower DEG count in our study reflects lower susceptibility to sample-, technical-, and regional-specific biases. In addition plaque stability, when defined solely by FC thickness as a key predictor of clinical outcomes, seemed to not be primarily driven by transcriptional differences at the whole-tissue level. We could further not detect any DEGs by comparing plaques from symptomatic (n=67) vs. asymptomatic (n=162) patients (FDR threshold < 0.05; minimum FDR in our analysis 0.13).

To resolve cell type-specific gene expression in carotid plaques, we performed scRNA-Seq alongside bulk RNA-Seq on advanced plaque and adjacent early lesion control tissues, integrating data from two studies (Paloschi et al.26 for the first 10 patients and patients 11-17 in this study). The analysis confirmed the expected cellular heterogeneity, identifying ECs, VSMCs, fibroblasts, macrophages, dendritic cells, mast cells, plasma cells, B cells, NK cells, and T cells (Fig. 2D, E)27,28, underscoring the complex immune and stromal landscape of carotid artery disease.

Comparative analysis of scRNA-Seq data identified cell type-specific differentially expressed genes undetected by bulk RNA-Seq in plaque vs. adjacent control tissues (Supplementary Fig. 1A-D), including ALOX5 and HIF1A in IL10high macrophages. Interestingly, stable vs. unstable lesion comparisons revealed an even greater number of DEGs (Supplementary Fig. 2A-K), such as CXCR4 and HSPA5 in proliferating T cells, both markers of T cell infiltration in atherosclerosis29,30. Additionally, CXCL2 was enriched in ECs, while CCL18 was altered in C1Q+ macrophages, both linked to atherosclerosis31,32.

Microfluidics-based capture biases in combination with enzymatic and mechanical digestion is a potential issue with scRNA-Seq, leading to an overrepresentation of immune cells. Indeed, our scRNA-Seq revealed high mean macrophage (24%) and T cell fractions (29%), together exceeding VSMCs (22%) and ECs (15%; Fig. 2F). To correct for this, we applied a deconvolution algorithm33 to estimate cell type proportions in bulk RNA-Seq data using scRNA-Seq-derived labels (Fig. 2G). Deconvolution determined VSMCs as most abundant in plaques (62%) and controls (81%), while macrophages were enriched in plaques (28%) compared to controls (12%; Fig. 2H). Other immune cells, including T and B cells, were rare, comprising less than 0.2% in both plaques and control tissues upon deconvolution. In addition to this deconvolution approach, Xenium-based spatial transcriptomics (designed with an equal amount of probes for all cell types) similar to bulk RNA-Seq confirmed that the observed immune cell dominance (Refs. 34–36) solely occurred in our scRNA-Seq study. Of note, these findings reflect the inclusion of the VSMC-rich media from our CEA specimens (Fig. 1A), a structure that has likely not been analyzed in previous studies that solely focused on the intimal (NC and FC) lesion content.

Spatial component complements transcriptomic landscape
To localise single cell expression signatures within the plaques, we applied single cell spatial transcriptomics, using the Xenium platform on 12 advanced plaques and 4 early lesions/controls. Four patients provided paired samples, meaning that one CEA specimen from the same carotid artery was divided into early and advanced lesions (Supplementary Fig. 3), while eight provided advanced plaques only. Two panels containing probes for 548 genes were designed. Panel 1 identified cell types based on scRNA-Seq data, and panel 2 captured atherosclerotic disease-associated transcripts (Supplementary Data 1). Cell segmentation and transcript assignment were performed using Baysor37, as recommended by Salas et al38. Testing multiple input and parameter combinations, we selected the optimal segmentation parameters based on Baysor’s cell assignment confidence and detected lesion cell types (see “Methods” for details). Manual annotation identified three hierarchical levels: high- and low-level cell types representing primary identities and subtypes (Fig. 3C-D, Supplementary Fig. 4A, B), and low-level substates capturing finer transcriptional differences (Supplementary Fig. 5–8).

Spatial transcriptomics revealed robust cellular profiles across gene panels. Final datasets comprised 2000–17,500 cells per sample for panel 1 (Figs. 3A) and 1000–12,000 cells for panel 2 (Supplementary Fig. 4C), with both panels detecting similar cell types (panel 1: Fig. 3C, D; panel 2: Supplementary Fig. 4A, B), mostly VSMCs and macrophages (Figs. 3B, and Supplementary Fig. 4D). Panel 2 yielded fewer transcript reads per sample, likely reflecting its focus on disease-associated genes rather than cell-defining transcripts (Supplementary Fig. 9A vs. Supplementary Fig. 9B). These findings underscore the value of integrating spatial transcriptomics to refine unbiased cellular and disease-related insights in plaque analysis.

Subregional Plaque Analysis Reveals Spatial Organization and Cellular Composition
To assess plaque morphology and spatial organization, we performed haematoxylin and eosin (HE) staining, identifying the media, intima, NC, and FC (Fig. 4A, B; and Supplementary Fig. 10, Supplementary Fig. 11). This approach provided a spatial framework for downstream analysis. Plaques consist of media, the outermost layer of CEA plaques, the NC as actual plaque region containing lipids and dead cell debris, as well as the FC separating the lumen from the NC. Plaques sometimes display lumen with thrombi and blood cells, while controls consist only of media, intima, and lumen. However, no tunica adventitia that could be the source of fibroblasts, neurons or lymphatic ECs14 is excavated during CEA. A cardiovascular pathologist annotated subregions, which were then converted into boundaries for polygons for area calculations (Fig. 4C). Structural distortions in some samples required approximation of subregion borders.

Clustering of normalized and absolute subregional areas distinguished plaques from controls, which lacked FC and NC regions and showed moderate variability in intima and media areas (Fig. 3E; Supplementary Fig. 4E). Plaques formed three clusters: one with ~30% normalized lumen area and NC below 50% (patients 1, 4, 6, 7), another dominated by NC (~70%) with a narrow lumen (<7%) (patients 5, 9, 10), and an intermediate group with NC ~50% and lumen 10–20% (patients 2, 3, 8, 11, 12). Correlation analysis revealed a strong inverse relationship between normalised NC and lumen areas (Supplementary Fig. 12B), which could imply a constant sum, however, this was not observed (Supplementary Fig. 13). In contrast, absolute area correlations showed no association between NC and FC or intima, suggesting NC expansion reduces the relative proportions of other regions (Supplementary Fig. 12A). A weak negative trend between absolute NC and lumen areas suggests NC growth narrows the lumen without affecting intima and media size.

Associations between plaque morphology, patient metadata, and medical conditions showed distinct trends (Supplementary Fig. 14–17). Patient age inversely correlated with FC absolute area, height correlated positively with total plaque area (Supplementary Fig. 14A, 15C), and hypertensive patients presented with increased normalized media areas (Supplementary Fig. 17A). Plaque stability and symptoms (e.g., stroke or transient ischemic attacks) showed no correlation with subregion morphology.

Clustering based on relative cell distributions across the four main subregions (media, intima, FC, NC) differentiated plaques from controls (Fig. 3F). Controls formed a distinct cluster with a consistent media-to-intima cell ratio (~65% to 35%), except for patient 4 (45% intima). Plaques clustered into two groups: one with most cells around the NC (patients 3, 4, 5, 8, 9, 11, 12) and another with higher media cell fractions and sparsely populated NCs (patients 1, 2, 6, 7, 10). Comparison with normalized region clustering (Fig. 3E) revealed that some plaques with enlarged NCs retained proportional cell fractions (patients 3, 5, 8, 9), while others had disproportionately low NC cell content (patients 1, 7, 10), suggesting advanced necrosis, as confirmed by DAPI-stained spatial transcriptomics (Supplementary Fig. 18).

Control samples consisted mainly of VSMCs with few macrophages in the neointima, while plaques exhibited diverse cell types (Fig. 4A, B; and Supplementary Fig. 10 and Supplementary Fig. 11). In plaques, VSMCs localized to the media and FC, macrophages accumulated around the NC and in the FC, while ECs lined the lumen-but interestingly also appeared intraplaque near the NC where they co-localized with immune cells. In this area, intraplaque haemorrhages of leaky neovessels facilitate immune cell infiltration, promoting subsequent inflammation. Immune cells, apart from macrophages, predominantly accumulated in neovessel-rich regions of advanced plaques, highlighting the link between neovascularization and immune activation (panel 1: Supplementary PDF 1, panel 2: Supplementary PDF 2).

Morphological plaque clustering
We performed hierarchical clustering of plaques using the relative low-level substate cell fractions per sample. This grouped our plaques into four distinct morphological clusters with a unique composition of structural and cellular features in every cluster (panel 1: Fig. 5, panel 2: Supplementary Fig. 19B). For intelligibility we decided to term the two main groups of our plaque clusters “structured” (morphological clusters 1 & 2) and “chaotic” (morphological clusters 3 & 4).

Morphological cluster 1 consists of plaques of patients 7 and 10 (Fig. 6A). Both plaques are organized in around shaped contractile VSMC-rich media, a huge cell free NC and a stable, thick FC (with many “synthetic” fibroblast-like VSMCs), and shoulders rich in HMOX1+ and TREM2+ macrophages. Morphological cluster 2 consists of 4 plaques (patients 1, 2, 6 and 11) that are all VSMC-rich (contractile, inflamed, and fibroblast-like) with panel 1 (Fig. 6B), while macrophages and other immune cells appear underrepresented. These plaques have huge, rather cell free NCs and overall, a similar morphology as morphological cluster 1. Panel 2 shows similar cell type compositions across these morphological clusters. Cluster 1 is characterised by TREM2+ macrophages and VSMC substates (contractile and fibroblast-like), while cluster 2 predominantly consists of different VSMC substates (Supplementary Fig. 19B).

Morphological cluster 3 consists of plaques from patients 4, 8 and 9 (Fig. 6C). They are less structured, showing two NCs on each side of the lumen (potentially one NC in another plane of the plaque; patient 4 and 8) or one large surrounding NC (patient 9), that are cell-rich with many macrophages and fibroblast-like VSMCs. The media is not as present as in morphological clusters 1 and 2 (Fig. 6A, B). With panel 2, these plaques are characterised by macrophages, fibroblast-like and inflamed VSMCs.

Morphological cluster 4 is even more chaotic (Fig. 6D). It is marked by numerous inflamed and fibroblast-like VSMCs in panel 1, as well as C1Q+ macrophages and T-cells in panel 2, infiltrating the cell-rich NC. C1Q+ macrophages are characterised by high SELENOP expression, an antioxidative gene, as well as CD14, which apart from its role in inflammation mediates apoptotic cell clearance39, enhancing phagocytosis and efferocytosis34. Upregulation of CCL8 in these cells is indicative of a “M2-like” polarization40. In plaques from patients 3, 5, and 12, the media is either not present at all or very thin, which could be due to minor differences in surgical removal of plaques during CEA (Supplementary Fig. 10A, and Supplementary Fig. 11B).

To explore regional differences within plaques, we performed hierarchical clustering using relative low-level substate fractions calculated per subregions in each sample. FC-clustering identified four main clusters (Supplementary Fig. 20A), a macrophage-rich (patients 1, 4, 7, 8), an EC-rich (patients 3, 11), one high in contractile VSMCs (patients 6, 10), and one high in fibroblast-like VSMCs (patients 2, 9, 12).

Intima substate fractions clustering (Supplementary Fig. 20B) formed two major clusters: a macrophage-rich with fibroblast-like VSMCs (patients 4, 5, 7, 12) and a VSMC rich (contractile and fibroblast-like; patients 1, 2, 3, 6, 8, 9, 11). The media substate fraction clustering split the plaques into two groups (Supplementary Fig. 20C), one with high fibroblast-like and low contractile VSMC fractions (patients 4, 5, 9, 12), and vice versa (patients 1, 2, 3, 6, 7, 8, 10, 11). Interestingly, samples from patients 3 and 8 had previously organized into the morphological clusters 3 and 4 (Fig. 6C, D), however their media substate fractions appeared more similar to the media of contractile VSMC substate-dominated plaques.

Clustering with NC substate fractions (Supplementary Fig. 20D) identified a group rich in various TREM2hi macrophage substates (patients 4, 8, 9, 10), encompassing plaques of morphological cluster 3, and another group with high fractions of fibroblast-like VSMCs and C1Q macrophages, containing almost all plaques of the other morphological clusters (patients 1, 2, 3, 5, 6, 7, 11, 12). Similar results with slight variations were observed with panel 2 substate clustering (Supplementary Fig. 21A–D). These results show that despite similar morphology and composition, plaques still exhibit significant regional heterogeneity.

Smoking status impacts carotid plaque progression
We noticed a separation pattern with smoking status (“Smoking ever”) in the plaque clustering (Fig. 5). Plaques from current or former smokers grouped into morphological clusters 1 and 2 (except for patient 1), whereas morphological clusters 3 and 4 included non-smoking patients only. This was driven by higher contractile VSMC cell fractions in smokers, seen in both, low-level cell types (Supplementary Fig. 22A–E, Supplementary Fig. 22K–N) and substates (Supplementary Fig. 22F–J, Supplementary Fig. 22O–U). However, we did not observe reproducible differences in VSMC cell fractions across smoking status with either bulk- or scRNA-Seq which reflects rigor in assessing reproducibility of our findings (Supplementary Fig. 23A, B).

Smoking status was associated with the most DEGs across metadata variables in our spatial transcriptomics approach (Supplementary Fig. 24, Supplementary Fig. 25, Supplementary Data 2, Supplementary Data 3). To validate these findings with more robust data, we conducted differential expression analysis with our bulk RNA-Seq dataset across plaque samples (n=99 “smoking” vs. n=106 “non-smoking”; Fig. 7A, and Supplementary Data 4), identifying 44 overlapping DEGs between spatial transcriptomics and bulk RNA-Seq with concordant (48%) (Fig. 7A, and Supplementary Fig. 26, Supplementary Fig. 27, Supplementary Data 5) and discordant (52%) (Supplementary Fig. 19A, Supplementary Fig. 26, Supplementary Fig. 27) log fold change values.

ATM was concordantly upregulated in smokers. This gene has previously been linked to lung cancer development by increasing the risk in non-smokers compared to smokers41. Tissue fibrosis related genes BMPR2 and TGFBR1, associated with smoking in lung fibrosis42,43, were also upregulated in smokers. Inflammatory genes IL6 and IFNG were downregulated in plaques of smokers, while ADAMTS7, downregulated in smokers, was reported to be lower in current smokers vs. non and ex-smokers in human CEA specimens44. Fig. 7B shows the spatial distribution of ATM and TGFBR1 as representative genes in two plaques (patients 7 and 10).

Although plaques rendered as stable were enriched in morphological clusters 1 and 4, no obvious cell substate fractions explained this clustering pattern. In addition, no other patient metadata variables showed discernible patterns of distribution across any of the morphological clusters (Fig. 5, and Supplementary Fig. 19B).

Cellular neighbourhoods determine the most prominent plaque features
To fully utilize the potential of our combined single cell and spatial transcriptomics data, we next investigated the composition of our samples on the cellular scale. Our goal was to identify similar local cell neighbourhoods across samples that could reveal local processes impacting plaque development and destabilization.

After defining the local neighbourhood of a cell (counting nearby low-level substates within a 100μm radius), we detected 21 neighbourhood clusters in panel 1 and 26 in panel 2 (see “Methods” section for details). We averaged the local neighbourhoods of cells within each neighbourhood cluster, creating local neighbourhood signatures (Supplementary Fig. 28A, B). We next focused our analysis on neighbourhood clusters 0, 2, 3, 4, 8 and 11 of panel 1 due to their interesting local signatures (Fig. 7C).

Neighbourhood clusters 2 and 4 were identified as “macrophage neighbourhoods”, with different macrophage substates residing in close proximity to each other. Cluster 2 mainly consists of substate C1Q+, a known enhancer of efferocytosis45, and the substate of HMOX1+ involved in heme degradation34. Cluster 4 is dominated by TREM2hi macrophages35,46, lacking overexpression of inflammatory genes34. Substate Mac_TREM2hi_CCL18hi, a slightly transitioned state of classical TREM2hi macrophages, localises to both clusters, but is more prevalent in cluster 2.

In plaques of patients 4, 8 (Figs. 8A), 7, and 10 (Supplementary Fig. 29A), most of the four macrophage substates reside in the crucial, rupture-prone shoulder region within FCs, displaying an interesting common spatial pattern from lumen towards NC. C1Q+ and HMOX1+ macrophages are located subluminally, followed by Mac_TREM2hi_CCL18 and finally TREM2hi macrophages (Fig. 8C, and Supplementary Fig. 29B). Spatial correlation analysis underpinned this and showed significant colocalization of C1Q+ and HMOX1+ macrophages near the lumen. Neither of them showed spatial association with TREM2hi cells in any other lesion area, a pattern that is consistent in all four plaques (Supplementary Fig. 30, and Supplementary Fig. 31).

Interestingly, the colocalization pattern of C1Q+, HMOX1+ and TREM2hi appears in other samples, but not necessarily in the FC (Supplementary Fig. 32–34).

The distinct spatial distributions of macrophage substates suggest differential roles in plaque development and progression. Subluminal HMOX1+ macrophages, beyond their role in heme degradation, highly express the antioxidant SELENOP, as well as the known lipid handling mediator APOE (Fig. 8B, and Supplementary Fig. 35E). TREM2hi cells, located deep in the shoulder regions of FCs, express SPP1 at very high levels, as well as the ECM regulators FN1, VIM, MMP9, and MMP12 (Supplementary Fig. 35D, E), suggesting their relevance in fibrosis47 in addition to lipid handling. These TREM2hi macrophages seem likely involved in plaque shoulder remodelling within the FC48. A decreasing gradient of expression in CCL18, CTSD, and LYZ was observed from HMOX1+ to TREM2hi macrophages, while SPP1, MMP9, and MMP12 show the opposite trend (Fig. 8B, and Supplementary Fig. 35A). Based on our spatial and DEG results, we hypothesize a gradual transdifferentiation of HMOX1+ into TREM2hi with Mac_TREM2hi_CCL18 as an intermediate, transient state. We assessed this with velocity analysis using TFvelo49, which showed increasing mean pseudotime along subtypes Mac_C1Q → Mac_HMOX1 → Mac_TREM2hi_CCL18hi (panel 1) → Mac_TREM2hi (Supplementary Figs. 36, 37), accompanied by a lumen-to-NC gradient in patients 4, 7, 8, and 10 (Supplementary PDF 3, 4). Partition-based graph abstraction (PAGA) supported the same progression (Supplementary Figs. 38C, 39C), though the link between Mac_TREM2hi_CCL18hi and Mac_TREM2hi was weak in panel 1. Despite not being fully concordant across all subtypes, complementary analyses of our scRNA-Seq data with two orthogonal methods (TFvelo, scVelo50) yielded consistent trajectories from HMOX1+ to TREM2hi cells and lower pseudotime values for HMOX1+ relative to TREM2hi (Supplementary Figs. 40; Supplementary Fig. 41), suggesting a possible transdifferentiation route that requires future experimental validation.

Neighbourhood cluster 8 appears specific to the control tissue of patient 4 (Fig. 8A), an early lesion with dawning infiltration of macrophages in the intima layer. We observed a subluminal accumulation of HMOX1+ macrophages, with neighbouring VSMCs and EC_2 cells, which display an ECM modulation phenotype by highly expressing ITLN1, CDH11 and OMD51,52. Interestingly, the subluminal C1Q+ and HMOX1+ cells are not spatially separated from TREM2hi cells (Supplementary Fig. 32), suggesting that this separation occurs later during plaque advancement.

Neighbourhood cluster 11, termed “EC neovessel” neighbourhood, was uncovered only in two plaques (patients 4, 8). It consists of HMOX1+ macrophages and EC_1 substates, with the THBS1+/SERPINE+ substate dominating the cluster (Fig. 9A). This appears to be an area of neovascularization, as THBS1 and SERPINE1 are closely linked to angiogenesis53 and EC dysfunction54. Interestingly, this neighbourhood is located in shoulder regions of plaques, rather than luminally, suggesting a potential entry point for immune cells into the lesion core through leaky neovessels formed by pro-angiogenic EC_1 cells. These cells upregulate ACKR1, recently suggested to play a role in human plaque neovascularization and immune cell trafficking55.

Neighbourhood clusters 0 and 3 were termed “VSMC neighbourhoods”, as they mainly consist of contractile VSMCs along with “synthetic” fibroblast-like and inflamed VSMC substates (Fig. 9B), detected in the media of plaques (patient 7, 10: Fig. 9B; patient 1, 3, 11: Supplementary Fig. 42) and control (patient 3, Fig. 9B). Only very few “classical” fibroblasts could be identified overall. As mentioned above, the surgically removed CEA specimens lack the adventitial layer where these cells typically reside in human arteries56.

Cluster 3 presents a high abundance of “classical” contractile VSMCs (substates MYH11hi and COL3A1low) and is situated on the outer ring of the media. This cluster seems to represent remnants of the healthy arterial wall with high contractile VSMC density. Cluster 0 is positioned more inwards facing the NC and FC. It contains similar transitioning VSMC substates as cluster 3, but lower contractile VSMC counts. This cluster represents areas of media with less densely packed structures, based on the lower average substate counts in its signature (Fig. 7C). This suggests that the media not only consists of contractile VSMCs but offers a gradual shift from contractile to modulated cell states towards the lesion core.

A recent publication identified plasma and B-cell rich areas, interspersed with VSMCs and macrophages as plaque tertiary lymphoid organs (PTLOs)12. In our samples, spatial analysis revealed significant plasma cell enrichment in the plaques of patients 6 and 7 (Supplementary Fig. 43A, B; Supplementary Fig. 44A, B). In addition, the plasma cell aggregation domain in the NC of patient 6 was identified as cluster 13 (Supplementary Fig. 45), enriched in VSMCs, macrophages and plasma cells (Supplementary Fig. 28A). Although suggestive of PTLOs, only two (MZB1, XBP1) out of four PTLO marker genes12 in our panels showed high expression in plasma cells (Supplementary Fig. 43C, D; Supplementary Fig. 44C, D), highlighting the need for further clarification with more specialised gene panels.

By performing neighbourhood clustering for panel 2, we detected similar neighbourhoods with similar features in the same plaques, with neighbourhood 0 and 1 being “VSMC neighbourhoods”, neighbourhoods 2, 5 and 8 being “macrophage neighbourhoods”, neighbourhood 10 being the “EC_Neovessel neighbourhood”, and neighbourhood 13 the “plasma-cell neighbourhood” (Supplementary Fig. 28B, C).

Methods
Humansamplecollection
HumansampleswereobtainedfrompatientsdiagnosedwithcarotidstenosisandsourcedfromtheMunichVascularBiobank.ThecollectionprocessadheredtoTheethicalprinciplesoftheDeclarationofHelsinkiandwasapprovedbythelocalethicscommitteeattheTechnicalUniversityofMunich(approvalnumbers:2799/10and2023-297).Informedwrittenconsentwasobtainedfromallparticipants,andclinicaldatawereretrievedfromelectronicpatientrecords.
Plaquesampleswerecollectedfrompatientsundergoingcarotidendarterectomy(CEA),thesurgicalremovalofplaquebuildupincarotidarteries.Immediatelyafterremovalfromtheintraoperativesitus,plaquesamplesweretransferredtopre-chilledPBS(Merck,Germany)fortransportation.
Thecarotidplaquesarethencutintoseveralpieces.Specifically,forthisstudy,weselectedthemostadvanced(i.e.,mostunstable)partoftheplaqueandwheneverpossiblealsotheleastdiseasedpartoftheplaqueforanalysis(SupplementaryFig.3).Thisuniquesetupallowsforstatisticallypowerfulpairedanalysesoftwopiecesofthesametissuefromthesameindividualpatient.

Histologicalanalysis
Tissuesampleswerefixedin4%paraformaldehyde(PFA)for24handsubsequentlydecalcifiedforfivedaysusingEntkalkersoftSOLVAGREEN®(CarlROTH,Karlsruhe,Germany).Followingdecalcification,thespecimenswereembeddedinparaffinandstoredasformalin-fixedparaffin-embedded(FFPE)blocks.Paraffin-embeddedsections(2μmthick)weremountedonglassslides(MenzelSuperFrost,76×26×1mm,FisherScientific,Schwerte,Germany).
Hematoxylin-eosin(HE)staining(ethanolic eosin Y solution, Mayer’s acidic hemalum solution, Waldeck, Münster, Germany)wasperformedaccordingtothemanufacturer’sprotocols.SlideweremountedusingPertex(Histolabproducts,Askim,Sweden)asthemountingmedium,withglasscoverslips(24×50mm,Engelbrecht,Edermünde,Germany).AllhistologicalsectionswerescannedusingtheAperioAT2(Leica,Wetzlar,Germany),andslideanalysiswasconductedusingtheopen-sourcesoftware72.ThehistologicalclassificationofcarotidplaqueswasperformedfollowingtheguidelinesestablishedbytheAmericanHeartAssociation(AHA73–76)andtheRedgrave&Rothwellcriteria(OxfordPlaqueStudy)18forplaquestability.

TotalRNAsequencingofatheroscleroticplaques
TissuesamplesforBulkRNASequencingweretransferredtochilledRNALater(ThermoScientific,Germany),andstoredat-80°Cawaitingfurtherprocessing.Total(bulk)RNAwasextractedfromhumancarotidarteryspecimensprovidedbytheMunichVascularBiobankandgeneratedlibrariesthatweresequenced(RNA-seq)viatheIlluminaNovaSeqplatform(NovaSeq6000,Illumina,SanDiego,USA).RNAwaspurifiedusingpoly-T-oligo-attachedmagneticbeadsandtheTruSeqstrandedtotalRNAkitandTruSeqstrandedmRNAkit(Illumina,SanDiego,USA)wereusedtopreparetheshort-readsequencinglibraries.RNAandlibraryintegritywasassessedviaQubit(ThermoFisher,USA),aswellasreal-timePCRforquantificationandtapestation(Agilent,USA)forsizedistributiondetection.
StrandedRNA-SeqwasconductedusingtheIlluminaNovaSeq6000platform,with101basepair(bp)paired-endreads.Sequencingdepthrangedfrom20millionto140millionreadsperample,dependingonRNAquality.TheRNA-Seqreadsweredemultiplexedandalignedtothehg19genomeassembly(UCSCGenomeBrowserbuild)usingSTAR(v2.7.0a).ReadquantificationwasperformedwithHTSeq-count(v0.6.0)todeterminethenumberofreadsmappingtoannotatedgenes.GeneswithanFPKM(FragmentsPerKilobaseoftranscriptperMillionmappedreads)belowthe95thpercentileof1werefilteredoutasinsufficientlyexpressed.
Afterfiltering,14,572genesweredetectedacrossatotalof297samplesderivedfrom252patients.Thedatasetcomprisestranscriptomesfrom232advancedplaquesand65earlylesioncontrols.Amongtheplaques,97wereclassifiedasstableand126asunstable,with9remainingunclassifiedbasedonourhistologicalanalysis.For42patients,pairedsamplesofplaqueandcontrollesionswereavailable,allowingforpaireddifferentialexpressionanalysis.
Batch effectscorrectionwasperformedusingcomBat-seq,adjustingfortheinfluenceofthe“batch”variable.DESeq2(v1.34.0)wasemployedfornormalizationanddifferentialgeneexpression(DEG)analysis.DESeq2normalizedreadcountsbyestimatingsizefactorstoadjustfordifferencesinsequencingdepthacrosssamples.Adesignformulawasconstructedtomodelthecountdatausingnegativebinomialgeneralizedlinearmodels(GLMs).Theformulaincorporatedage,gender,anddiseasestatusascovariates.Forpairedanalysis,individualIDsanddiseasestatuswereincludedtoaccountforwithin-samplecorrelations.TheWaldtestwasusedtoevaluatesignificantdifferencesingeneexpressionbetweenpatientsandcontrols,whilecontrollingfortheeffectsofcovariates.DEGs wereidentifiedusingtheresultsfunctioninDESeq2,withBenjamini-Hochberg-adjustedp-values(falsediscoveryrate,FDR)toaccountformultipletesting.GeneswithanFDR<0.05and|log2foldchange>1|wereconsideredstatisticallysignificant.Forthesmoker-non-smokercomparison,differentialexpressionwasperformedwithpyDESeq2(DESeq2framework)usinganegativebinomialGLMwithdesign~Smoking_ever.Theanalysiswasrestrictedtoplaquesamplesonly.Smokersandformersmokerswerecollapsedintoasingle‘Smokers’groupandcontrastedagainst‘Nosmokers’(reference).Noadditionalcovariates(age,gender)wereincluded.Waldstatisticswereusedandtwo-sidedp-valueswereadjustedusingtheBenjamini-Hochbergprocedure.StatisticalsignificancewasdefinedasFDR<0.05and|log2foldchange|≥0.5.ThefullresultsforallassessedgenesareprovidedinSupplementaryData4,whileSupplementaryData5containsthesubsetofgenesalsoincludedinourXeniumpanels.

SinglecellRNAsequencingofhumancarotidplaques
ForsinglecellRNAsequencingofhumancarotidarteries,sampleswereharvestedduringCEAinourDepartmentforVascularandEndovascularSurgery.Intotal17plaquesweremincedanddigestedusingtheMultiTissueDissociationKit2(MiltenyiBiotech,130-110-203)andthe37C_Mulit_Gprogram(GentleMACSDissociator;MiltenyiBiotech,130-093-235)accordingtothemanufacturer’sinstructions.ThecellsuspensionwascleanedupandfinallyresuspendedinPBS+0,04%BSA.Cellswereloadedintoa10xGenomicsmicrofluidicsChipGandencapsulatedwithbarcodedoligo-dT-containinggelbeads.Gel-Beads-in-emulsion(GEM)cleanup,cDNAamplificationand3´GeneExpressionLibraryConstructionwasperformedaccordingtothemanufacturer´sinstructions(CG000204RevD).Librariesfromindividualpatientsamples1-10(batch1)weremultiplexedintoonelanebeforesequencingonanIlluminaNovaSeq6000instrument.Onewholeflow-cellwasoccupiedandusedforsequencing.ThisdatahasbeendepositedinNCBI’sGeneExpressionOmnibus77andisaccessiblethroughGEOSeriesaccessionnumberGSE247238.Librariesfromindividualpatientsamples11-17(batch2,GSE294291)wereloadedonanIlluminaNovaSeqwith2×150paired-endkitsandsequencedwithNovogene(Martinsried,Germany).AnadditionalscRNAseqdatasetofhumanplaqueswasdownloadedfromGEOrepository(GSE159677)51.
Alignment.ForthealignmentacustomreferencegenomewasbuiltusingtheGencodeGRCh38-p14-V44primaryassemblyFASTAfileandmatchingGencodeprimaryassemblyannotation.gtffile.Forthisweusedthemkreffunctionof10x’scellrangerversion7.1.0,includingthefollowingbiotypes:“protein_coding”,“lncRNA”“antisense”,“IG_LV_gene”,“IG_V_gene”,“IG_V_pseudogene”,“IG_D_gene”,“IG_J_gene”,“IG_J_pseudogene”,“IG_C_gene”,“IG_C_pseudogene”,“TR_V_gene”,,“TR_V_pseudogene”,“TR_D_gene”,“TR_J_gene”,“TR_J_pseudogene”,“TR_C_gene”.Alignmentofthereadswasperformedwithcellranger’scountfunctionwithdefaultparameters,whichincludedintronicreads.
Preprocessing steps.BackgroundRNAcorrectionwasperformedwithcellbender(version0.3.0)78,dropletswithacellprobabilitylargerthan0.5predictedbycellbenderwereretainedforanalysis.Single-cellanalysiswasperformedwithscanpy,version1.9.879.Cellswithhigherythrocytemarkergeneexpression(HBB(HemoglobinSubunitBeta),HBA1(HemoglobinSubunitAlpha1),HBA2(HemoglobinSubunitAlpha2))wereremovedbeforerunningdoubletdetectionwithscDblFinder(version1.12.0)80.Qualitycontrolentailedretainingcellswithfollowingparameters:genecountpercellbetween200and55000,numberofgenespercellbetween200and8000,scDblFinderdoubletscorelowerorequalthan0.5,mitochondrialreadratiolowerorequalthan0.2.Genesexpressedinatleast5cellswerekeptinthefinaldataset,log-normalisationwascarriedoutbyapplyingthe“acosh”transformationtothecellbendercorrectedcountvalues.
Batchcorrectionandcelltypeannotation.WeplannedtoannotatecelltypesmanuallyinasimilarfashionaspreviouslydescribedinthecaseoftheXeniumdata:detectingclustersafterbatchcorrectionandannotatingthembasedontheirdifferentiallyexpressedgenes.Forthisbothlinear(Harmony,version1.1.0)81,andnonlinearautoencoderbased(scVI82withscvi-tools1.1.2;scANVIwithscarchesversion0.6.183)batchcorrectionmethodswereappliedonourdataset,butnoneofthemwereabletocompletelyremoveallofthem,eventuallyconfoundingthefurtherstepsofclusteringandDGEanalysis.Thereforeweadoptedadifferentapproach:wedecidedtointegrateourdatawithalargerplaquescRNA-seqdatasetofasimilarexperimentalsetuptooursandtransferitscelltypelabels,hypothesisingthatcelltypeswerealsosimilaracrossthetwodatasets.WeselectedtheplaquescRNA-seqdatasetpublishedbyAlsaighetal51,whichhadasimilarexperimentalsetup(matchedatheroscleroticcore-proximaladjacentsamplesfromsamepatients)toours(matcheddiseased-controlsamplesfromsamepatients).However,ascelltypeannotationswerenotpublishedalongsidethedataset,wehadtoannotateitmanually.ThedatasetisaccessibleundertheaccessionnameGSE159677intheGeneExpressionOmnibus(GEO)repository(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE159677).
Celltypeannotation-GSE159677dataset.WedownloadedthefilteredcountmatrixfromtheGEOrepository,ranscDblFinderfordoubletdetectionandkeptcellswithfollowingparameters:genecountpercellbetween400and55000,numberofgenespercellbetween350and5000,scDblFinderdoubletscorelowerorequalthan0.5,mitochondrialreadratiolowerorequalthan0.1.Log-normalisationagainoccurredby“acosh”transformation.BatchcorrectionwascarriedoutwithscVI82,usingthetop4000highlyvariablegenes,correctingfor“patient”asbatcheffect.
ManualannotationoftheGSE159677datasetwasperformedinasimilarfashio totheXeniumdataset:coarseLeiden-clusteringanddifferentialgeneexpression(DGE)betweenthecoarseclusterswithWilcoxon-methodprovideduswithmarkergeneoftheclusters,whichwerecomparedtomarkergene collectedfromliteraturetodeterminehigh-levelcelltypes.Todetectcellsubtypes,weperformedadditionalLeidenclusteringanddifferentialgeneexpression(DGE)analysiswithinthehigh-levelcelltypes,usingapseudobulkapproachwithDESeq2(version1.38.3)84,controllingfor“patient”asafactorandcreatingpseudobulkcountsforeachclusterbyaggregatingcellscomingfromthesamesample.DESeq2wasappliedwithdefaultparameters,returningtwo-sidedWaldtestp-values,whichwerethencorrectedformultipletestingusingtheBenjamini-Hochbergmethodandreportedasfalsediscoveryrates(FDR).Clusterswereannotatedbasedontheirmarkergeneandcellsubtypemarkergene fromliterature.
Celltypeannotation-integrationandcelltypelabeltransferfromGSE159677withsubsequentrefinement.ForthiswechoosescPoli(scarchesversion0.6.1)85,anautoencoderbasedmodeldesignedtoregressoutbatcheffectsacrosydatasetswhileretainingbiologicalinformationofcelltypespresent.FirstareferencemodelwascreatedbytrainingscPoliontheGSE159677data,usingourcellsubtypelevelannotationsascelltypelabelsandcontrollingforfactors:“patient”,”-batch”,”condition”,”sequencingmethod”.Then,wefurthertrainedthereferencemodelonourdatacontrollingforthesamefactors,resultinginanintegratedlatentrepresentationofthetwodatasets.Ofnote:basedonthemethodssectionofAlsaighet.al’spublicationwedeemedtheirsampleswithconditionlabel“necroticcore”equaltoour“plaque”,butwewerenotsureif“proximaladjacent”meantaproximalhealthyorlessaffecteddiseasedvascularsegment,thereforeweused3conditionlabelsfortheintegration:“healthy”,”plaque”,and“adjacent”.Celltype labels weretransferredtoourdatasetbythefullytrainedscPolimodelalongwithapredictionuncertaintyscoreforeachcellrangingfrom0to1,showinggroupsofcelltypestransferredwithbothlowandhighuncertaintylevels.InspectingtheintegratedlatentrepresentationrevealedcellswithhighpredictionuncertaintyscorestodisplaylowerintegrationlevelswiththeGSE159677dataset.Thislackofintegrationforcertaincellgroupsmaystemfromtechnicalfactors,suchasunknownbatcheffectsorscPoli’smodelarchitecture,orfromlimitedbiologicalvariabilityinthecellstatesoftheGSE159677dataset,giventhat,despitethehighcellcount,thedataisderivedfromonlythreepatients.Weperformedtheintegrationwithtwoothermodelsaswell(scVIandscANVI),buttheintegratedlatentrepresentationsshowedworseintegrationcomparedtoscPoli.However,thesubsetofthescPolilatentrepresentationcontainingonlyourcellsdemonstratedamoreeffectivebatcheffectremovalthantheoriginalsingle-datasetattemptswithHarmony,scVIandscANVI,thereforewederivedthefinalcelltypesbyperformingthesamemanualannotationpipelineonourdataasdescribedfortheGSE159677data,usingthescPolirepresentationasbatchcorrectedlatentrepresentation.DESeq2wasrunwithdefaultsettingsasbefore(two-sidedWaldtestwithBenjamini-Hochbergcorrection).DGEresultsacrosshigh-levelcelltypesandlow-levelcelltypesarecontainedinSupplementaryData6andSupplementaryData7respectively.
Differentialgeneexpressionanalysis.Fordifferentialgeneexpression(DGE)analysisweagainutilisedthepseudobulkapproach:wecreatedpseudobulksamplesbyaggregatingcellsofthesamesamplewithinlow-levelcelltypesandsubsequentlycarriedoutDGEanalysisformultipleconditionvariablesusingDESeq2.Weincludedcellsfromallsampleconditions(plaqueandcontrol)intheDGEanalysis.Theeffectsof“patient”,“batch”and“condition”werecontrolledforinDESeq2,exceptforcomparingacrossconditions,whereonly“patient”and“batch”wereincludedascovariates.DESeq2wasagainrunwithdefaultparameters,producingtwo-sidedWaldtestp-valuesthatweresubsequentlycorrectedformultipletestingusingtheBenjamini-Hochbergmethod.DGEresultscanbefoundinSupplementaryData8.
RNAvelocityanalysisofmacrophages.Asourspatialtranscriptomicdatasetrevealedintriguingspatialandexpressionalpatternsinmacrophagesubtypes,wesetouttovalidatethisfindinginourscRNA-seqdata.Todothis,weinvestigatedthetranscriptionaltrajectoriesofmacrophagesubtypesbyperformingRNAvelocityanalysisusingtherecentlypublishedTFvelo49method,whichutilisestranscriptionalregulationtocalculateRNAvelocitytrajectories.Wefirstsubsettedthedatatoonlycontainmacrophages(2,805cells)andranthepreprocessingstepwithparametersrecommendedbytheauthors,usingthescPoli-UMAPrepresentationasalow-dimensionalrepresentationandacoshnormalisedcounts.Forthecalculationofvelocitystreamplots,weonlyincludedgeneswithalargerlatenttime-pseudotimeSpearmancorrelationof0.4.Pseudotimeandvelocitystreamplots(SupplementaryFig.40A,B)werecreatedusingthescPoli-UMAPrepresentation.Inadditiontothecellular-leveltrajectoryinformationprovidedbythestreamplots,wealsoaimedtoassesstrajectoriesatthecelltypelevel.Tothisend,weappliedthePAGA(partition-basedgraphabstraction)algorithm86usingtheinferredTFvelopseudotimeaspriorsandthemacrophagelow-levelcelltypesasnodegroups.Theresultinggraph,displayingvelocity-directededgesbetweencelltypes,isshowninSupplementaryFig.40C.WealsoplottedtheTFvelopseudotimedistributionofthethreemacrophagesubtypes(Mac_C1Q,Mac_HMOX1,andMac_TREM2hi;SupplementaryFig.40D),whichshowedorderedspatialdistributionpatternstheXeniumdata,tobettervisualizetheirrelativeorderingalongtheTFvelopseudotimeaxis.
To complement TFvelo with an orthogonal RNA velocity approach,wealsoappliedasplicingcount-basedmethodtoinfercelltrajectories.SplicedandunsplicedexpressioncountswereobtainedusingVelocyto(version0.17)50onthe.bamfilesandfilteredbarcodesgeneratedduringalignmentwithCellRanger’scountfunction.
Forthesplicingcount-basedanalysis,weusedscVelo(version0.3.2)87.Wefirstcreatedasubsetofourdatacontainingonlymacrophagecells(1,837cells)andcalculatedfirst-andsecond-ordermomentswithscvelo.pp.moments()usingthescPoli representation,withthenumberofneighbourssetto30.Next,weestimatedvelocityvectorswithscVelo’sdynamicalmodel.InferredpseudotimeandvelocitystreamplotsweregeneratedusingthescPoli-UMAPrepresentation(SupplementaryFig.41A,B).APAGAgraphbasedonpseudotimeinferredbyscVeloisshowninSupplementaryFig.41C,togetherwiththescVelopseudotimedistributionacrossthethreemacrophagesubtypes(SupplementaryFig.41D).

Spatialtranscriptomicsofhumancarotidplaques
Xeniumsamplepreparationandrun.Tissuesectioningandsamplepreprocessingwereperformedaccordingtothemanufacturer’sprotocol.Briefly,5μmFFPEsectionsweremountedonXeniumslidesanddriedtooptimizeadherence.Subsequently,thesampleswentthrougha3-dayprotocolincludingmultiplestepssuchasdeparaffinization,decrosslinking,probehybridization,ligation,andamplification.Thepre-processedsampleswereanalyzedusingtheXeniumAnalyzerwiththeinstrumentsoftwareversion1.5.1.2or1.7.6.0andanalysisversion1.5.0.3or1.7.1.0.Subsequently,H&Estainingofthetissuesectionswasperformedaccordingtothemanufacturer’sprotocolforpost-Xeniumstains.
Panel Information.Targetedspatialtranscriptomics(Xenium,10XGenomics)wasrunwithtwofullycustom-madepanels:“9MHNNK-FullycustomhArterypanelforTUM”with314individualgenes(panel1)and“PXD8AD-FullycustomhArterypanelforTUM”with345individualgenes(panel2).Geneswerechosenbasedoncellularanddiseaseassociatedpathways(apoptosis,DNAdamage,efferocytosis,CADassociatedgenes,cellcycle,inflammasome,macrophagebiology,senescence,unfoldedproteinresponse)aswellasDEGs,markergeneandmarkerproteinsfromplaque27,58,88andGWAS89datasets,includingourownscSeqdataset(GSE247238,usedin26,55,90)(SupplementaryData1).Thepanelsconsistofacorelistpresentinbothpanels(comprising111centralgenesforcelltypeannotationsandkeyprocesses)aswellasindividualgenesforthepanels,resultinginacell-defining(panel1)anddiseaseassociated(panel2)panel.Intotal548uniquegeneshavebeensequencedwithbothpanels.
Sequencingwithpanel1andpanel2geneswasperformedonconsecutiveFFPEsectionsforeachsample,therebyintroducinga5µmshiftinthemorphologyofthemountedXeniumslidesacrossthepanels.
CellsegmentationandtranscriptassignmentwithBaysor.XeniumOnboardAnalyzer(XOA)outputsthe3-dimensionalcoordinatesofthetranscriptsandbordercoordinatesofnucleibyperformingnuclearsegmentationusingthehigh-definitionDAPI-stainedtissueslides.Furthermore,astandardcellsegmentationwasalsoprovided,thisinvolvedtheexpansionofthenuclearsegmentationbordersbyaconstantnumberandassignedwithinbordertranscriptstotherespectivecells.AbenchmarkingpaperbySalasetal38.showedthatthisapproachcouldleadtomisassignmentoftranscriptsfromneighbouringcells,andrecommendedbestpracticesforcellsegmentationandtranscriptassignmentofXeniumdata.Toensurethebestpossiblecellsegmentation,wehavefollowedtheirbestpracticespipelineandperformedourowncellsegmentationandtranscriptassignmentwithBaysor,version0.6.237.Baysorusesspatialtranscriptomicdatatoinfercellborders,optionallytakingimagedataofnuclearorcellstainingintoaccountinadditiontothetranscriptcoordinates.
WeconsideredfourdatacombinationsasinputtoBaysor:spatialtranscriptomicdataalone,spatialtranscriptomicdatapairedwith10xXOAnuclearsegmentationmasks,andspatialtranscriptomicdatawithsegmentationmaskscreatedfromtheDAPI(4’,6-Diamidino-2-phenylindol)stainedslideswitheitherThe‘’Cytoplasm’orThe‘Nucleus’modelofCellpose(version2.2.2)91.ThemostsensitiveparameterofBaysoris“scale”,whichcorrespondstotheapproximateradiusofthecellsinthesample,alongwith“scale-std”,whichsetsthestandarddeviationofcellradiusacrosscells.Asplaquesnormallycontaincelltypeswithdifferentshapesandsizes(i.espindle-likevascularsmoothmusclecells92,orroundmacrophages93,94),wehavesetthe“scale-std”parameterto50%insteadofthedefault25%,andranBaysorwith3increasing“scale”parameters(5,10,15μm)foreachofthefourinputcombinations,therebyperformingcellsegmentationswith12differentdatainput-scalecombinationsforeachsample.Thiswaywewereabletochoosethebestperforminginputdata-scaleparametercombinationforeachsamplelater,asthereweresignificanttechnicalandbiologicalvariationsbetweensamples.Ofnote,Baysorprovidesthepossibilitytoestimatethe“scale”parameterforeachsample,howeveraswesawlargevariationintheestimated“scale”parameteracrossthesamples,wedecidedtousethe3fixvalues(5,10,15μm)approach.Asper10xGenomicsrecommendation(https://www.10xgenomics.com/support/software/xenium-onboard-analysis/latest/analysis/xoa-output-archive-data),wefilteredoutlow-qualitytranscriptreads(Q-Score<20)beforerunningcellsegmentation.
Baysoroutputcellboundaryandcellcentrecoordinatesoftheidentifiedcells,aswellasthetranscriptsassignedtothesecellsalongwithconfidencemetricsforthecellassignments.Wekepttranscriptswithassignmentconfidence>0.5,ametricthatreflectedBaysor’sconfidenceinassigningthetranscripttothecorrectcell.AcountmatrixforeachsamplewasconstructedbysummarisingthetranscriptsforeachcellinthesampleandsavedasanAnnDataobject,withtherowscorrespondingtocellsandcolumnstotranscripts.Forthesubsequentprocessingstepsweusedscanpy79.
Selection of best cellsegmentation per sample for each “scale” parameter.WereliedonBaysor’scellassignmentconfidencevaluestoselectthebestcellsegmentationforeachsample.Withinasample,wecalculatedtheaveragetranscriptassignmentconfidenceofthecellsforeachdatainputcombinationandplotteditsdistributionviaakerneldensityestimation(KDE)plot.Themajorityofcellshadanaveragetranscriptassignmentconfidenceintherangeof0.75-1,thereforewecalculatedtheareaunderthecurveoftheKDEcurvesinthatrange(AUC0.75)andselectedthedatainputcombinationwiththehighestAUC0.75valueasthebestcellsegmentationforthatsample.
However,asthe“scale”parameterhadalargeeffectonthenumberandsizeofthesegmentedcells,wedecidedtoperformtheselectionforinputcombinationswithineach“scale”parameter,resultinginonebestcellsegmentationper“scale”parameterforeachsample.Thenextsteps(qualitycontrolandnormalisation)wereperformedontheAnnDataobjectsderivedfromthebestcellsegmentations,i.e.,3AnnDataobjectspersample,oneforeach“scale”parameter.
Quality control, normalisation and batch correction.Wefilteredcellswith:lessthan10non-zero-countgenespercell;cellswithloweraverageassignmentconfidencethan0.75;cellswithlower“max_cluster_frac”valuethan0.9,aslow“max_cluster_frac”valuesoftenbelongtodoublets;cellswithlower“lifespan”valuethan50(lower“lifespan”valuescorrespondtonoise).Normalisationwasperformedbydividingtheindividualcellcountsbytheirindividualcellarea,similartoSalasetal’sbenchmarkingpaper38.Afterfiltering,wehavemergedthedataofthesamplesper“scale”parameter,resultinginthreecompositeAnnDataobjects.Asthesampleswerescreenedwithtwodifferentpanelsofgenes,weobtainedthreeAnnDataobjectspergenepanelrespectively,totallinginsixcompositeAnnDataobjects.
Similartosingle-cellRNA-sequencingdata,thecellcountsfollowedannegativebinomialdistribution,thereforewehavedeemeditappropriatetousescVIforbatchcorrection,amethoddevelopedforscRNA-seq82.Wecorrectedfortheeffectof“patient”,asthisseemedtocausethemainbatcheffectsacrosssamplesbasedontheUMAPscreatedwithrawexpressioncounts.BatchcorrectionwasperformedforeachcompositeAnnDataobjectseparately.
Celltypeandcellsubstateannotation.InitiallydataintegrationandcelltypelabeltransferfromanannotatedplaquescRNA-seqdataset(TabulaSapiens,vasculardataset95)wereattemptedwithvariousmethods(scVI82,scANVI83,scPoli85).However,asgenerallythenumberofgenesandtranscriptcountspercellweresubstantiallylowerinourdatacomparedtotheTabulaSapiensdataset,noneofthemethodsweresuccessfulatintegratingthem.Instead,weperformedmanualcelltypeannotationbyclusteringourdataandannotatingtheclustersbasedontheirdifferentiallyexpressedgenes.
First,todetecthigh-levelcelltypesreflectingmaincellidentities,weperformedcoarseLeiden-clusteringandperformeddifferentialgeneexpression(DGE)analysisina“onevs.all”mannerusingtheWilcoxon-methodasimplementedinscanpy(two-sidedWilcoxonrank-sumtest,withBenjamini-Hochbergcorrection).Forthedetectionoflow-levelcelltypes,correspondingtopreviouslyknowncellsubtypes,furtherLeiden-clusteringandDGEanalysiswithinhigh-levelcelltypeswasconducted.Toaccountforthebiologicalandtechnicalvariabilityacrosssamplesandpatients96,weemployedapseudobulkapproachwithDESeq284astheDGEmethodinthelatterstep.Withineachhigh-levelcelltypesubcluster,pseudobulksampleswerecreatedbyaggregatingcellspersample.Ontopof“patient”wealsocontrolledfor“condition”,assomepatientshadsamplesfrombothconditions,“control”and“plaque”.Low-levelcelltypesweredeterminedbycomparingtheclustermarkergene(two-sidedWaldtest,withBenjamini-Hochbergcorrection)andcellsubtypemarkergene,collectedfromliterature34,88,97,98. Of note, despite not having C1Q or HMOX1 in our gene panels, we assigned the labels Mac_C1Q and Mac_HMOX1basedonoverlappinggenesinourpanelsandthetop20highlyexpressedgenesidentifiedforthesecelltypesinthestudybyDibetal34.Differentialgeneexpressionanalysiswasperformedusingbothapproaches,thetwo-sidedWilcoxonrank-sumtestandthetwo-sidedWaldtestwithDESeq2onpseudobulkcounts(controllingfor“patient”and“condition”),withBenjamini-Hochbergcorrectionappliedacrossfinalhigh-levelcelltypes(SupplementaryData9,SupplementaryData10)andlow-levelcelltypes(SupplementaryData11,SupplementaryData12).
Interestingly,duringthehigh-levelcelltypeannotationweconsistentlydetectedacellcluster,whichwashallmarkedbylowtranscriptcountspercell(SupplementaryFig.46D,SupplementaryFig.47D),highexpressionofmultipleimmunecell-relatedgenes,alongwithgenesofothermajorcelltypes.Basedonthesefeatureswedecidedtoannotatethecellclusteras“Inflammatorymixedcells”.Wehypothesizethatthesecellscouldreflectcellfragments,whosecytoplasmwasonlypartlycontainedintheplaneofthehistologysample,althoughtheycanalsobearesultofincorrectcellsegmentation.Althoughitispossiblethatthesecellsrepresentdoublets,weconsiderthisunlikely.PotentialdoubletsarisingfromcellsegmentationwereexcludedusingBaysor’s“max_cluster_frac”parameter,asrecommendedbytheauthors,andthelowtranscriptcountsobservedinthesecellsfurtherargueagainstthedoubletinterpretation.
We furtherrefinedthelow-levelcelltypesintolow-levelcellsubstatesbyLeiden-clusteringwithinlow-levelcelltypesandperformingDGEanalysisusingthepreviouslydescribedpseudobulkapproachwithDESeq2,controllingforfactors“patient”and“condition”again.Annotationofthelow-levelcellsubstateswascarriedoutmanuallybyinspectingthetopdifferentiallyexpressedgeneswithinlow-levelcelltypes.
CelltypeandcellsubstateannotationwasperformedseparatelyforeachcompositeAnnDataobject.
Annotation of sample subregions.Theannotationofsamplesubregionsisahistologicalanalysis,usingestablishedOxfordPlaqueStudies18–20andAHAcriteria73–76andwasperformedbyanexperiencedcardiovascularpathologist.ThebordersoftheannotatedsubregionswereinitiallydeterminedontheHEstainedslidesbyapathologyexpert,thentransferredtotheDAPI-stainedslidesmanuallybydrawingthedeterminedsubregionbordersbyhandinXeniumExplorer.BordercoordinateswereextractedfromXeniumExplorerandconvertedtopolygons,whichthenwereusedtoassignindividualcellstotheirrespectivesubregionsandtocalculatethesubregionareas.Normalisedsamplesubregionareaswerecalculatedasthepercentageofthewholesamplearea.Someofthesamplestoreorpartlydisintegratedduringthepreparationprocess(controlsamplesofpatient1&2,plaquesamplesofpatient3,4,7,8),inthesecasesonlyanapproximationofthedamagedregion’sareawaspossible.AllsampleswithregionsannotatedaredisplayedinSupplementaryPDF5.
SelectionoffinaldatasetbasedonBaysor’s“scale”parameter.Thepreviousstepsresultedin3×2annotateddatasets:onedatasetpereachofthethreeBaysor“scale”parameters,forbothgenepanels.Inthefollowingwedescribetheinsightsregardingtheimpactofthescale parameteronthecellsegmentation.ComparisonsinSupplementaryFig.48A,BandSupplementaryFig.49A,Bwerecarriedoutbyapplyingpairedtwo-sidedt-testswithBenjamini-Hochbergcorrection,whereasforthecomparisonsinSupplementaryFig.48C,DandSupplementaryFig.49C,Dunpairedtwo-sidedt-testswereused.
- “Scale” of 5μm:Segmentationwithscaleparameterof5μmresultedinahighernumberofcells(SupplementaryFig.48E,SupplementaryFig.49E)withlowertranscriptandcountpercell(SupplementaryFig.46,SupplementaryFig.47),especiallyforvascularsmoothmusclecells(VSMCs),butalsoasignificantlyhigherpercentageoftranscriptswereconsideredasnoisecomparedtosegmentationswiththelargerscaleparameters(SupplementaryFig.48A,SupplementaryFig.49A.WefailedtodetectVSMC-modulinbothpanels(SupplementaryFig.50B,SupplementaryFig.51B),inthepanel2datasetcellsubtypesMac_C1QandMac_S100A8+clusteredtogetherinamixedcluster,whereasMac_HMOX1couldnotbeannotatedatall(SupplementaryFig.51B).Plottingthesegmentedindividualcellswiththelargestcellsizesshowedmostlyonecellnucleuspercellandnoobviousduplicates(SupplementaryFig.52,SupplementaryFig.53).
- “Scale” of 15μm:Thelowestpercentageoftranscriptsweredeemedasbackgroundnoisecomparedtotheothertwoscaleparameters(SupplementaryFig.48A,SupplementaryFig.49A),alongwithtranscriptandgenecountspercellbeingreachingthehighestrange(SupplementaryFig.46,SupplementaryFig.47).Thesamelow-levelcelltypeswerefoundacrossbothpanels,amongstcellsubtypeVSMC_modul,whichwasmissingwithascaleof5(SupplementaryFig.50F,SupplementaryFig.51F).However,thecellsizedistributionofsomeVSMCs(VSMC_fb-like,VSMC_contr,VSMC_inflamed)andMacrophagesubtypes(Mac_TREM2hi,Mac_HMOX1,Mac_C1Q)segmentedwithscaleparameter15showedunrealisticallylargecellsizes(>2000μm2)forthequantile75%andabove(SupplementaryFig.48F,SupplementaryFig.49F),alongwithhighermedianvaluescomparedtosmallerscaleparameters.Individualassessmentoftheselargecellsclearlyshowedmultiplecellssegmentedtogetherasone,especiallyforVSMCandformacrophagesubtypesaswell(SupplementaryFig.54,SupplementaryFig.55).
- “Scale” of 10μm:Noisepercentagewasinbetweenthenoisepercentageofscales5μm and15μm(SupplementaryFig.48A,SupplementaryFig.49A).WewereabletodetectcellsubtypeVSMC_modul,andtheannotatedlow-levelcelltypeswereconsistentwiththecellsubtypesfoundwithscaleparameter15,exceptforMac_ITGAXhi,whichwasnotdetectedinpanel2(SupplementaryFig.50D,SupplementaryFig.51D).TheVSMCandmacrophagecellsubtypeshavingextremelylargecellswithscaleparameter15alsoshowedhighercellsizevaluesinthequantile95%andabove,howevermostofthesecellsweresmallerthan1500μm2(SupplementaryFig.48F,SupplementaryFig.49F).InspectionofindividualVSMCsubtypecellswiththehighest5%cellsizesrevealedthattheywereoftenmultiplespindleshapedcells,segmentedtogetherinoneroundcell.Amongthelargestmacrophageshowever,alongwithduplicateswesawanincreasednumberoflarge,roundmacrophages,thatwerecorrectlysegmentedbyBaysor(SupplementaryFig.56,SupplementaryFig.57).
Althoughsignificantdifferencespresentedinthepercentageoftranscriptslabelledasnoiseacrossthethreeparameters(SupplementaryFig.48A,SupplementaryFig.49A),wedidnotobservestrongsignificantdifferencesinthepercentageoftranscriptssurvivingfiltering(SupplementaryFig.48B,SupplementaryFig.49B).Takentogether,thescaleparameterof5μmwasmoresuccessfulinsegmentingthespindleshapedVSMCcellsbutseemedtofailatcapturingthelargermacrophagesubtypes,whereasscale15μmtendedtoresultin2ormorecellssegmentedasoneforspecificmacrophageandVSMCsubtypeswithlargersizes.Therefore,wedecidedtoselectthedatasegmentedwiththescaleparameterof10μm,whichconstitutedamiddlegroundbetweensegmentinglargerroundmacrophagesandelongated,spindleshapedVSMCs.
SinceBaysorcanalsoestimatethe“scale”parameterforeachsample,individualvaluescanbeusedforsegmentingmultiplesamples.Werecommendcheckingforoutliersintheestimatedvaluesandconsideringthetrade-offbetweenMacrophageandVSMCcellsbasedontheresearchquestion.
HE-stainedslidesofindividualsamples,alongwiththespatialdistributionofcellsfromthefinaldataset,colouredbylow-levelcelltypesandsubregionalassignments,arepresentedinhighqualityformatinSupplementaryPDF1&2.
Assessmentofunassignednoisetranscriptlevelsacrossconditionsandsubregions.WealsoinvestigatedthedifferenceweobservedbetweenplaquesandcontrolsamplesinthelevelofnoisetranscriptsthatwereunassignedtocellsbyBaysor(SupplementaryFig.48C,D,SupplementaryFig.49C,D),utilisingthesubregionalannotationpolygons.Wehypothesizedthatthepresenceofnecroticcoresmayexplainthehigherpercentageofnoiseinplaques,astheygenerallytakeuplargeareasofplaques(SupplementaryFig.58A),andcontainvastcell-poor,celldebrisfilleddomains,whichmayrepresentachallengeforcellsegmentation.NCscontainedonaveragethemosttranscriptsinplaques,bothassignedtocellsaswellasclassifiedasnoise(SupplementaryFig.58B,SupplementaryFig.58C),alsodisplayingthehighestrelativecontributiontotheoverallnoisedetectedamongtheplaquesubregions(SupplementaryFig.58D,bluedots).Thiseffectemergesasaconsequenceoftheirlargesize,possiblycontributingtothedifferenceinnoisepercentagesacrossplaquesandcontrolsamplesseenearlier.Seeingthislargeimpactofthenecroticcoreonnoise,wewerewondering,ifthisisonlycausedbytheNChoardingmostofthetranscripts,ordoesBaysorassignrelativelymoretranscriptsasnoiseintheNCascomparedtoothersubregions?Toassessthis,wecalculatedtheratioofnoisetranscriptswithineachsubregionandfoundNChavingsignificantlyhigherrelativepercentageofnoisecomparedtofibrouscapandmediasubregions,althoughnottointima(SupplementaryFig.58E,bluedots).Thissuggeststhat,inadditiontoitssize,thecompositionoftheNCmayalsocontributetonoise.Wealsocomparedtherelativenoiseofcell-richsubregionsintimaandmediabetweenconditionsandfoundcontrolsamplesdisplayinglowernoisepercentagesintheirmediacomparedtoplaques(SupplementaryFig.58E).AllcomparisonspresentedinSupplementaryFig.58E(withinplaqueandacrossplaque-controlsubregions)wereperformedbyapplyingunpairedtwo-sidedMann-Whitney-Utests,andcorrectedformultipletestingwiththeBenjamini-Hochbergmethod,yieldingfalsediscoveryrates(FDR).Theseobserveddifferencesintherelativenoiseratioswithinsubregionsmaybecausedbyunknowntechnicalfactorsinthepipeline(i.e.,Xeniumreader,Baysor),orcouldalsoarisefromplaquesgenerallycontainingmorecelldebris,eveninregionswithnoobviousnecrosis,likethemedia.Analysiswithalargersamplesizecouldfurtherconfirmthesefindings.
Differentialgeneexpressionanalysis.Fordifferentialgeneexpression(DGE)analysisweagainutilisedthepseudobulkapproach:wecreatedpseudobulksamplesbyaggregatingcellsofthesamesamplewithinlow-levelcelltypesandsubsequentlycarriedoutDGEanalysisformultipleconditionvariablesusingDESeq2,alwayscontrollingfor“patient”asfactor.DifferentialgeneexpressionwasperformedwithDESeq2’sdefaultparameters,yieldingrawtwo-sidedp-valuesfromtheWaldtest,whichwerethencorrectedformultipletestingusingtheBenjamini-Hochbergmethodandreportedasfalsediscoveryrates(FDR).
Cells originatingfromintraluminalthrombiwereexcludedfromDGEanalysis.Touncomplicateouranalysis,weusedsubsetsofdataaccordingtoourscientificquestions:forthecomparison“Plaquevs.Control”onlysamplesfromPatients1-4wereincludedastheyhadmatching“Control”and“Plaque”samples,allothercomparisonswereperformedacrossthe12“Plaque”samples(Patient1-12).ForPatients1-4withavailablecontrolandplaquesamples,wealsoperformedDGEanalysisincludingtheirmedia-derivedcellsonly,duetothemediabeingpresentandproposedtoberelativelysimilarinbothcontrolandplaquesamples.DGEresultscanbefoundinSupplementaryData2,SupplementaryData3andSupplementaryData13andSupplementaryData14.
Clusteringofsamples-morphologicalclusters.First,wecalculatedthecellfractionpercentageofthedifferentcelltypesubstatesnormalisedwithinthesample.Thisenabledustocapturefinerdifferencesincellularcompositionacrossallplaques.Next,weperformedhierarchicalclusteringofthesamplesbasedonthesenormalizedcelltypefractions,usingEuclideandistanceandWard’s-linkagemethod.Thisisanunsupervisedclusteringmethod,anddoesnotrelyonapredefinedstatisticaltestfordetermininggroupassignments.Cellsoriginatingfromintraluminalthrombiwereexcludedfromboththecellfractioncalculationsandclusteringsteps.
Thefinalnumberofclusters(4)wasselectedbasedonclusterrobustness(Fig.5)andfurthersupportedbyvisualinspectionofthesamples(Fig.6).BasedonthedendrogramstructureshownatthetopoftheclustermapinFig.5,selectingfourclustersprovidedabiologicallymeaningfulandrobustseparationbetweengroupsofsamples,whileatthesametime,sampleswithineachclusterdisplayedconsistentpatternsintheirspatialmorphology(Fig.6)asdescribedinthe“Morphologicalplaqueclustering”subsectioninResults.Theseconsiderationstogetherguidedourchoiceoffourclusters.
Weonlyincludedplaquesamplesintheclustering,aswesawmoresignificantheterogeneityamongthesesamplescomparedtothecontrols.Weusedthecelltypesubstatesderivedfrompanel1genes,asthispanelresultedin~20%morecellsthanpanel2,possiblygivingushigherresolutionlookintotheplaquecompositions.
Wealsoperformedhierarchicalclusteringofsamplescontainingintraluminalthrombi,basedonthelow-levelsubstatecellfractionswithincellsdetectedineachthrombus(SupplementaryFig.59).Comparedwiththeoveralldataset,thrombicontainedfewercellsubstates,showedlessvariationacrosssamples,andwerepredominantlycomposedofneutrophils,inflammatorymixedcells,andCD8effectorT-cells.Theendothelialcellsdetectedwithinintraluminalthrombimostlikelyoriginatedfromintimallayerlacerationsintroducedduringsamplepreparation.
Correlationanalyses.Associationsbetweenregionareasandpatientmetadata(SupplementaryFig.12,SupplementaryFig.14–17),cellfractionoflow-levelcellsubstateandmetadata(SupplementaryFig.22,SupplementaryFig.23),aswellascorrelationbetweencellsegmentationmetrics(SupplementaryFig.9)wereassessedtoderivepossiblebiologicalinsightsandtoidentifypotentialconfoundingfactors.ForassociationsbetweennumericalvariablesinvolvingonlyplaquesampleswefittedalinearregressionandreporttheR-squaredvalue,therawtwo-sidedp-valueofthelinearfitandafterBenjamini-Hochbergcorrectionformultipletesting,wealsoreportthecorrectedp-valueofthefitasthefalsediscoveryrate(FDR)(SupplementaryFig.12,SupplementaryFig.14–17,SupplementaryFig.22,SupplementaryFig.23).Inthecaseofcellsegmentationmetricassociations,wecomparedallsamplesincludingsomeplaquesandcontrolscomingfromthesamepatients,thereforewefittedalinearmixedmodelwithpatientshavingrandominterceptstohaveanunbiasedestimateoftheassociation(SupplementaryFig.9).Herewereporttherawtwo-sidedp-valuesandtheBenjamini-Hochbergcorrectedp-valuesasFDR.Toinvestigatetheassociationsbetweennumericalandcategoricalvariablesweconductednon-parametricWilcoxon-testsandreporttheirrawtwo-sidedp-valuesalongwiththeBenjamini-Hochbergcorrectedtwo-sidedp-valuesasFDR(SupplementaryFig.12,SupplementaryFig.14-17,SupplementaryFig.22,SupplementaryFig.23).
Cellneighbourhoodclustering.Wehypothesizedthatalthoughthesamplesshowsubstantialheterogeneityintheirmorphology,assessingthelocalneighbourhoodsofcellscouldrevealsimilarspatialdomainsacrosssamples.Suchdomainsmightreflectsharedbiologicalprocessesorfunctions,potentiallyleadingtoabetterunderstandingofplaqueformation.Unlikehistologicalsubregions,whichrepresentwell-establishedanatomicalcompartmentsofplaques(i.e.,lumen,intima,media,fibrouscap,necroticcore),wesetouttoidentifyinanunbiasedmannerpatternsinlocalcellneighbourhoods.Weappliedunsupervisedclusteringoftheseneighbourhoodstogenerateclustersrepresentingspatialdomainsthatarecharacterisedbysimilarcellularcomposition.Becausetheapproachreliesonlocalneighbourhoodsandunsupervisedclustering,theresultingclustersareindependentofpredefinedsubregionboundaries.Todothis,foreachcellwedefineditscellneighbourhoodasacirclewitha100µmradiusarounditusingsquidpy’s(1.4.1)99spatial_neighbors()function,andcountedthenumberofeachcelltypesubstatewithintheneighbourhood.Thisresultedinanneighbourhoodmatrixwheretherowscorrespondtoallcellsacrossthesamplesandthecolumnscorrespondtoallcelltypesubstatesinthedata.Cellneighbourhoodmatricesweredevisedforeachpanelofgenesseparately.Wesettledon100µmasourneighbourhoodradiusbecauselowerradiusestendedtoresultincellswith0neighbours,whilehigherradiusneighbourhoodsstartedtocapturetheuniquecompositionoftheindividualsamples.InthefinalstepweperformedLeiden-clusteringofthecellneighbourhoodmatricesusingscanpy’spp.neighbors()andtl.leiden()functions,resultingincellneighbourhoodclusters(SupplementaryFig.60,SupplementaryFig.61).Toprovidestatisticalsupportfortheclusters,weperformedpermutationtestingusingneighbourhoodenrichmentscore,ametricdefinedasthetotalnumberofuniquecell-celladjacencyedgesbetweensubstatepairswithinacluster.Wehypothesisedthatifourclustersindeedcapturedistinctspatialpatterns,bycalculatingtheneighbourhoodenrichmentscore,wecouldidentifywhichcellsubstatepairsaresignificantlyenrichedineachother’sneighbourhoodswithinacluster.
First,wecreatedabackgrounddistributionofcellneighbourhoodmatricesbyrandomlyshufflingcellsubstatelabels10,000times.Becausessamplesexhibitedlargedifferencesincellsubstatecomposition,weperformedlabelshufflingwithinsamples.Thiswasnecessarysinceshufflingacrosssampleswouldhaveviolatedtheexchangeabilityassumptionofpermutationtesting100andcouldhaveledtoinflatedp-values.Next,forneighbourhoodenrichmentanalysis,weconstructedspatialneighbourhoodgraphsbyconnectingcellswithina100µmradius,witheachedgerepresentingapairofcellswithinthisdistance.Foreachneighbourhoodcluster,werestrictedthegraphtothecellsassignedtothatclusterandcalculatedneighbourhoodenrichmentasthetotalnumberofuniqueundirectedadjacencyedgesconnectingsubstatepairswithinthecluster-specificgraph.Thiswasagaincomputedforboththeobservedandrandomlyshuffledsubstates,resultinginoneenrichmentscorematrixperclusterfortheobservedcaseand10,000matricesperclusterforthepermutationcase.Empiricalleftandrightp-valueswerecalculatedusingtheMonte-Carlop-valuemethodasdescribedbyPhipsonetal101:p_left=l+1/N+1 (1); p_right=r+1/N+1 (2) wherep_leftandp_rightaretheleftandrightempiricalp-values,listhenumberoftimesenrichmentscorefromthepermutationwaslowerthantheenrichmentscoreoftheobservedcase,risthenumberoftimesenrichmentscorefromthepermutationexceededtheenrichmentscoreoftheobservedcase,andNisthenumberofpermutations.
Inaddition,Agresti-Coullconfidenceintervals(CI)102oftheempiricalleftandrightp-valueswerealsocalculated.Two-sidedempiricalp-valuesandtheirAgresti-Coullconfidenceintervalswerecalculatedusingtheformulaappliedinscipy.stats.permutation_test()function(https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.permutation_test.html):p_two-sided=min(2×p_left,2×p_right) (3); CI_two-sided=min(2×CI_left,2×CI_right) (4)
Rawempiricaltwo-sidedp-valueswerethencorrectedformultipletestingbyapplyingtheBenjamini-Hochbergcorrectionmethod,andreportedasfalsediscoveryrate(FDR).SignificancewasdefinedasFDR<0.05.Permutationtestresultsforbothpanel1(SupplementaryPDF6,SupplementaryData15)andpanel2(SupplementaryPDF7,SupplementaryData15)neighbourhoodclustersidentifiedsignificantlyenrichedordepletedsubstatepairsineveryclustercomparedtotherandombackground.Basedontheseresultswebelieveourneighbourhoodclusteringapproachwasabletoidentifyspatialdomainswithsignificantlydifferentcompositionalpatterns,evenacrossmultiplesamples.Highresolutionspatialscatterplotsofallsampleswithcellscolouredbylow-levelcelltypesandneighbourhoodclustersarecontainedinSupplementaryPDF8forpanel1genes,andinSupplementaryPDF9forpanel2genes.Toshowcasethemacrophageneighbourhoods,weselectedplaquesfrompatients4and8forinclusioninthemainfigures,asthesecontainedmacrophageclusterswiththemostillustrativespatialpatterns.Inaddition,thesesamplesalsocontainedEC_neovesselneighbourhoods,whicharealsohighlightedinthemaintext.ToillustrateVSMCneighbourhoods,weselectedsampleswithacell-denseandintactmedia(plaquesfrompatients7and10,andthecontrolfrompatient3)forthemainfigures.
Gene expressiondifferencesofmacrophagesubstatesinneighbourhoodclusters2and4.AsthespatialdistributionofthefourmacrophagesubstatesMac_C1Q,Mac_HMOX1,Mac_TREM2hi_CCL18hi,Mac_TREM2hiwithinthepanel1neighbourhoodclusters2and4showedasomewhatorderedpatternfromlumentothenecroticcoreinspecificplaques(Patients4,7,8,10),wehypothesizedthatagradualexpressiontransitionmayexistinthedirectionofMac_HMOX1/Mac_C1Q→Mac_TREM2hi_CCL18hi→Mac_TREM2hi,astheycellexpressionsarealsosimilar(SupplementaryFig.35B).Therefore,wesetouttoidentifygenescharacterisingthesesubstatesandalsotoassessthesubtledifferencesintheirexpressionlevels,aimingtoderivebiologicalinsights.Wefilteredthedatatoincludecellsfromthespecifiedplaques(Patients4,7,8,10)thatbelongtothethreemacrophagesubstatesandPanel1neighbourhoodclusters2and4.Toaccountforthesequencingdepthdifferencesacrosssamples,weadditionallymin-maxnormalisedthearea-normalisedexpressionofcellswithinsamples,resultinginexpressionlevelsbetween0and1.Nextweaveragedthemin-maxnormalisedexpressionsofallgeneswithincellsofeachsubstate,andselected20geneswiththehighestvarianceinmeanexpressions.WeclusteredthegenesbytheZ-scorestandardisedmeanexpressionswithinthesubstate(SupplementaryFig.35A),todetectgeneswitheventualgradientsinexpressionlevelsacrosssubstates,alsoplottingthegeneexpressionsofindividualcells(Fig.8B).Toteststatisticalsignificanceofgeneexpressionleveldifferenceswesawinthepreviousplots,weperformedpairwisedifferentialgeneexpression(DGE)analysisbetweenthethreehypothesized“endpoint”substatesMac_C1Q,Mac_HMOX1andMac_TREM2hiinoursubcohort(SupplementaryFig.35C-E,SupplementaryData16).DGEwasperformedinthesammannerasdescribedintheprevioussection“Differentialgeneexpressionanalysis”.
Spatialcorrelation ofmacrophagesubstates.Toassessthestatisticalsignificanceofthespatialdistributionpatternsseeninthethree“endpoint”macrophagesubstates(Mac_C1Q,Mac_HMOX1,Mac_TREM2hi),wecomputedthespatialcorrelationbetweenpairsofthesesubstates,usingtheLocalBivariateMoran’sIstatistic103.Foreachpairofsubstates,wefirstbinarizedcellidentitiesacrossallspatialcoordinates(1=presence,0=absence),wethenconstructedaK-nearestneighbours(KNN)spatialweightsmatrixusingEuclideandistancesbetweenallcells(k=30),withrow-standardizationapplied.LocalBivariateMoran’sIwascalculatedforeachcell,measuringtheassociationbetweenthepresenceofonecelltypeandthespatialdistributionofanotherinitsneighbourhood.Therawone-sidedempiricalp-valuesoflocalassociationsignificancewerecorrectedformultipletestingusingtheBenjamini-Hochbergmethodpersubstatepairpersample,resultinginfalsediscoveryrates(FDR).Cellswithsignificantspatialcorrelation(FDR<0.05)werevisualizedbymappingsignificantspatialclusters(e.g.,bothsubstateshigh,onesubstatehigh-onesubstatelow,orbothlow)ontospatialscatterplotsofsamples(SupplementaryFig.30-34).
Ofnote,macrophagecellswithpanel2datasetexpressedmoreTFs(totalnumberofTFs:12;meannumberofTFs/gene:5.3),thanpanel1(totalnumberofTFs:5;meannumberofTFs/gene:1).
Statistics&reproducibility
Nostatisticalmethodwasusedtopredeterminesamplesize.Nodatawereexcludedfromtheanalyses.Theexperimentswerenotrandomized.TheInvestigatorswerenotblindedtoallocationduringexperimentsandoutcomeassessment.
Reportingsummary
FurtherinformationonresearchdesignisavailableintheNaturePortfolioReportingSummarylinkedtothisarticle.
Dataavailability
DatafromourownscRNA-SeqandsinglecellspatialtranscriptomicsexperimentshavebeendepositedinNCBI’sGeneExpressionOmnibus.Thepreviouslygenerateddatasets,thatarealreadypublished,areaccessiblethroughGEOSeriesaccessionnumberGSE247238(scRNA-Seq,patients1-10)26,55,90,105orinthismanuscriptinSupplementaryData17(rawbulk-seqcountmatrix:https://zenodo.org/records/17454209)andSupplementaryData1826,90,106–108.Newlyacquireddatasetsforthisstudyare:GSE294291(scRNA-Seq,patients11–17)andGSE294466(Xeniumspatialtranscriptomics).AnadditionalexternalscRNA-Seqdatasetofhumancarotidplaques51usingasimilarset-upofplaquevs.adjacentcontrolwasdownloadedfromtheGEOrepository(GSE159677).ApreprocessedandmanuallyannotatedversionoftheGSE159677dataset,alongwiththebatch-correctedbulkRNA-seqcountmatrixusedinthisstudy,canbedownloadedfromZenodo:https://zenodo.org/records/17521879.WealsoprovidethefinalprocessedversionsofourscRNA-seqandXeniumdatasetsasresourcestothecommunity.scRNA-seqdataareavailableat(https://zenodo.org/records/17520898);Xeniumdatasets,alongwithsubregionannotationandsampleareasareavailableat:(https://zenodo.org/records/17526248).InteractiveSupplementaryFig.7and8,alongwiththeirlegends canbedownloadedfromZenodo:(https://zenodo.org/records/17523511).
Codeavailability
Codeisavailableat:(https://github.com/MendenLab/Carotid_artery_plaque).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors integrated spatial, single cell and bulk transcriptomics of carotid plaque, identifying inflammatory smooth muscle cell subtypes and localizing regions of neovascularization and macrophage activity. The authors applied a clever deconvolution algorithm to correct for immune cell bias and estimate true cell type proportions in bulkRNA-Seq data using scRNA-Seq-derived labels. The focus is on VSMCs and macrophages. There are several new and interesting findings:\n1. TREM2 macrophages are not intermingled with C1Q macrophages, which sit below the endothelium\n2. The authors convincingly discovered that immune cells are overrepresented in scRNA-Seq experiments\n3. The authors dicovered that there was no difference between symptomatic and asymptomatic plaques. In many studies, this is the only distinction between plaques provided\nSpecific comments\n1. Previous work has provided systematic pathological classifications (Pasterkamp, Virmanu). Please discuss how your findings relate to these.\n2. Note that endarerectomies do not provide access to the adventitia, where many immune cells reside. This may explain why the immune cell content is so low.\n3. Recently,TLOs were discovered in human plaque1 using methods similar to yours. Did you see any indications of TLOs?\n4. Define CEA (probably the surgery).\n5. The samples were macroscopically assessed and dissected into advanced atherosclerotic (plaque) and adjacent healthy/early lesion (control) tissue. It is not clear that the “control” is really healthy.\n6. Since some samples had thrombus: what cells were in the thrombus?\n7. FC thickness set at 200 μm: This is not easy to measure, because the sections may be at different angles, and cap thickness is not uniform. How reliable is this measure?\n8. PI3K has many isoforms. Which ones are DE?\n9. CD79A is not an immune regulator, but a B cell marker\n10. Page 14: It is not clear which datasets are new (unpublished). Clearly distinguish between new and reanalyzed data.\n11. Figure 2H: The immune cells are invisible. Either make an extra panel with a different scale, or use a log scale.\n12. Figure 5A: numbers too small (not legible)\n13. Figure 6A: show directionality (which genes were up in smokers?)\n14. Page 4: “unpaired (n=232, n=65) and paired (n=42) designs, as well as stable (n=97) and unstable (n=126) plaques”. Some may be overlapping. What was the total number of samples?\n15. The authors often talk about “carotid plaques”, but much of the resected material is not plaque but media. Maybe call them endarterectomy specimens.\n16. SMCs, ECs and fibroblasts show up much better in single nuclei sequencing. Many studies used this approach, mitigating the bias. Please mention this and cite some studies.\n17. The authors write that “VSMC-rich media from our CEA specimens (Fig 1A), a structure that has likely not been analysed in previous studies“. Was there a difference in gene expression in the media of adjacent normal vs atherosclerotic artery tissue? Paired analysis should be possible.\n18. “Four patients provided paired samples”. Were these bilateral endarterectomies, with one side worse than the other? If so, what was the time elapsed between the two surgeries?\n19. Page 6: ”Correlation analysis revealed a strong inverse relationship between NC and lumen” suggests that the sum of the two may be constant and, by inference, the fibrous cap area % must also be constant.\n20. The definition of controls is unclear. A sample with 45% intima is not a control. It seems that the authors call samples without NC controls.\n21. It is interesting that there was no difference between symptomatic and asymptomatic plaques. This should be emphasized. In many studies, this is the only distinction provided.\n22. Some patients showed “two NCs on each side of the lumen”. This is interesting, but not conclusive. They may be connected in another plane.\n23. In several scRNA-Seq studies and in at least 2 reviews, about half of the foam cells are SMC-derived and half are macrophages. There are two sets of macrophage foam cells and at least 3 other macrophage sets that are not foamy. The information obtained in the present study should allow these assignments.\n24. Some of the correlations with clinical and demographic parameters are not statistically analyzed. Also, mentioning a single SNP in the smoking section is not appropriate. SNPs were not tested in this study.\n25. Please define the difference between neighborhoods and subregions.\n26. Figure 2E: the choice of colors is unfortunate. Impossible to see the cluster calling.\n27. Figure 3D: what are “inflammatory mixed cells”? Doublets?\n28. Figure 3E: Age seems to be a driver for necrotic core. This should be emphasized.\n29. Figure 5: What statistical test was used to assign the morphological clusters? How do we know there should be 5 clusters? Or is this arbitrary?\n30. Figure 7C: How were patients 4 and 8 selected?", "reviewer_2": "Overall Summary:\nPauli et al. present their investigation seeking to integrate multiple transcriptomic methodologies (bulk RNA seq, scRNAseq, and spatial transcriptomics) as a means to molecularly define carotid atherosclerotic disease. Using arterial and atheroma samples collected during CEA, the investigators attempt to elucidate a spatially-resolved molecular architecture underlying plaque development. The authors highlight unique plaque morphologies linked to different cell type compositions, disease heterogeneity (i.e. symptomatology), and clinical demographics (i.e. smoking, etc) with inference towards plaque development and stability. Notably, they uncover spatial associations between HMOX1+ macrophages as they potentially transition to foam cells (TREM2+).\nWhile the manuscript is well written and demonstrates technical rigor, the most important limitation is a lack of methodologic novelty or conceptual impact to our molecular understanding of atherosclerosis. Notably, HMOX1 and TREM2 macrophages are well established in the progression of atherosclerosis (PMID: 39731912 – HMOX1, 38974464 – TREM2) and modulation of macrophage polarity (PMID: 30282830 – HMOX1, 38646596 – TREM2). Furthermore, other spatial transcriptomic investigations of arterial atherosclerosis have recently been published (PMID: 39234691 and 39715784). Considering these limitations, the current manuscript offers incremental advancement.\nOther Criticisms include:\n1. Association with symptomatology / stability is uncertain. One area of potential significance was the assertion that the current findings molecularly define symptomatic (vs. asymptomatic) carotid artery stenosis. As presented, however, the data do not clearly support these claims. Figure 5 provides an example of this missed opportunity. A clearer presentation of the results would have been based upon symptom status as opposed to various morphologic clusters.\n2. Surprising lack of DEG in bulk RNAseq. Figure 2c shows relatively few differences in DEG across stable and unstable plaques. This is in contrast to current literature (39234691), yet no explanation is given.\n3. Underpowered study with significant demographic heterogeneity. It is difficult to interpret some of the results due to relatively few patients (N=4 Ctl) and variable clinical demographics.\n4. Biased spatial transcriptomic approach based off Panels 1 and 2. Finally, the spatial transcriptomic analysis used a custom (pre-defined) gene set (~500) for cellular and disease status identification. While this approach is helpful in resolving cellular spatial location, it falls short of providing true spatially-resolved, single cell transcriptomic analysis. This undercuts the novelty.", "reviewer_3": "Summary\nThis manuscript presents a comprehensive spatially resolved transcriptomic analysis of human carotid atherosclerotic plaques. The authors integrate bulk RNA-Seq, single-cell RNA-Seq (scRNA-Seq), and single-cell spatial transcriptomics (Xenium platform) on paired samples from patients undergoing carotid endarterectomy (CEA). Their multimodal approach reveals spatial microenvironments, including distinct VSMC and macrophage subtypes associated with plaque vulnerability, progression, and smoking status. They propose a model of subluminal HMOX1⁺ macrophage transdifferentiation into lipid-handling TREM2⁺ cells. A novel morphological clustering of plaques into “structured” and “chaotic” types is introduced, offering mechanistic insights into spatial immune-VSMC interplay and plaque destabilization.\nMajor Comments\n1. Justification of Technology Choice: The authors employ the Xenium platform for spatial transcriptomics but do not justify its selection over other platforms like 10x Visium HD. A justification regarding resolution, sensitivity, and gene coverage tradeoffs is needed.\n2. Immune Cell Overrepresentation Not Corrected: Despite deconvolution, the overrepresentation of immune cells in scRNA-Seq is not adequately addressed in the spatial analysis. This may bias subregion interpretations.\n3. Smoking Association Requires Caution: The observed association with smoking is only partially supported by bulk and scRNA-Seq data. Acknowledgment of limited reproducibility is needed.\n4. Trajectory Analysis of Macrophages: The proposed HMOX1⁺ to TREM2⁺ trajectory should be validated with orthogonal methods beyond TFvelo and visualized more explicitly.\n5. Subregion Definition and Reproducibility: Manual annotation of regions may introduce bias. Providing inter-rater reliability or uncertainty metrics would enhance robustness.\n6. Sample Size for Spatial Analysis: The n=12 for spatial transcriptomics limits generalizability. Interpretations of clustering and risk associations should be more cautiously framed.\n7. Neighbourhood Cluster Analysis: Statistical support (e.g., permutation tests, spatial proximity significance) for proposed functional zones is missing and should be included.\n8. Clinical Utility and Predictive Value: The translational relevance of findings should be discussed. Can DEGs or spatial patterns serve as prognostic markers for plaque rupture?\nMinor Comments\n1. Terminology Clarification: Terms like 'structured' vs. 'chaotic' could be more quantitatively or mechanistically defined.\n2. Figure Readability: Dense color legends in Figures 6–8 make interpretation difficult. Consider enhancing clarity.\n3. Figure 2 caption: “-log10 value” => “-log10 p-value”"}, "author_response": "Point-by-point response\nWe would like to thank the reviewers and editors for the assessment of our submitted\nmanuscript. We carefully addressed all of the reviewers’ concerns with new data\nanalysis as well as clarifications, which we believe has enabled us to submit a\nsubstantially improved version of our study.\nReviewer #1 (Remarks to the Author):\nThe authors integrated spatial, single cell and bulk transcriptomics of carotid plaque,\nidentifying inflammatory smooth muscle cell subtypes and localizing regions of\nneovascularization and macrophage activity. The authors applied a clever\ndeconvolution algorithm to correct for immune cell bias and estimate true cell type\nproportions in bulkRNA-Seq data using scRNA-Seq-derived labels. The focus is on\nVSMCs and macrophages. There are several new and interesting findings:\n1. TREM2 macrophages are not intermingled with C1Q macrophages, which sit below\nthe endothelium\n2. The authors discovered that immune cells are overrepresented in\nscRNA-Seq experiments\n3. The authors discovered that there was no difference between symptomatic and\nasymptomatic plaques. In many studies, this is the only distinction between plaques\nprovided\nSpecific comments\n1. Previous work has provided systematic pathological classifications (Pasterkamp,\nVirmani). Please discuss how your findings relate to these.\nWe thank the reviewer for this important comment on variations of assessing plaque\nvulnerability statuses. This is obviously a rather debated question, and we believe that\nall current approaches have their advantages but also disadvantages. Pioneers in the\nfield, like the mentioned groups led by Renu Virmani and Gerad Pasterkamp\nrespectively, established their plaque phenotyping based on the available vascular\nspecimens their labs and institutes had and partly still have access to. For the carotid\nplaque specimens in our study, we utilize the Rothwell/Redgrave criteria of the Oxford\nPlaque Study, being explained in more detail below. We made sure to more clearly\nstate that characterization using this criterion for carotid artery disease in our revised\nmanuscript.\nWork by Virmani and colleagues was in large parts performed and influenced\nby postmortem autopsy specimens from different vascular beds, with a focus on\ncoronary artery tissue specimens. According to these autopsy studies, the most\ncommon type of vulnerable plaques are so called thin-cap fibroatheroma (TCFA) with\nthe following key features: fibrous cap thickness <65µm (for coronary arteries and not\ncarotids as used in our study) located over a lipid/cholesterol crystal-rich necrotic core,\nhigh presence of T cell and macrophage infiltrations in close proximity to the fibrous\ncap (plaque shoulder), low numbers of smooth muscle cells, and areas of intraplaque\nbleedings (“Pathology of the Vulnerable Plaque,” 2006).\nThe work by Gerard Pasterkamp offers a broader, multifactorial definition of\nvulnerable atherosclerotic plaques, emphasizing not only histological features but also\nbiological activity, mechanical stress, as well as systemic factors. In studies by his\ngroup, mostly tissue specimens from carotid endarterectomies (CEA) are utilized and\ninvestigated predominantly based on patient symptoms (asymptomatic or symptomatic\nfor stroke, transient ischemic attack, and amaurosis fugax; (Langley et al., 2017; Mocci\net al., 2024)), and more recently based on transcriptomic expression patterns (Mokry\net al., 2022).\nA third approach, which was used in our current study stems from Drs. Rothwell\nand Redgrave and the Oxford Plaque Studies. Here, a crucial fibrous cap thickness of\n>/< 200µm is established exclusively for carotid artery disease, taking into account\ncarotid imaging modalities and clinical outcomes from patients with stroke (Lovett et\nal., 2004; J. N. Redgrave et al., 2008; J. N. E. Redgrave et al., 2006). Based on our\nexperience and work published by us in previous studies (Bashore et al., 2024; Fasolo\net al., 2021), this carotid artery specific assessment of plaque vulnerability suits our\napproach of combining histomorphological plaque characterization, integrative multi-\ntranscriptomic profiling (with bulk, single cell, and spatial transcriptomics) and\ncorrelation with clinical patient information (symptoms, risk factors, accompanying\ndiseases) best.\nWe added the following statement in the Results section to clearly state our approach\nat the beginning of the manuscript:\n“Human carotid plaques from CEA patients were selected based on conclusive\nmorphology, including remnant media and macroscopically assessed and\ndissected into advanced atherosclerotic (plaque) and adjacent early lesion\n(control) tissue. Plaques are characterized based on the Oxford Plaque Studies\n(Lovett et al., 2004; J. N. Redgrave et al., 2008; J. N. E. Redgrave et al., 2006).\nBoth pieces of this paired unique setup underwent bulk- and scRNA-Seq.”\nWe furthermore clarified the Methods section in this regard:\n“The histological classification of carotid plaques was performed following the\nguidelines established by the American Heart Association (AHA(Stary, 2000;\nStary et al., 1992, 1994, 1995)) and the Redgrave & Rothwell criteria (Oxford\nPlaque Study)(J. N. Redgrave et al., 2008) for plaque stability.”\n2. Note that endarterectomies do not provide access to the adventitia, where many\nimmune cells reside. This may explain why the immune cell content is so low.\nThe reviewer is correct that plaques from carotid endarterectomies (CEA) lack\nadventitia, as the surgical removal of the lesion must exclude the adventitia remaining\nin the patient. Thus, CEA plaques contain the media layer, which leads to a higher\nSMC content in general and therefore as suggested by the reviewer likely a lower\nabundance of immune cells. Of note, previous prominent studies have used FACS\nsorting of CD45+ cells to analyze human carotid lesions, which fueled the assumption\nthat human plaques mainly consist of immune cells (Depuydt et al., 2023; Dib et al.,\n2023; Fernandez et al., 2019). This is the main reason why we performed single cell-\nbulk RNA-Seq deconvolution (utilizing BayesPrism) in our manuscript, enabling us to\novercome the immune cell bias of applying single cell RNA-Seq in human plaques.\nImportantly, our finding is confirmed by our own spatial transcriptomics (Xenium)\nanalysis, in concordance with very recent literature (Bleckwehl et al., 2025; Traeuble\net al., 2024).\nWe clarified this matter in our Introduction:\n“The intima-media layer, rich in VSMCs and fibroblasts, represents the pre-\natherosclerotic arterial wall. The adventitia is not excised during the surgical\nprocedure and is retained in the patient. All three substructures…”\n3. Recently,TLOs were discovered in human plaque1 using methods similar to yours.\nDid you see any indications of TLOs?\nstructures in human atherosclerotic plaques. Nat Cardiovasc Res,\ndoi:10.1038/s44161-025-00639-9 (2025).\nThe TLOs described in the Lai et al. paper were aggregations of lymphoid cells,\ncharacterised by plasma and B cells in particular. As our current study focuses on\nintegration of bulk, single cell and spatial transcriptomics in an unbiased fashion (not\ntargeting specific structures), our Xenium probe panels did not necessarily include\nsufficient marker genes for a more detailed identification of lymphoid cell subtypes. In\naddition, plasma and B cell marker gene counts are comparatively low in relation to\nmacrophage and SMC marker gene counts.\nBut to directly address the reviewer’s comment, we carefully examined whether\nTLOs could still be identified in our plaques. For this, we calculated the spatial\nautocorrelation of plasma cells (Gerdis-Ord G statistic) and assessed the expression\nof four genes overlapping between our panels and the PTLO marker genes reported\nby Lai et al. (CD19, CD3D, MZB1, and XBP1). Plaques from patients 6 and 7\nconsistently exhibited significant plasma cell autocorrelation across both the panel 1\n(Suppl. Fig. 43A,B, see below) and panel 2 (Suppl. Fig. 44A,B, see below) datasets.\nIn addition, upon revisiting our neighbourhood cluster analysis, we found cluster 13 in\nboth panel 1 and panel 2 to be enriched in plasma cells (Suppl. Fig. 28A,B). These\nclusters were detected only in the plaque of patient 6 and showed strong spatial\noverlap with the plasma cell aggregation in the necrotic core of this sample (Suppl.\nFig. 45, see below).\nHowever, among the four marker genes, only MZB1 and XBP1 showed\nsignificant expression within plasma cells in plaques from patient 6 and 7 (Suppl. Fig.\n43C,D; Suppl. Fig. 44C,D).\nTaken together, the spatial patterns of plasma cells may suggest the presence\nof PTLOs in the plaques of patients 6 and 7, however, these cells expressed only two\nof the four PTLO marker genes. This highlights that, given the limited ability of our\ncurrent gene panels to resolve lymphoid cell subtypes, further clarification of these\nstructures will require future analyses with more dedicated gene panels.\nTo incorporate this analysis, we expanded our Methods section:\n“Spatial autocorrelation of plasma-cells and plaque tertiary lymphoid organ\n(PTLO) marker gene expression\nLai et al. identified aggregations of lymphoid cells in plaques, especially\nenriched in B and plasma cells, as plaque tertiary lymphoid organs (PTLOs) in\na recent publication (Lai et al., 2025). Motivated by these findings, we sought to\nassess the presence of such domains in our samples by calculating the spatial\nautocorrelation of lymphoid cells. For this, we used the Getis-Ord G statistic,\nsuited for “hotspot” analysis in spatial datasets (Ord & Getis, 1995). As our\nXenium panels did not include sufficient marker genes for a more detailed\nidentification of lymphoid cell subtypes (i.e. B-cells), we were only able to test\nthe autocorrelation of plasma-cells.\nFirst, a spatial weight matrix was calculated for each sample using libpysal’s\nDistanceBand() function with a radius of 2500 µm. This larger threshold was\nchosen to allow detection of broader-scale clustering patterns and to reduce\nspurious local pairings, which were apparent at lower thresholds due to the\nrelatively low number of plasma cells. Local Getis-Ord G statistics were then\ncalculated using the G_Local() function from the esda package, with binarized\ncell identities across all spatial coordinates (1 = plasma cell, 0 = non–plasma\ncell). To generate an empirical background distribution, the permutations\nparameter was set to 999. Statistical significance of autocorrelation was\ndetermined by applying Benjamini–Hochberg correction to the empirical p-\nvalues, yielding false discovery rate (FDR) values.\nIn addition, we assessed the expression of four PTLO marker genes present in\nour panels (CD19, CD3D, MZB1, XBP1). Within each sample, normalised\nexpression of the marker genes were compared between plasma-cells and all\nother cell types, using the Wilcoxon rank-sum test. The resulting p-values were\ncorrected with the Benjamini–Hochberg method.\nOur analysis identified two plaques (patients 6 and 7) with consistent and\nsignificant plasma cell autocorrelation across both panels. Results obtained with\npanel 1 genes are shown in Suppl. Fig. 43, and those with panel 2 genes in\nSuppl. Fig. 44.”\nAlong with the Results section:\n“A recent publication identified plasma and B-cell rich areas, interspersed with\nVSMCs and macrophages as plaque tertiary lymphoid organs (PTLOs) (Lai et\nal., 2025). In our samples, spatial analysis revealed significant plasma cell\nenrichment in the plaques of patients 6 and 7 (Suppl. Fig. 43A,B; Suppl. Fig.\n44A,B). In addition, the plasma cell aggregation domain in the NC of patient 6\nwas identified as cluster 13 (Suppl. Fig. 45), enriched in VSMCs, macrophages\nand plasma cells (Suppl. Fig. 28A). Although suggestive of PTLOs, only two\n(MZB1, XBP1) out of four PTLO marker genes (Lai et al., 2025) in our panels\nshowed high expression in plasma cells (Suppl. Fig. 43C,D; Suppl. Fig.\n44C,D), highlighting the need for further clarification with more specialised gene\npanels.\nBy performing neighbourhood clustering for panel 2, we detected similar\nneighbourhoods with similar features in the same plaques, with neighbourhood\n0 and 1 being “VSMC neighbourhoods”, neighbourhoods 2, 5 and 8 being\n“macrophage neighbourhoods”, and neighbourhood 10 being the\n“EC_Neovessel neighbourhood”, and neighbourhood 13 the “plasma cell\nneighbourhood” (Suppl. Fig. 28B,C).“\nWe adapted the Introduction and mention the suggested publication, together with\nother spatial transcriptomics publications in human atherosclerotic plaques:\n“Sequencing technologies, including bulkRNA-Seq, single cell (sc) RNA-Seq,\nand spatial transcriptomics, have helped tremendously to better determine\nmolecular alterations and cellular composition in human diseases (Bleckwehl et\nal., 2025; Campos et al., 2025; Gastanadui et al., 2024; Lai et al., 2025; Longo\net al., 2021; Sun et al., 2023)”\n4. Define CEA (probably the surgery).\nWe apologize for this mistake. Carotid endarterectomy (CEA) is defined as the surgical\nremoval of plaque buildup in carotid arteries. We adapted the Introduction as well as\nthe Methods section accordingly:\nIntroduction:\n“Atherosclerotic plaques from patients undergoing carotid endarterectomy\n(CEA), the surgical removal of plaque buildup in carotid arteries, consist of three\nmajor subregions with distinct cellular and molecular features.”\nMethods:\n“Plaque samples were collected from patients undergoing carotid\nendarterectomy (CEA), the surgical removal of plaque buildup in carotid\narteries.”\n5. The samples were macroscopically assessed and dissected into advanced\natherosclerotic (plaque) and adjacent healthy/early lesion (control) tissue. It is not clear\nthat the “control” is really healthy.\nWe thank the reviewer for this comment. It is correct that the control pieces in this\ncurrent study do not reflect a fully healthy artery, as this tissue still stems from elderly\npatients suffering from carotid artery disease. The ‘healthier’/early or lesser diseased\nlesions are located adjacent to the highly stenotic carotid artery tissue. Controls in our\nwork are defined as the least diseased piece of tissue from the patient, with no signs\nof necrotic core formation (which is the defining hallmark of an advanced carotid artery\nplaque). Intimal thickening as well as early signs of fatty streaks - and in some cases\nimmune cell infiltrates might still be present, as this initiation process occurs already\nvery early in life and is highly present in all (in particular Western) populations (Finn et\nal., 2010; Masawa et al., 1994; Pasterkamp et al., 1998; Velican, 1969; Virmani et al.,\n2000).\nWe carefully screened our revised manuscript again and removed the word “healthy”\nfrom the Results section (1x), as well as from the Methods section (1x) when explaining\nand discussing sample selection and comparison thereof.\nResults:\n“Human carotid plaques from CEA patients were selected based on conclusive\nmorphology, including remnant media and macroscopically assessed and\ndissected into advanced atherosclerotic (plaque) and adjacent tissue\nhealthy/from the same patient (matched paired diseased-control samples).\nPlaques were characterized for stability based on the Oxford Plaque Studies.\nBoth pieces of this unique setup (statistically powerful paired analysis possible)\nunderwent bulk- and scRNA-Seq. In addition, formalin-fixed paraffin (FFPE)-\nembedded samples of neighbouring pieces—maintaining the same paired\nplaque vs. control configuration—were processed for histomorphological and\nspatial transcriptomics analyses\nMethods:\n“We selected the plaque scRNA-seq dataset published by Alsaigh et al(Alsaigh\net al., 2022), which had a similar experimental setup (matched atherosclerotic\ncore - proximal adjacent samples from same patients) to ours (matched\ndiseased - healthy control samples from same patients).”\nFurthermore, we emphasized this unique setup of early/control and advanced lesion\nfrom the same individual patient that allows for statistical extremely powerful analyses\nin the same paragraph of the Results section:\n“Both pieces of this paired unique setup underwent bulk- and scRNA-Seq. Both\npieces of this unique setup (statistically powerful paired analysis is possible)\nunderwent bulk- and scRNA-Seq. In addition, formalin-fixed paraffin (FFPE)-\nembedded samples of neighbouring pieces—maintaining the same paired\nplaque vs. control configuration—were processed for histomorphological and\nspatial transcriptomics analyses (Fig. 1).”\nAdditionally, we clarified our definition of the early lesion in the Methods section:\n“Immediately after removal from the intraoperative situs, plaque samples were\ntransferred to pre-chilled PBS (Merck, Germany) for transportation. The carotid\nplaques are then cut into several pieces. Specifically, for this study, we selected\nthe most advanced (i.e., most unstable) part of the plaque and whenever\npossible also the least diseased part of the plaque for analysis (Suppl. Fig. 3).\nThis unique setup allows for statistically powerful paired analyses of two pieces\nof the same tissue from the same individual patient.”\nSpecifically for this purpose, we also created a new supplemental figure - Suppl. Fig.\n3 - that illustrates the explanations from above:\nFigure Legend: “Suppl. Fig. 3: Carotid Artery Plaque tissue sectioning and\ncharacterisation routine. (A) shows the whole plaque upon surgical removal\n(CEA). (B) shows the sectioned plaque getting processed for molecular and\nhistomorphological characterization. The arrows indicate the most advanced\nand earliest piece of tissue in this example. (C) Histological analysis is\nperformed according to the American Heart Association (AHA) classification\nguidelines(Stary, 2000; Stary et al., 1992, 1994, 1995) using different\nhistological stainings (HE and EvG). FC thickness is evaluated based on the\nOxford Plaque Study ((J. N. Redgrave et al., 2008)) and measured using the\nAxioScan Software (Leica) or QuPath. Here, all four measured positions of the\nFC are below 200µm, therefore this carotid plaque was classified as unstable.”\n6. Since some samples had thrombus: what cells were in the thrombus?\nThis is an important and interesting comment by the reviewer. We excluded cells\noriginating from intraluminal thrombi from our clustering and DGE analyses. This was\nhowever not mentioned in the previously submitted version of our manuscript. We now\nexpanded the respective sections in our Methods part:\n“For differential gene expression (DGE) analysis we again utilised the\npseudobulk approach: we created pseudobulk samples by aggregating cells of\nthe same sample within low-level cell types and subsequently carried out DGE\nanalysis for multiple condition variables using DESeq2, always controlling for\n“patient” as factor. Cells originating from intraluminal thrombi were excluded\nfrom DGE analysis.”\n“First, we calculated the cell fraction percentage of the different cell type\nsubstates normalised within the sample, then performed hierarchical clustering\nwith Euclidean distance and Ward’s-linkage method. Cells originating from\nintraluminal thrombi were excluded from both the cell fraction calculations and\nclustering steps. ”\nFurther, in the revised manuscript, we included the analysis of the five plaques and\none control sample that contained cells detected in thrombi. We also expanded the\nMethods section accordingly and added a new supplementary figure (Suppl. Fig. 59):\n“We also performed hierarchical clustering of samples containing intraluminal\nthrombi, based on the low-level substate cell fractions within cells detected in\neach thrombus (Suppl. Fig. 59). Compared with the overall dataset, thrombi\ncontained fewer cell substates, showed less variation across samples, and were\npredominantly composed of neutrophils, inflammatory mixed cells, and CD8\neffector T-cells. The endothelial cells detected within intraluminal thrombi most\nlikely originated from intimal layer lacerations introduced during sample\npreparation.”\n7. FC thickness set at 200 μm: This is not easy to measure, because the sections may\nbe at different angles, and cap thickness is not uniform. How reliable is this measure?\nWe acknowledge the reviewer’s concern regarding the challenges in measuring fibrous\ncap (FC) thickness due to sectioning angles and intrinsic variability in cap morphology.\nOur human samples originate from patients undergoing carotid endarterectomy, and\nall tissues are processed under a standardized biobank protocol that has been in place\nfor >15 years. We made sure to add extra information into the revised manuscript to\nbetter explain how this is performed.\nWe added a detailed description of the FC measurement to Suppl. Fig. 3\n(Suppl. Fig. 3C), that can be found in your Question #5.\nPlease also see Question #5 for additional details about how we addressed\nchanges regarding our sample setup in the text for your clarification.\n8. PI3K has many isoforms. Which ones are DE?\nAs we used short read sequencing to obtain bulk- and scRNA-Seq data, we cannot\nmake distinctions on gene isoforms.\nWe have clarified this in the Methods section of our revised manuscript:\n“RNA was purified using poly-T-oligo-attached magnetic beads and the TruSeq\nstranded total RNA kit and TruSeq stranded mRNA kit (Illumina, San Diego,\nUSA) were used to prepare the short-read sequencing libraries.”\n9. CD79A is not an immune regulator, but a B cell marker\nWe thank the reviewer for highlighting this mistake. We obviously agree and have\nchanged the text accordingly:\n“Key DEGs that were increased in advanced plaque vs. early lesion controls\nincluded immune regulators like FLT4 (Schmeisser et al., 2006) or B-cell\nmarkers like CD79A (Ma et al., 2021) and…”\n10. Page 14: It is not clear which datasets are new (unpublished). Clearly distinguish\nbetween new and reanalyzed data.\nWe have clarified the text in our revised manuscript accordingly:\n“Data from our own scRNA-Seq and single cell spatial transcriptomics\nexperiments have been deposited in NCBI's Gene Expression Omnibus. The\npreviously generated datasets, that are already published, are accessible\nthrough GEO Series accession number GSE247238 (scRNA-Seq, patients 1-\n10) (Fidler et al., 2024; Paloschi et al., 2024; Parma et al., 2025; Pourteymour\net al., 2024) or in this manuscript in Suppl. Table 6 and Suppl. Table 7 (bulk\nseq data) (Bonfiglio et al., 2025; Fidler et al., 2024; Paloschi et al., 2024; Tan et\nal., 2025; Wang et al., 2024). Newly acquired datasets for this study are:\nGSE294291 (scRNA-Seq, patients 11-17) and GSE294466 (Xenium spatial\ntranscriptomics).”\n11. Figure 2H: The immune cells are invisible. Either make an extra panel with a\ndifferent scale, or use a log scale.\nWe have adapted Figure 2H to match the style of the previous figures, and split it into\ntwo subplots: Figure 2H and Figure 2I. This allows the reader to put the visualized\nresults to determine differences of ‘Control vs. Plaque’ specimens.\n12. Figure 5A: numbers too small (not legible)\nWe thank the reviewer for this suggestion. To improve readability, we have split the\noriginal Figure 5 into two separate figures. Revised Figure 5 now contains only the\nenlarged cluster map with improved readability, while revised Figure 6 presents the\nscatterplots of the morphological clusters. Figure numbering has been updated\naccordingly throughout the manuscript and figure captions.\n13. Figure 6A: show directionality (which genes were up in smokers?)\nWe thank the reviewer for the suggestion. We updated and revised the volcano plot\nwith a legend clearly indicating directionality of DEGs in smokers.\n14. Page 4: “unpaired (n=232, n=65) and paired (n=42) designs, as well as stable\n(n=97) and unstable (n=126) plaques”. Some may be overlapping. What was the total\nnumber of samples?\nWe revised the text in the Method section accordingly and added the following\nstatement for clarification:\n“After filtering, 14,572 genes were detected across a total of 297 samples\nderived from 252 patients. The dataset comprises transcriptomes from 232\nadvanced plaques and 65 early lesion controls. Among the plaques, 97 were\nclassified as stable and 126 as unstable, with 9 remaining unclassified based\non our histological analysis. For 42 patients, paired samples of plaque and\ncontrol lesions were available, allowing for paired differential expression\nanalysis.”\nFurthermore, we clarified this in the Results section as well:\n“...classified by FC thickness above or below 200µm. The advanced plaques\nwere further classified into stable (n=97) and unstable (n=126) lesions based\non FC thickness above or below 200µm18. Nine advanced plaque specimens\n(of n=232) were rendered non-classifiable, and not included into the stability-\nbased comparison.”\n15. The authors often talk about “carotid plaques”, but much of the resected material\nis not plaque but media. Maybe call them endarterectomy specimens.\nWe understand and partly agree with the reviewer that the word ‘plaque’ might be\nmisleading - at least for some of the specimens being analyzed. However, the media\nis never the main component for any of the CEA specimens being assessed (see HE\nstains shown in Figure 4, Suppl. Fig. 11 and Suppl. Fig. 12 in the manuscript). We\nas well as others have termed these endarterectomy specimens carotid plaques, as\nthe indication for surgical removal is the stenotic plaque and the occurrence of a lipid-\nrich lesion core, which is always present in all our advanced plaque specimens from\nCEA.\nWe added the term “endarterectomy specimens” to the Discussion to acknowledge\nthe reviewer’s point:\n“Understanding plaque the morphology of these endarterectomy specimens is\ncritical to improve clinical decision making for patients at risk for stroke or\nmyocardial infarction”\n16. SMCs, ECs and fibroblasts show up much better in single nuclei sequencing. Many\nstudies used this approach, mitigating the bias. Please mention this and cite some\nstudies.\nWe thank the reviewer for this comment. It is indeed true that single nuclei sequencing\nusually overcomes the immune cell bias getting introduced with single cell RNA-Seq.\nThis is due to the methodology (freezing down and/or fixation of tissue before nuclei\nisolation) as utilized in previous studies (Chou et al., 2022; Song et al., 2024). Also,\nnuclei have similar sizes in different cell types and usually do not lead to clogging of\nthe microfluidics system used for single cell RNA-Seq. Therefore, single nuclei\nsequencing can indeed reveal more realistic proportionality of cell types.\nIn this current study, we however decided to overcome this bias by introducing\nsc-bulk RNA-Seq deconvolution to carotid plaques. In addition to deconvolution, spatial\ntranscriptomics using Xenium is unbiased when it comes to cell type proportions, as\nthe tissues are fixed and all RNA of the plaque in situ equally accessible. Of note,\nsingle nuclei sequencing has the disadvantage of losing all cytoplasmic RNAs, which\nare important regulators of cell function. In our belief, the integration of all three\ntranscriptomic sequencing methods proved to be a valid approach to better determine\nrealistic cell type proportions in human atherosclerotic plaques.\nBased on the reviewer’s comment, we added a sentence in the Discussion\nsection of our revised manuscript:\n“Alternatively, this bias could be overcome by performing single nuclei\nsequencing, with the downside of losing cytoplasmic RNAs (Chou et al., 2022;\nSong et al., 2024).”\n17. The authors write that “VSMC-rich media from our CEA specimens (Fig 1A), a\nstructure that has likely not been analysed in previous studies“. Was there a difference\nin gene expression in the media of adjacent normal vs atherosclerotic artery tissue?\nPaired analysis should be possible.\nWe thank the reviewer for this very relevant suggestion. We indeed had the same\nthought and addressed this issue in a spatial proteomics study that is available as a\npreprint (Sinha et al., 2025) and currently in revision elsewhere. In this study, laser\nmicrodissection of the different subregions (media, NC, FC) was performed, and the\nproteomes of these subregions were analyzed individually. We did not follow this\napproach with spatial transcriptomics, as in particular the necrotic core of advanced\nplaques is rather poorly conserved for RNA expression analysis.\nHowever, we extended our work to include the proposed paired DGE analysis of media\ncells. In the comparison of “Plaque vs. Control”, we identified three DEGs, ITLN1,\nPECAM1, and CD34, all of which were downregulated in media EC_2 cells from\nplaques compared to controls. Importantly, no other cell type exhibited any significant\nchanges in gene expression between “Plaque vs. Control”. The results of this analysis,\nalong with comparisons of additional condition variables within medial cells, are\nprovided in Suppl. Table 15 and Suppl. Table 16. We have also added a description\nof this analysis to the corresponding section in the Methods:\n“To uncomplicate our analysis, we used subsets of data according to our\nscientific questions: for the comparison “Plaque vs. Control” only samples from\nPatients 1-4 were included as they had matching “Control” and “Plaque”\nsamples, all other comparisons were performed across the 12 “Plaque” samples\n(Patient 1-12). For Patients 1-4 with available control and plaque samples, we\nalso performed DGE analysis including their media-derived cells only, due to\nthe media being present and proposed to be relatively similar in both control\nand plaque samples. DGE results can be found in Suppl. Table 2, Suppl. Table\n3 and Suppl. Table 15 and Suppl. Table 16.”\n18. “Four patients provided paired samples”. Were these bilateral endarterectomies,\nwith one side worse than the other? If so, what was the time elapsed between the two\nsurgeries?\nWe thank the reviewer for a chance to clarify our sampling approach. Please also see\nQuestions #5 and #7 by the same reviewer for additional explanation about how we\naddressed changes regarding our sample setup in the text for your clarification.\nCEA-derived plaque specimens are cut into several pieces. The most advanced\nand diseased plaque piece, as well as the healthiest/less diseased piece from the\nsame carotid artery are taken for downstream analysis to capture the dynamic of\nathero-progression.\nWe changed the Methods section accordingly:\n“Four patients provided paired samples, meaning that one CEA specimen from\nthe same carotid artery was divided into early and advanced lesions (Suppl.\nFig. 3), while eight provided advanced plaques only.”\n19. Page 6: ”Correlation analysis revealed a strong inverse relationship between NC\nand lumen” suggests that the sum of the two may be constant and, by inference, the\nfibrous cap area % must also be constant.\nWe appreciate the comment of the reviewer. While we observed a strong inverse\ncorrelation between normalised NC and lumen areas, we didn’t find the sum of the two\nsubregion areas, whether absolute or normalised, to be constant as also shown\npreviously (Schoenhagen et al. 2003; Glagov et al. 1987; Schoenhagen et al. 2000)\n.We have expanded the manuscript with a new figure showing these values plotted\nagainst each other (Suppl. Fig. 12, see below), both of which display considerable\nvariation, and clarified the corresponding section in the Results:\n“Correlation analysis revealed a strong inverse relationship between normalised\nNC and lumen areas (Suppl. Fig. 12B), which could imply a constant sum,\nhowever, this was not observed (Suppl. Fig. 12). In contrast, absolute area\ncorrelations showed no association between NC and FC or intima, suggesting\nNC expansion reduces the relative proportions of other regions (Suppl. Fig.\n12A).”\n20. The definition of controls is unclear. A sample with 45% intima is not a control. It\nseems that the authors call samples without NC controls.\nWe thank the reviewer for a chance to clarify this. Please also see Questions #5,\n#7\nand #18 by the same reviewer for additional explanation about how we addressed\nchanges regarding our sample setup in the text for more clarity.\nIntimal thickening is a natural phenomenon that occurs in almost all populations\n(especially Western) early in life. Control pieces in our study are indeed samples with\nno necrotic core, but showing early signs of intimal thickening and little to none immune\ncell infiltration. We carefully checked our manuscript again to make sure that we are\nnot using the term “healthy” control, as these tissue specimens are clearly not\ncompletely healthy or undiseased arteries (Finn et al., 2010; Masawa et al., 1994;\nPasterkamp et al., 1998; Velican, 1969; Virmani et al., 2000).\n21. It is interesting that there was no difference between symptomatic and\nasymptomatic plaques. This should be emphasized. In many studies, this is the only\ndistinction provided.\nThe reviewer is absolutely correct that most studies solely base their characterization\nand analysis on comparing symptomatic vs. asymptomatic patients. These studies\nmight lack the thorough histomorphological characterization provided in this study.\nFurthermore, studies that identify differences in symptomatic vs. asymptomatic\npatients either focus on much larger datasets (for example 632 bulk RNA-seq samples\n(Mokry et al., 2022)) or solely on differences in immune cell composition (Fernandez\net al., 2019). Spatial transcriptomics using Xenium with limited gene probe sets are\nlikely unsuited to provide enough statistical power to distinguish symptomatic from\nasymptomatic patients. However, one has to also consider in this context that linking\npatient symptoms to plaque phenotypes (stable or unstable) is extremely challenging,\nas the correlation is weak or even completely absent (Libby et al., 2019). This is due\nto many factors being capable of inducing stroke or stroke-like symptoms, including\nparoxysmal atrial fibrillation with embolization, microcirculatory disruption,\nmicrovascular deformations, as well as hemorrhagic insults (Lavallée et al., 2023). In\naddition, the heterogeneous timing in symptomatic patients from clinical event to CEA\nsurgery - and in many cases the sudden initiation of pharmacological treatment\nregimen (e.g., high-dose statins) can rapidly alter plaque morphology and cellular\ncomposition right before lesion removal.\nImportantly, among all spatial transcriptomics studies of human plaques, none\nreported differences in symptomatic vs. asymptomatic patients. All previous studies\neither do not report symptoms at all, or solely use plaques from symptomatic patients.\nFor the reviewer’s convenience, we have summarized all human plaque studies using\nspatial transcriptomics approaches in the table below:\nSun et al, Bleckwehl et Lai et al, Gastanadui Campos et Pauli et al,\nJACC, 2023 al, NCVR, NCVR, 2025 et al, ATVB, al, EMBO under\n2024 2024 Mol. Med., review,\n2025 2025\ntissue Carotid Carotid + Carotid Coronary Coronary Carotid\n[artery] Coronary\n(for spatial only\ncoronary)\nbulk seq n=163 - n=149 - - n=297\n[sample PTLOs only (assumptions\nsize] (no full tissue) from spatial)\nscSeq - n=13 n=9 - - n=17\n[patients] (previous (assumptions\nstudies; no own from spatial)\ndata)\nspatial VISIUM VISIUM StereoSeq GeoMX GeoMX & Xenium\nmethod (FFPE; probe- (FFPE; probe- (fresh frozen CosMX\nbased whole based whole tissue; ~ (~ Xenium)\ntranscriptome) transcriptome) VISIUM HD,\ncapture based)\nspatial n=6 n=12 n=10 n=10 n=9 n=16\nsamples (no (2 controls) (no controls) (no controls) (2 controls) (4 controls)\n[samples; controls)\ntherein\ncontrols]\nsymptomati no report in no report report in no report no report report in\nc/ bulk seq; PTLO bulk spatial\nasymptoma spatial all seq;\ntic status symptomatic spatial all\nsymptomatic\nWe added this finding to our revised Results section:\n“We could not detect any DEGs by comparing advanced plaques of\nsymptomatic (n=67) vs. asymptomatic (n=162) patients (FDR threshold < 0.05;\nminimum FDR in our analysis 0.13).”\n22. Some patients showed “two NCs on each side of the lumen”. This is interesting,\nbut not conclusive. They may be connected in another plane.\nThis is correct, as it cannot be excluded that in another plaque section only one large\nNC is present. In all spatial transcriptomics approaches, only snapshots of plaque\nmorphology are captured and described.\nWe changed the text in our revised manuscript accordingly:\n“They are less structured, showing two NCs on each side of the lumen (\n(potentially one NC in another plane of the plaque; patient 4 and 8) or one large\nsurrounding NC (patient 9)...”\n23. In several scRNA-Seq studies and in at least 2 reviews, about half of the foam cells\nare SMC-derived and half are macrophages. There are two sets of macrophage foam\ncells and at least 3 other macrophage sets that are not foamy. The information obtained\nin the present study should allow these assignments.\nIn human data, it is unfortunately impossible to make conclusive statements on cellular\norigin. All work that confidently confirms SMC origin of macrophages (or macrophage-\nlike cells) stem from mouse lineage tracing studies (Allahverdian et al., 2018; Harman\n& Jørgensen, 2019; Y. Li et al., 2021; Shankman et al., 2015; Wang et al., 2019). Both\nmacrophages and SMCs are highly plastic cells that acquire features of one another\nin a disease context like atherosclerosis. Thus, our ‘human-only’ study does not allow\nfor such interpretation. Currently, others are working on a holistic integration of several\nlineage tracing studies (Sharma et al., 2024). Progress has also been made on\ntranslating (deconvolution and label transfer) murine fate mapping into human\npathology (D. Y. Li et al., 2025).\n24. Some of the correlations with clinical and demographic parameters are not\nstatistically analyzed.\nThank you for your remark. We have plotted the associations between patient\ndescriptors (plaque morphology, patient metadata, and medical conditions) in Suppl.\nFig. 11 and 12; however, not all parameters were statistically analyzed. Peripheral\nartery disease (PAD) was included in the differential gene expression analysis of the\nspatial transcriptomic data (Suppl. Tables 2 and 3) but was erroneously omitted from\nthese supplementary figures. In addition, the captions of the two plots incorrectly\nsuggested that COPD was one of the assessed clinical parameters. A previous version\nof the figure did include COPD, but we removed it since only one patient had COPD,\nrendering statistical analysis unfeasible. We apologize for these errors and have\nupdated both figures to include PAD as an assessed parameter and corrected their\ncaptions accordingly. In addition, we split them into four figures to improve readability\n(see below), and corrected their references in the main text:\nAssociations between plaque morphology, patient metadata, and medical\nconditions showed distinct trends (Suppl. Fig. 11, Suppl. Fig. 12 Suppl. Fig.\n14-17). Patient age inversely correlated with FC absolute area, height correlated\npositively with total plaque area (Suppl. Fig. 11A, 11F Suppl. Fig. 14A, 15C),\nand hypertensive patients presented with increased normalized media areas\n(Suppl. Fig. 12D Suppl. Fig. 17A). Plaque stability and symptoms (e.g., stroke\nor transient ischemic attacks) showed no correlation with subregion\nmorphology.\nAlso, mentioning a single SNP in the smoking section is not appropriate. SNPs were\nnot tested in this study.\nYes, the reviewer is correct: SNPs were not analyzed in this study. We changed the\ntext in our revised manuscript accordingly:\n“... ATM was concordantly upregulated in smokers. This gene has previously\nbeen linked to lung cancer development by its SNP rs189037 (G>A), increasing the\nrisk in non-smokers compared to smokers. …”\n25. Please define the difference between neighborhoods and subregions.\nThank you for your comment and the chance for clarification. We realised that using\nthe word “subregions” in the Methods section describing the neighbourhood cluster\nanalysis was unfortunate and confusing. We apologise for the mistake, and expanded\nthe corresponding section in our Methods to further clarify the goal of neighbourhood\nclustering, as well as the differences between the neighbourhood clusters and\nsubregions:\n“We hypothesized that although the samples show substantial heterogeneity in\ntheir morphology, assessing the local neighbourhoods of cells could reveal\nsimilar subregions spatial domains across samples. Such domains might reflect\nshared biological processes or functions, potentially leading to a better\nunderstanding of plaque formation. Unlike histological subregions, which\nrepresent well-established anatomical compartments of plaques (i.e., lumen,\nintima, media, fibrous cap, necrotic core), we set out to identify in an unbiased\nmanner patterns in local cell neighbourhoods. We applied unsupervised\nclustering of these neighbourhoods to generate clusters representing spatial\ndomains that are characterised by similar cellular composition. Because the\napproach relies on local neighbourhoods and unsupervised clustering, the\nresulting clusters are independent of predefined subregion boundaries.”\nWe clarified the fact that subregions are based on established histological\ncharacteristics in comparison to the otherwise computational heavy manuscript in the\nMethods section:\n“The annotation of sample subregions is a histological analysis, using\nestablished Oxford Plaque Studies18–20 and AHA criteria78–81, and was\nperformed by an experienced cardiovascular pathologist.”\n26. Figure 2E: the choice of colors is unfortunate. Impossible to see the cluster calling.\nWe thank the reviewer for this valuable comment. To improve visibility in Figure 2E,\nwe increased the dot size in the UMAP plots and selected a color palette with higher\ncontrast to better distinguish cell clusters.\nPlease see Question #11 for the improved Figure 2.\n27. Figure 3D: what are “inflammatory mixed cells”? Doublets?\nWe thank the reviewer for this comment and the chance to further clarify this aspect.\nA detailed rationale behind the annotation of “Inflammatory mixed cells” is included in\nthe revised Methods section. To further clarify our hypothesis regarding the nature of\nthese cells and their possible relationship to doublets, we have expanded this section\nof the revised manuscript. For your convenience, we include the revised section with\nthe additional clarification highlighted below:\n“Based on these features we decided to annotate the cell cluster as\n“Inflammatory mixed cells”. We hypothesize that these cells could reflect cell\nfragments, whose cytoplasm was only partly contained in the plane of the\nhistology sample, although they can also be a result of incorrect cell\nsegmentation. Although it is possible that these cells represent doublets, we\nconsider this unlikely. Potential doublets arising from cell segmentation were\nexcluded using Baysor’s “max_cluster_frac” parameter, as recommended by\nthe authors, and the low transcript counts observed in these cells further argue\nagainst the doublet interpretation.”\n28. Figure 3E: Age seems to be a driver for necrotic core. This should be emphasized.\nWe thank the reviewer for this observation. We assessed the relationship between\nnormalised necrotic core area and patient age, which revealed a weak positive trend.\nHowever, as this association was not statistically significant (Suppl. Fig. 12E, please\nsee Question #24 for the figure), we did not include it in the main manuscript.\nAssessment in a larger patient cohort seems necessary to determine whether this trend\nrepresents a genuine biological association.\n29. Figure 5: What statistical test was used to assign the morphological clusters?\nWe thank the reviewer for this comment, and have revised and expanded the Methods\nsection to provide additional details on the clustering approach and its implementation:\n“First, we calculated the cell fraction percentage of the different cell type\nsubstates normalised within the sample, then performed hierarchical clustering\nwith Euclidean distance and Ward’s-linkage method.\nFirst, we calculated the cell fraction percentage of the different cell type\nsubstates normalised within the sample. This enabled us to capture finer\ndifferences in cellular composition across all plaques. Next, we performed\nhierarchical clustering of the samples based on these normalized cell type\nfractions, using Euclidean distance and Ward’s-linkage method. This is an\nunsupervised clustering method and does not rely on a predefined statistical\ntest for determining group assignments.”\nHow do we know there should be 5 clusters? Or is this arbitrary?\nWe also revised the paragraph describing how the final number of clusters was chosen:\n“We agreed on the final number of clusters (4) based on the clustermap and the\nindividual inspection of the samples. …The final number of clusters (4) was\nselected based on cluster robustness (Fig. 5) and further supported by visual\ninspection of the samples (Fig. 6). Based on the dendrogram structure shown\nat the top of the clustermap in Fig. 5, selecting four clusters provided a\nbiologically meaningful and robust separation between groups of samples, while\nat the same time, samples within each cluster displayed consistent patterns in\ntheir spatial morphology (Fig. 6) as described in the “Morphological plaque\nclustering” subsection in Results. These considerations together guided our\nchoice of four clusters.”\n30. Figure 7C: How were patients 4 and 8 selected?\nWe thank the reviewer for the comment and clarified the rationale for our selection of\nsamples included in the main figures. We have added a paragraph to the Methods\nsection:\n“...To showcase the macrophage neighbourhoods, we selected plaques from\npatients 4 and 8 for inclusion in the main figures, as these contained\nmacrophage clusters with the most illustrative spatial patterns. In addition, these\nsamples also contained EC_neovessel neighbourhoods, which are also\nhighlighted in the main text. To illustrate VSMC neighbourhoods, we selected\nsamples with a cell-dense and intact media (plaques from patients 7 and 10,\nand the control from patient 3) for the main figures.”\nReviewer #2 (Remarks to the Author):\nOverall Summary:\nPauli et al. present their investigation seeking to integrate multiple transcriptomic\nmethodologies (bulk RNA seq, scRNAseq, and spatial transcriptomics) as a means to\nmolecularly define carotid atherosclerotic disease. Using arterial and atheroma\nsamples collected during CEA, the investigators attempt to elucidate a spatially-\nresolved molecular architecture underlying plaque development. The authors highlight\nunique plaque morphologies linked to different cell type compositions, disease\nheterogeneity (i.e. symptomatology), and clinical demographics (i.e. smoking, etc) with\ninference towards plaque development and stability. Notably, they uncover spatial\nassociations between HMOX1+ macrophages as they potentially transition to foam\ncells (TREM2+).\nWhile the manuscript is well written and demonstrates technical rigor, the most\nimportant limitation is a lack of methodologic novelty or conceptual impact to our\nmolecular understanding of atherosclerosis. Notably, HMOX1 and TREM2\nmacrophages are well established in the progression of atherosclerosis (PMID:\n39731912 – HMOX1, 38974464 – TREM2) and modulation of macrophage polarity\n(PMID: 30282830 – HMOX1, 38646596 – TREM2). Furthermore, other spatial\ntranscriptomic investigations of arterial atherosclerosis have recently been published\n(PMID: 39234691 and 39715784). Considering these limitations, the current\nmanuscript offers incremental advancement.\nWe appreciate the assessment by the reviewer of our current work. We have tried in\nthe revision to clarify the intention, novelty, and relevance for atherosclerosis and\ncarotid artery disease of our study. We fully agree with the reviewer on his comments\nregarding macrophage biology, their transition and modulation in human\natherosclerosis. We still believe that our current study adds value to the existing\nmacrophage knowledge by spatially locating (with single cell resolution) their subsets\nin early and advanced lesions. Transition from TREM2 high to CCL18 high for example\ncan be specifically located to the highly relevant plaque shoulder region of advanced\nlesions. Unstable carotid artery plaques become prone to rupture in this specific area\nwhere this macrophage transition occurs.\nOur study is not limited to macrophage biology, as it determines endothelial cell\n(EC) heterogeneity (angiogenic and EndoMT-like), spatial distribution of these EC\nsubtypes within the plaque, as well as their highly specific expression patterns in early\nvs. advanced lesions. Angiogenic ECs colocalize with immune cells in areas of the\nnecrotic core, likely providing an entry gate for immune cells that destabilize the plaque.\nSimilar to ECs, various smooth muscle cell (SMC) subclusters surfaced in our study,\nsuggesting specific spatial differentiation patterns between early and advanced carotid\nplaques.\nOur study provides further novelty and disease relevance through the\ncomparison of early vs. advanced lesions in a paired manner (meaning that the plaque\nand control lesion stem from the exact same patient), which has not been analyzed\nwith (multi-)transcriptomics before. We created a new Supplemental Figure 3 to explain\nour approach in sampling and our unique setup of plaque vs. control from the same\npatient better:\nIn addition to being the first study ever to perform Xenium spatial transcriptomics on\nhuman carotid artery plaques, the integration of all three transcriptomic sequencing\nmodalities enables us to provide a comprehensive view on the dynamic of\natherosclerotic lesion progression and destabilization (via comparing early to\nadvanced lesions from the same patient).\nThe table provided below compares all performed studies to date (six in total,\nincluding our own study reviewed here) using varying spatial transcriptomic\ntechnologies. The reviewer mentioned two spatial transcriptomics studies (PMID:\n39234691 and 39715784) that solely analyzed coronary arteries, and not plaques from\ncarotids. Of note, none of the other four publications provide a full evaluation and\nintegration of bulk, single cell and spatial transcriptomics. Our study offers by far the\nlargest bulk RNA-Seq dataset with 297 samples, as well as the richest single-cell RNA-\nSeq dataset (n=17 patients) and the highest number of samples analyzed with spatial\ntranscriptomics (n=16; including 4 controls/early lesions) of all five studies. We are\nfurthermore the only study that can compare plaques from symptomatic vs.\nasymptomatic patients in spatial transcriptomics, as all other studies either do not\nreport symptom status - or only use plaques from symptomatic individuals:\nSun et al, Bleckwehl et Lai et al, Gastanadui Campos et Pauli et al,\nJACC, 2023 al, NCVR, NCVR, 2025 et al, ATVB, al, EMBO under\n2024 2024 Mol. Med., review,\n2025 2025\ntissue Carotid Carotid + Carotid Coronary Coronary Carotid\n[artery] Coronary\n(for spatial only\ncoronary)\nbulk seq n=163 - n=149 - - n=297\n[sample PTLOs only (assumptions\nsize] (no full tissue) from spatial)\nscSeq - n=13 n=9 - - n=17\n[patients] (previous (assumptions\nstudies; no own from spatial)\ndata)\nspatial VISIUM VISIUM StereoSeq GeoMX GeoMX & Xenium\nmethod (FFPE; probe- (FFPE; probe- (fresh frozen CosMX\nbased whole based whole tissue; ~ (~ Xenium)\ntranscriptome) transcriptome) VISIUM HD,\ncapture based)\nspatial n=6 n=12 n=10 n=10 n=9 n=16\nsamples (no (2 controls) (no controls) (no controls) (2 controls) (4 controls)\n[samples; controls)\ntherein\ncontrols]\nsymptomati no report in no report report in no report no report report in\nc/ bulk seq; PTLO bulk spatial\nasymptoma spatial all seq;\ntic status symptomatic spatial all\nsymptomatic\nOther Criticisms include:\n1. Association with symptomatology / stability is uncertain. One area of potential\nsignificance was the assertion that the current findings molecularly define symptomatic\n(vs. asymptomatic) carotid artery stenosis. As presented, however, the data do not\nclearly support these claims. Figure 5 provides an example of this missed opportunity.\nA clearer presentation of the results would have been based upon symptom status as\nopposed to various morphologic clusters.\nWe thank the reviewer for the chance to clarify our approach in carotid artery disease\nand human CEA specimens. We do not primarily characterize plaques based on\npatients’ symptoms, despite symptoms of course being recorded and taken into\naccount as a key clinical parameter. Most studies analyzing plaques solely based on\npatient symptoms have the limitation that they lack histomorphological stratification of\ncarotid plaque biology and stability.\nAs aforementioned in response to Reviewer 1 (comment # 21), correlation of\npatient symptoms with plaque stability is weak and mostly absent in carotid artery\ndisease (Libby et al., 2019). This well-established missing link between patient\nsymptoms and plaque stability is a critical topic in stroke prevention research and\nclinical decision-making (Gupta et al., 2013; Saba et al., 2025; Underhill et al., 2010).\nIn asymptomatic patients, the degree of carotid artery stenosis and plaque burden\nhave traditionally guided treatment (e.g., surgery vs. medical therapy). But it has\nbecome increasingly clear that biological plaque characteristics - not just the severity\nof the narrowing caused by the plaque - play the most pivotal role in determining stroke\nrisk (Libby et al., 2019; Paraskevas & AbuRahma, 2024). The main reason behind the\nweak linkage of symptoms with stability in carotid artery disease relates to the fact that\na magnitude of factors can ultimately cause stroke and TIA (Lavallée et al., 2023). As\naforementioned, additional factors to the here studied unstable carotid plaques include\ncardiac embolization due to accompanying heart diseases (dilative cardiomyopathies,\natrial fibrillation), hemorrhagic events, vascular malformations, and microcirculatory\ndisruption (Lavallée et al., 2023). In addition, the heterogeneous timing in symptomatic\npatients from clinical event to CEA surgery - and in many cases the sudden initiation\nof pharmacological treatment regimen (e.g., high-dose statins) can rapidly alter plaque\nmorphology and cellular composition right before lesion removal potentially skewing\nresults and interpretations thereof.\nThis missing link separates carotid artery disease from coronary artery disease,\nwhere the identification of the culprit lesion within the coronary system causing\nmyocardial ischemia can be easily distinguished (via CT, MRI, angiographical or\nintravascular imaging, electrocardiogram/ECG, and finally histopathology; (Kubo et al.,\n2007)).\nIn response to the reviewer’s comment, we assessed the association between\nplaque stability and patients’ symptoms. We pooled data across all methods for 252\npatients in total (Xenium n=12 plaques, scRNA-seq n=17 plaques, bulk RNA-seq\nn=223 plaques). This was feasible because there were no overlapping patients across\nthe different datasets. We then constructed a contingency table based on the\ncombined data and performed Fisher’s exact test, which indicated no significant\nassociation between plaque stability and symptomatic status (odds ratio = 1.144, p-\nvalue = 0.67). This confirms the above-described, well-established missing link\nbetween symptoms and plaque stability in carotid artery disease.\nDue to this, we and others (Bashore et al., 2024; Dib et al., 2023; Fidler et al.,\n2024; P. Singh et al., 2024; Sun et al., 2023) believe that it is highly relevant to\ncharacterize plaques histomorphologically, and to use plaque biology (e.g., fibrous cap\nthickness) as the main category to understand changes in lesion stability. Future\ndevelopments into AI-driven digital pathology will hopefully provide enhanced tools to\ndiscriminate different plaque phenotypes from histopathology. But these approaches\nare still under development and have yet to be implemented into CEA biobank routines.\nStill, the comparison of early vs. advanced plaque in our study creates an additional\nuseful resource for this previously lesser understood context of human\natheroprogression.\n2. Surprising lack of DEG in bulk RNAseq. Figure 2c shows relatively few differences\nin DEG across stable and unstable plaques. This is in contrast to current literature\n(39234691), yet no explanation is given.\nWe appreciate the comment by the reviewer and the chance to put our results into\ncontext with previous studies. The mentioned study (PMID: 39234691) was based on\na small sample size (5 stable vs. 5 unstable plaques) and was potentially biased by\nselective regional sampling, limiting the generalizability of its findings and increasing\nthe chance of false positive findings. In contrast, our analysis of transcriptomes from\n126 unstable vs. 97 stable plaques provides a more robust assessment, with greater\npower and reduced susceptibility to sample-, technical-, and region-specific biases.\nThe revised manuscript´s Result section has been updated to follow the suggestion of\nthe reviewer to comment on the low count of DEGs when comparing stable and\nunstable plaque:\n“While a previous study (Gastanadui et al., 2024) identified 107 DEGs between\nstable and unstable coronary plaques (n=10 in total), our bulk transcriptomes of\n223 carotid plaques identified only four DEGs, including CXCL5 and PI3K (Fig.\n2C). The lower DEG count in our study reflects lower susceptibility to sample-,\ntechnical-, and regional-specific biases. In addition, plaque stability, when\ndefined solely by FC thickness as a key predictor of clinical outcomes, seemed\nto not be primarily driven by transcriptional differences at the whole-tissue level.\nWe could further not detect any DEGs by comparing plaques from symptomatic\n(n=67) vs. asymptomatic (n=162) patients (FDR threshold < 0.05; minimum\nFDR in our analysis 0.13).”\n3. Underpowered study with significant demographic heterogeneity. It is difficult to\ninterpret some of the results due to relatively few patients (N=4 Ctl) and variable clinical\ndemographics.\nWe understand the concern of the reviewer. It is indeed true that more samples,\nespecially more controls would be needed to draw stronger conclusions. Still, our\ndataset is the largest spatial transcriptomics dataset of human plaques. Apart from the\nstudies by (Bleckwehl et al., 2025) and (Campos et al., 2025), our study is the only one\nthat utilizes controls, and therefore the only study that shows control samples for\ncarotid arteries (as Bleckwehl et al. as well as Campos et al. solely studied coronary\narteries). In comparison to all other spatial transcriptomics studies, we further provide\nthe highest number of additional (and integrated) bulk and single cell RNA-Seq data\nfrom the same CEA patient cohort now getting released with this manuscript.\nFurthermore, we are the first ones to perform a targeted approach of spatial\ntranscriptomics, including long non-coding RNAs, which are highly relevant markers\nfor cell substate detection, while playing significant roles in plaque progression and\nvulnerability.\nFor details in study design and provided data of all current spatial\ntranscriptomics studies, please refer to the Table provided in response to the\nintroductory comment above.\n4. Biased spatial transcriptomic approach based off Panels 1 and 2. Finally, the spatial\ntranscriptomic analysis used a custom (pre-defined) gene set (~500) for cellular and\ndisease status identification. While this approach is helpful in resolving cellular spatial\nlocation, it falls short of providing true spatially-resolved, single cell transcriptomic\nanalysis. This undercuts the novelty.\nIn our current manuscript, spatial transcriptomics was not used as a discovery platform\n(for which the two Xenium panels would indeed be suboptimal given the other\nalternatives). Our aim was to spatially locate and integrate (with the highest available\nresolution) cell types, their substates, and genes that surfaced from our bulk and single\ncell RNA-Seq studies. Thus, the Xenium panels are based on markers from single cell\nRNA-Seq data, enabling us for the first time to locate cells and genes in the same\ntissues (comparing early vs. advanced lesions). Furthermore, Xenium allowed us to\nlocate non-coding RNAs in FFPE tissue currently not included in Visium. As\naforementioned, long non-coding RNAs play a significant role in determining cell\nsubstates based on their level of cell specificity. Visium (without HD; only available\noption at the time we initiated our study in August 2023) loses about 60-70% of the\ntissue area due to the spot anatomy. One spot has the size of 55µm, leading to more\nthan 1 cell (rather 5-10 cells) being sequenced.\nOf note, our manuscript introduces Baysor to achieve single cell spatial\nresolution during the analysis of human atherosclerotic plaques. Compared to Visium\n(HD and no HD), probe-based Xenium provides subcellular resolution (~0.5 µm/pixel),\nwhereas even Visium HD ‘only’ provides a resolution of ~5µm/pixel. For a complex\ndisease like carotid artery disease, the highest possible resolution seems relevant for\ncell type detection and intercellular communication due to the complex and\nheterogeneous nature of these human tissue specimens.\nReviewer #3 (Remarks to the Author):\nSummary\nThis manuscript presents a comprehensive spatially resolved transcriptomic analysis\nof human carotid atherosclerotic plaques. The authors integrate bulk RNA-Seq, single-\ncell RNA-Seq (scRNA-Seq), and single-cell spatial transcriptomics (Xenium platform)\non paired samples from patients undergoing carotid endarterectomy (CEA). Their\nmultimodal approach reveals spatial microenvironments, including distinct VSMC and\nmacrophage subtypes associated with plaque vulnerability, progression, and smoking\nstatus. They propose a model of subluminal HMOX1⁺ macrophage transdifferentiation\ninto lipid-handling TREM2⁺ cells. A novel morphological clustering of plaques into\n“structured” and “chaotic” types is introduced, offering mechanistic insights into spatial\nimmune-VSMC interplay and plaque destabilization.\nMajor Comments\n1. Justification of Technology Choice: The authors employ the Xenium platform for\nspatial transcriptomics but do not justify its selection over other platforms like 10x\nVisium HD. A justification regarding resolution, sensitivity, and gene coverage tradeoffs\nis needed.\nVisium HD was released in March 2024, whereas our Xenium runs started in August\n2023. And Visium (without HD; only available option back then) loses about 60-70% of\nthe tissue area due to the spot anatomy. Also one spot has the size of 55µm, leading\nto 5-10 cells getting simultaneously sequenced. Another reason why we chose Xenium\nover Visium is the fact that Visium for FFPE is a probe based assay that can only\nestablish data for protein coding genes. We also wanted to measure long-non coding\nRNAs, as they are important markers of cellular substates. Visium for fresh tissue (not\nfixed) did not work in our hands despite multiple trials using our plaque specimens due\nto tissue detachment from the slide.\nOur goal was to utilize Xenium, as it reaches subcellular resolution\n(~0.5µm/pixel) compared to Visium’s 55µm spot diameters. Of note, even the new\nVisium HD platform ‘only’ provides a resolution of about ~5µm/pixel, so 10x lower\ncompared to Xenium.\nThe sensitivity (transcript detection/cell) is reported to be higher with Xenium\n(up to 5000 transcripts/cell) compared to Visium and Visium HD. Of note, our\natherosclerotic plaques do not reach this sensitivity due to certain tissue specifics\n(matrix and lipid deposition, calcifications and scar tissue, multiple apoptotic cells and\ndebris thereof).\nGene coverage is obviously higher in Visium compared to Xenium, as Visium\nspans 18,000 genes. This appears ideal for discovery-based approaches. However,\nthe focus of our study was the integration of bulk, single cell and spatial transcriptomics\nfrom the same tissue, which is why we followed a more targeted spatial transcriptomics\napproach using Xenium. Hence, the higher resolution of Xenium became our weapon\nof choice for the more integrative, and less discovery-driven approach.\n2. Immune Cell Overrepresentation Not Corrected: Despite deconvolution, the\noverrepresentation of immune cells in scRNA-Seq is not adequately addressed in the\nspatial analysis. This may bias subregion interpretations.\nDespite the same amount of marker genes per cell type (on average 5 genes/cell type)\nin our Xenium probe set, we were able to confirm the higher account of structural cells\ncompared to immune cells with spatial transcriptomics. Of importance, probes for\nmarker genes for Xenium were selected based on single cell RNA-Seq results (which\nprovided the immune cell overrepresentation). Furthermore, Xenium confirmed single\ncell-bulk RNA-Seq deconvolution (depicted in Fig. 1G and Fig. 2B).\nWe further clarified this aspect in our revised manuscript (Results) and added the\nfollowing sentence:\n“In addition to this deconvolution approach, Xenium-based spatial\ntranscriptomics (designed with an equal amount of probes for all cell types)\nsimilar to bulk RNA-Seq confirmed that the observed immune cell\ndominance(Depuydt et al., 2023; Dib et al., 2023; Fernandez et al., 2019) solely\noccurred in our scRNA-Seq study.”\n3. Smoking Association Requires Caution: The observed association with smoking is\nonly partially supported by bulk and scRNA-Seq data. Acknowledgment of limited\nreproducibility is needed.\nWe agree with the reviewer that more samples in Xenium would be needed to make\nmore cautious statements on the association between smoking and plaque\nmorphology. We therefore tried to validate our findings with bulk- and scRNA-Seq, but\nthis data did not fully support our findings.\nWe further clarify this now and acknowledge the limited reproducibility.\n“However, we did not observe reproducible differences in VSMC cell fractions\nacross smoking status with either bulk- or scRNA-Seq, which reflects rigor in\nassessing reproducibility of our findings (Suppl. Fig. 23A,B).”\n4. Trajectory Analysis of Macrophages: The proposed HMOX1⁺ to TREM2⁺ trajectory\nshould be validated with orthogonal methods beyond TFvelo and visualized more\nexplicitly.\nThank you for your valuable comment. Since TFvelo infers cell trajectories based on\ntranscription factor activity, we considered an orthogonal approach using methods that\nrely on spliced and unspliced mRNA counts. Because our original alignment excluded\nintronic reads, calculation of spliced/unspliced counts was not feasible. We therefore\nrepeated the alignment including intronic reads and subsequently obtained\nspliced/unspliced counts using Velocyto. We chose to infer trajectories with scVelo, an\nestablished splicing count-based method. In addition to pseudotime and velocity\nstreamplots, which display trajectories on a cellular level, we also aimed to visualise\nthe trajectories on the cell type level. For this, we created PAGA-graphs, which show\nvelocity-directed edges between cell types, and plotted pseudotime distribution within\nmacrophage subtypes of interest (Suppl. Fig. 41). As the realignment altered our count\nmatrix (details provided below), we also reran TFvelo on the dataset generated with\nintronic reads included (updated Suppl. Fig. 40). Furthermore, to provide an additional\nlayer of analysis, TFvelo was also applied to our Xenium datasets, as spliced/unspliced\ncounts cannot be derived with this methodology. The results are presented in Suppl.\nFig. 36-39 below, and due to space limitations, in Suppl. PDF 3,4.\nTFvelo analysis of the Xenium datasets revealed a clear ordering of the three\nmacrophage subtypes of interest along pseudotime, with mean pseudotime values\nincreasing in the direction Mac_C1Q → Mac_HMOX1 → Mac_TREM2hi in each\nsample (Suppl. Fig. 36, 37). This pattern was also evident in spatial scatterplots\ncoloured by pseudotime and subtype (Suppl. PDF 3,4). PAGA-graphs showed strong\nedges with the same directionality (Suppl. Fig. 38C, Suppl. Fig. 39C), although the\nlink between Mac_TREM2hi_CCL18hi and Mac_TREM2hi was weak with panel 1.\nWhile not fully concordant across all macrophage subtypes, both TFvelo and scVelo\nconsistently inferred trajectories from HMOX1 towards TREM2hi cells in velocity\nstream plots and PAGA graphs, and both showed lower mean pseudotime values for\nHMOX1+ cells compared to TREM2hi cells (Suppl. Fig. 40; Suppl. Fig. 41). These\nresults could suggest the existence of a transdifferentiation route from HMOX1+ toward\nTREM2hi subtypes, however further experimental proof is necessary.\nWe included the changes in the respective Results section of our revised manuscript:\n“…Based on our spatial and DEG results, we hypothesize a gradual\ntransdifferentiation of HMOX1+ into TREM2hi with Mac_TREM2hi_CCL18 as\nan intermediate, transient state. This hypothesis was corroborated by RNA\nvelocity analysis of our scRNA-Seq data, suggesting a transcriptional trajectory\nover time from C1Q+ subtype towards TREM2hi via HMOX1+ (Suppl. Fig. 31).\nWe assessed this with velocity analysis using TFvelo(J. Li et al., 2024), which\nshowed increasing mean pseudotime along subtypes Mac_C1Q →\nMac_HMOX1 → Mac_TREM2hi_CCL18hi (panel 1) → Mac_TREM2hi (Suppl.\nFig. 36, 37), accompanied by a lumen-to-NC gradient in patients 4, 7, 8, and 10\n(Suppl. PDF 3,4). Partition-based graph abstraction (PAGA) supported the\nsame progression (Suppl. Fig. 38C, 39C), though the link between\nMac_TREM2hi_CCL18hi and Mac_TREM2hi was weak in panel 1. Despite not\nbeing fully concordant across all subtypes, complementary analyses of our\nscRNA-Seq data with two orthogonal methods (TFvelo, scVelo (La Manno et\nal., 2018)) yielded consistent trajectories from HMOX1+ to TREM2hi cells and\nlower pseudotime values for HMOX1+ relative to TREM2hi (Suppl. Fig. 40;\nSuppl. Fig. 41), suggesting a possible transdifferentiation route that requires\nfuture experimental validation. …”\nWe also revised the corresponding section in the Discussion:\n“... In a subset of plaques C1Q+ and HMOX1+ reside close to the lumen, while\nTREM2hi macrophages border and harbour the NC. Interestingly, in a control\nsample, all three substates appear subluminally, without TREM2hi cells\nseparating. RNA velocity analysis revealed a transcriptional trajectory over time\nfrom C1Q+ subtype towards TREM2hi via HMOX1+. Of note, a similar HMOX1+\n→ TREM2hi trajectory, along with HMOX1+ and C1Q+, though occurring in\nopposite direction along with pseudotime was described before(Dib et al.,\n2023). We report a gradual transdifferentiation from HMOX1+ macrophages via\nan intermediate phenotype (Mac_TREM2hi_CCL18) that reduces CCL18\nexpression while increasing SPP1, MMP9, MMP12. Finally, these cells become\nTREM2hi, probably when reaching the crucial shoulder region of the FC. This\nsequence of cells and associated DEGs seem highly orchestrated, based on\nthe repetitive layering and positive correlation of HMOX1+ macrophages with\nTREM2hi. We hypothesise a gradual transdifferentiation from HMOX1+\nmacrophages via an intermediate phenotype (Mac_TREM2hi_CCL18hi) that\nreduces CCL18 expression while increasing SPP1, MMP9, MMP12. Finally,\nthese cells become TREM2hi, probably when reaching the crucial shoulder\nregion of the FC. This sequence of cells and associated DEGs seem highly\norchestrated, based on the repetitive layering and spatial correlation patterns of\nHMOX1+ macrophages with TREM2hi. Xenium RNA velocity analysis\ncorroborated this by revealing a transcriptional trajectory over time from C1Q+\nsubtype towards TREM2hi via HMOX1+, with the HMOX1+ → TREM2hi\ntransition also supported by our scRNA-seq data, as previously reported (Dib et\nal., 2023). Although intriguing, these findings require further experimental\nvalidation. …”\nWe included the following changes in the respective Methods sections:\n“... Alignment of the reads was performed with cellranger’s count function with\ndefault parameters, which included intronic reads.except for setting\n“include_introns” parameter to “false” to exclude intronic reads.“\n“...For the calculation of velocity stream plots, we only included genes with a\nlarger latent time - pseudotime Spearman correlation of 0.4. Pseudotime and\nvelocity streamplots (Suppl. Fig. 40A,B) were created using the scPoli-UMAP\nrepresentation.\nIn addition to the cellular-level trajectory information provided by the\nstreamplots, we also aimed to assess trajectories at the cell type level. To this\nend, we applied the PAGA (partition-based graph abstraction) algorithm (Wolf\net al., 2019) using the inferred TFvelo pseudotime as priors and the macrophage\nlow-level cell types as node groups. The resulting graph, displaying velocity-\ndirected edges between cell types, is shown in Suppl. Fig. 40C. We also plotted\nthe TFvelo pseudotime distribution of the three macrophage subtypes\n(Mac_C1Q, Mac_HMOX1, and Mac_TREM2hi; Suppl. Fig. 40D), which\nshowed ordered spatial distribution patterns in the Xenium data, to better\nvisualize their relative ordering along the TFvelo pseudotime axis.\nTo complement TFvelo with an orthogonal RNA velocity approach, we also\napplied a splicing count-based method to infer cell trajectories. Spliced and\nunspliced expression counts were obtained using Velocyto (La Manno et al.,\n2018) on the .bam files and filtered barcodes generated during alignment with\nCell Ranger’s count function.\nFor the splicing count-based analysis, we used scVelo (Bergen et al., 2020).\nWe first created a subset of our data containing only macrophage cells (1,837\ncells) and calculated first- and second-order moments with\nscvelo.pp.moments() using the scPoli representation, with the number of\nneighbours set to 30. Next, we estimated velocity vectors with scVelo’s\ndynamical model. Inferred pseudotime and velocity stream plots were\ngenerated using the scPoli-UMAP representation (Suppl. Fig. 41A,B). A PAGA\ngraph based on pseudotime inferred by scVelo is shown in Suppl. Fig. 41C,\ntogether with the scVelo pseudotime distribution across the three macrophage\nsubtypes (Suppl. Fig. 41D).”\n“RNA velocity analysis of macrophages in Xenium\nTo further explore the potential relationship between the spatial distribution and\nexpression patterns of macrophage subtypes observed in the Xenium data, we\nassessed their trajectories using RNA velocity analysis. For this purpose, we\napplied TFvelo (J. Li et al., 2024), a transcription factor (TF) based model, since\nspliced and unspliced RNA counts cannot be derived from Xenium due to its in\nsitu hybridization methodology.\nWe subsetted our dataset to include only macrophages (panel 1: 35,265 cells;\npanel 2: 28,747) and performed preprocessing with the parameters\nrecommended by the authors, using the batch-corrected scVI-UMAP\nrepresentation as the low-dimensional embedding and acosh-normalised\ncounts. We chose acosh normalisation because some cells had erroneously\nsmall area estimates from Baysor, which led to extreme values in area-\nnormalized counts. These outliers interfered with TFvelo and prevented reliable\ntrajectory inference. For the calculation of velocity stream plots, we only\nincluded genes with a larger latent time-pseudotime Spearman correlation of\n0.3. Pseudotime and velocity streamplots are shown in Suppl. Fig. 38A,B for\npanel 1, and in Suppl. Fig. 39A,B for panel 2. Similar to our trajectory analysis\nof our scRNA-seq data, we constructed PAGA graphs for both panels using\nTFvelo-inferred pseudotime as priors and macrophage low-level cell types as\nnode groups (panel 1: Suppl. Fig. 38C; panel 2: Suppl. Fig. 39C). For panel\n1, we included the cell substate Mac_TREM2hi_CCL18hi into the construction\nof the PAGA graph, as well as in the colouring of cells in the streamplot.\nLeveraging available cell coordinates, we plotted spatial scatterplots of the\nindividual samples with macrophages coloured by TFvelo pseudotime and cell\ntype respectively (Suppl. PDF 3,4).\nOf note, macrophage cells with panel 2 dataset expressed more TFs (total\nnumber of TFs: 12; mean number of TFs/gene: 5.3), than panel 1 (total number\nof TFs: 5; mean number of TFs/gene: 1).\nFor your convenience, we included the updated Suppl. Fig. 40 and newly created\nSuppl. Fig. 41, Suppl. Fig. 38, Suppl. Fig. 39, Suppl. Fig. 36, Suppl. Fig. 37 as\nshown below:\nIncluding intronic reads in the new alignment required rerunning our complete scRNA-\nseq pipeline. This resulted in more cells (without introns: 12,149 cells; with introns:\n14,641 cells) and more genes (without introns: 23,362 non-zero count genes; with\nintrons: 28,816 non-zero count genes), however, the overall cell type landscape did\nnot change significantly:\nCell fractions per sample also didn’t show significant changes. We updated the\ncorresponding section in the Results, as well as subplot Fig. 2F (see original, without\nand updated, with introns version below):\n“...Microfluidics-based capture biases in combination with enzymatic and\nmechanical digestion is a potential issue with scRNA-Seq, leading to an\noverrepresentation of immune cells. Indeed, our scRNA-Seq revealed high\nmean macrophage (27%24%) and T cell fractions (23%29%), together\nexceeding VSMCs (24%22%) and ECs (15%; Fig. 2F)....”\nIncluding intronic reads, however, resulted in a greater number of differentially\nexpressed genes across cell types. Consequently, we split the original Suppl. Fig. 1\ninto two figures (updated Suppl. Fig. 1 and Suppl. Fig. 2, see below), and revised the\nResults section to reflect these changes:\n“Comparative analysis of scRNA-Seq data identified cell type-specific\ndifferentially expressed genes undetected by bulkRNA-Seq in plaque vs.\nadjacent control tissues (Suppl. Fig. 1 G-H A-D), including ALOX5 IFI27 and\nHIF1A CTSB in IL10high macrophages fibroblast-like VSMCs. Interestingly,\nstable vs. unstable lesion comparisons revealed an even greater number of\nDEGs (Suppl. Fig. 1A -F X4 A-K), such as CXCR4 and HSPA5 in proliferating\nT cells, both markers of T cell infiltration in atherosclerosis(Dong et al., 2022;\nKumar et al., 2006). Additionally, CXCL2 was enriched in ECs, while CCL18\nCXCL3 was altered in C1Q+ macrophages, both linked to atherosclerosis (A.\nSingh et al., 2024; Yan et al., 2021).”\n5. Subregion Definition and Reproducibility: Manual annotation of regions may\nintroduce bias. Providing inter-rater reliability or uncertainty metrics would enhance\nrobustness.\nCurrently, there is no automated subregion annotation tool for human plaques\navailable. Manual subregion annotation is a well established method in the field (J. N.\nRedgrave et al., 2008; Stary, 2000; Stary et al., 1992, 1994, 1995). In our Xenium\napproach, as we only had DAPI images available (histological stainings were not\nembedded in the XOA at that time) and had to manually draw the subregions based\non a cardiovascular pathology expert’s annotation.\nWe edited our Methods section for clarification:\n“The borders of the annotated subregions were initially determined on the HE\nstained slides by a pathology expert, then transferred to the DAPI-stained slides\nmanually by manually drawing the determined subregion borders by hand in\nXenium Explorer.”\nDigital pathology, including automated subregion annotations, is the tool of the future,\nbut unfortunately not available yet for this study.\n6. Sample Size for Spatial Analysis: The n=12 for spatial transcriptomics limits\ngeneralizability. Interpretations of clustering and risk associations should be more\ncautiously framed.\nWe are very well aware that our sample size for spatial transcriptomics is limited but\nwant to stress once more that the overall concept of our study was not to use spatial\ntranscriptomics for discovery, but for integration and deconvolution with our in parallel\nperformed single cell- and bulkRNASeq study. In response to the reviewer’s\nsuggestion, we carefully assessed the text again for cautious framing and believe that\nwe address this appropriately in the resubmitted manuscript. As aforementioned to\nreviewers 1 & 2, compared with the other existing spatial datasets on human plaques,\nour effort provides the largest dataset to date:\nSun et al, Bleckwehl et Lai et al, Gastanadui Campos et Pauli et al,\nJACC, 2023 al, NCVR, NCVR, 2025 et al, ATVB, al, EMBO under\n2024 2024 Mol. Med., review,\n2025 2025\ntissue Carotid Carotid + Carotid Coronary Coronary Carotid\n[artery] Coronary\n(for spatial only\ncoronary)\nbulk seq n=163 - n=149 - - n=297\n[sample PTLOs only (assumptions\nsize] (no full tissue) from spatial)\nscSeq - n=13 n=9 - - n=17\n[patients] (previous (assumptions\nstudies; no own from spatial)\ndata)\nspatial VISIUM VISIUM StereoSeq GeoMX GeoMX & Xenium\nmethod (FFPE; probe- (FFPE; probe- (fresh frozen CosMX\nbased whole based whole tissue; ~ (~ Xenium)\ntranscriptome) transcriptome) VISIUM HD,\ncapture based)\nspatial n=6 n=12 n=10 n=10 n=9 n=16\nsamples (no (2 controls) (no controls) (no controls) (2 controls) (4 controls)\n[samples; controls)\ntherein\ncontrols]\nsymptomati no report in no report report in no report no report report in\nc/ bulk seq; PTLO bulk spatial\nasymptoma spatial all seq;\ntic status symptomatic spatial all\nsymptomatic\n7. Neighbourhood Cluster Analysis: Statistical support (e.g., permutation tests, spatial\nproximity significance) for proposed functional zones is missing and should be\nincluded.\nWe thank the reviewer for this valuable comment. To provide statistical support, we\nhave performed permutation tests by shuffling cell substate labels within samples\n10,000 times, to create a background distribution of cell neighbourhoods. We then\ncalculated neighbourhood enrichment score, a test metric defined as the total number\nof unique cell-cell adjacency edges between substate pairs within a cluster, using the\ngenerated background local neighbourhoods. Enrichment scores were calculated\nwithin cells of the original clusters, and the significance of the observed clusters was\nassessed by comparing the background distributions to the values obtained from the\nobserved neighbourhoods.\nPermutation tests of neighbourhood enrichment score detected significantly enriched\nor depleted cell substates in all clusters, both for panel 1 and 2. Due to space\nlimitations, the complete results of the enrichment permutation tests are provided in\nSuppl. Tables K1 and K2, while enrichment heatmaps showing significantly enriched\nor depleted substates are included in the newly created Suppl. PDFs D1 and D2. We\nhave expanded the Methods section accordingly (see revised text below):\n“...In the final step we performed Leiden-clustering of the cell neighbourhood\nmatrices using scanpy’s pp.neighbors() and tl.leiden() functions, resulting in cell\nneighbourhood clusters (Suppl. Fig. 60, Suppl. Fig. 61).\nTo provide statistical support for the clusters, we performed permutation testing\nusing neighbourhood enrichment score, a metric defined as the total number of\nunique cell-cell adjacency edges between substate pairs within a cluster. We\nhypothesised that if our clusters indeed capture distinct spatial patterns, by\ncalculating the neighbourhood enrichment score, we could identify which cell\nsubstate pairs are significantly enriched in each other’s neighbourhoods within\na cluster.\nFirst, we created a background distribution of cell neighbourhood matrices by\nrandomly shuffling cell substate labels 10,000 times. Because samples\nexhibited large differences in cell substate composition, we performed label\nshuffling within samples. This was necessary since shuffling across samples\nwould have violated the exchangeability assumption of permutation testing\n(Good, 2002) and could have led to inflated p-values.\nNext, for neighbourhood enrichment analysis, we constructed spatial\nneighbourhood graphs by connecting cells within a 100 µm radius, with each\nedge representing a pair of cells within this distance. For each neighbourhood\ncluster, we restricted the graph to the cells assigned to that cluster and\ncalculated neighbourhood enrichment as the total number of unique undirected\nadjacency edges connecting substate pairs within the cluster-specific graph.\nThis was again computed for both the observed and randomly shuffled\nsubstates, resulting in one enrichment score matrix per cluster for the observed\ncase and 10,000 matrices per cluster for the permutation case.\nEmpirical left and right p-values were calculated using the Monte-Carlo p-value\nmethod as described by Phipson et al (Phipson & Smyth, 2016):\n𝑙+1 𝑟+1\n𝑝 = ; 𝑝 = ;\n𝑙𝑒𝑓𝑡 𝑟𝑖𝑔ℎ𝑡\n𝑁+1 𝑁+1\nwhere p and p are the left and right empirical p-values, l is the number of\n𝑙𝑒𝑓𝑡 𝑟𝑖𝑔ℎ𝑡\ntimes enrichment score from the permutation was lower than the enrichment\nscore of the observed case, r is the number of times enrichment score from the\npermutation exceeded the enrichment score of the observed case, and N is the\nnumber of permutations. In addition, Agresti-Coull confidence intervals (CI)\n(Bozhko et al., 2021) of the empirical left and right p-values were also\ncalculated. Two-sided empirical p-values and their Agresti-Coull confidence\nintervals were calculated using the formula applied in\nscipy.stats.permutation_test() function (Scientific Python Forum, n.d.):\n𝑝 :𝑚𝑖𝑛(2×𝑝 ,2×𝑝 );\n𝑡𝑤𝑜−𝑠𝑖𝑑𝑒𝑑 𝑙𝑒𝑓𝑡 𝑟𝑖𝑔ℎ𝑡\n𝐶𝐼 :𝑚𝑖𝑛(2×𝐶𝐼 ,2×𝐶𝐼 ).\n𝑡𝑤𝑜−𝑠𝑖𝑑𝑒𝑑 𝑙𝑒𝑓𝑡 𝑟𝑖𝑔ℎ𝑡\nRaw empirical two-sided p-values were then corrected for multiple testing by\napplying the Benjamini-Hochberg correction method, and reported as false\ndiscovery rate (FDR).\nPermutation test results for both panel 1 (Suppl. PDF 6, Suppl. Table 17) and\npanel 2 (Suppl. PDF 7, Suppl. Table 17) neighbourhood clusters identified\nsignificantly enriched or depleted substate pairs in every cluster compared to\nthe random background. Based on these results we believe our neighbourhood\nclustering approach was able to identify spatial domains with significantly\ndifferent compositional patterns, even across multiple samples.\nHigh resolution spatial scatterplots of all samples with cells coloured by low-\nlevel cell types and neighbourhood clusters are contained in Suppl. PDF 8 for\npanel 1 genes, and in Suppl. PDF 9 for panel 2 genes.”\n8. Clinical Utility and Predictive Value: The translational relevance of findings should\nbe discussed. Can DEGs or spatial patterns serve as prognostic markers for plaque\nrupture?\nThis is an interesting thought and question by the reviewer. We indeed believe that the\ncombination of multiple transcriptomic profiling methods will help pave the way for\nnovel therapies as well as biomarkers of carotid plaque stability. What is currently\nlacking is data from longitudinal carotid disease and stroke studies, in which these\npotential markers could be tested for their prognostic value. Integration of additional\ndata information from histology (using digital pathology tools as mentioned multiple\ntimes in this point-by-point response) on plaque stability and even more promising\ncorrelation to imaging (in particular CTs from carotid arteries) could further enhance\ndecision making in patients with asymptomatic carotid artery disease. Our approach of\nassessing plaque instability with integrated multi-transcriptomics enhances our\nunderstanding of advanced plaque biology. Novel therapeutic targets can be linked,\nbased on our data to cell type and spatial location (and thus relevance), within\nadvanced plaques. Here, our dataset provides the most complete transcriptomic view\non early vs. advanced plaques to date.\nWe addressed this thought in our revised Discussion section:\n“Future studies on the induction and therapeutic targeting of cellular plasticity\nmay yield novel markers and interventions for plaque vulnerability, together with\ndigital pathology approaches and enhanced vascular imaging modalities.”\nMinor Comments\n1. Terminology Clarification: Terms like 'structured' vs. 'chaotic' could be more\nquantitatively or mechanistically defined.\nIn this work, we tried to use simplified and understandable terms. We discussed this\ninternally quite a lot, as it is a difficult task and challenging to choose additional or other\nterms to enhance the readers’ understanding of plaque morphology. We introduced\nthese terms for this study to highlight the complexity of human plaques. We want to\nshow that easy discriminations like “stable and unstable” or \"symptomatic and\nasymptomatic” (terms that are used in all other studies that deal with plaques in a\nspatial context (Bleckwehl et al., 2025; Gastanadui et al., 2024; Lai et al., 2025; Sun\net al., 2023) are not sufficient to capture their complexity. Fig. 5A, shows the number\nof cells in each cluster which serves as a quantitative definition, furthermore we\nspeculate about mechanisms based on spatial location.\nBased on the reviewer’s comment we adapted the Results section for clarification:\n“For intelligibility, we decided to term the two main groups of our plaque clusters\n“structured” (morphological clusters 1 & 2) or and “chaotic” (morphological\nclusters 3 & 4).”\n2. Figure Readability: Dense color legends in Figures 6–8 make interpretation difficult.\nConsider enhancing clarity.\nWe appreciate this comment. At the initial submission, the figures were uploaded at 96\ndpi due to a misunderstanding, which reduced their readability. To improve clarity, we\nhave increased the dot size in the spatial scatterplots of Figures 6–8 (renumbered as\nFigures 7–9) and re-uploaded them at 300 dpi.\n3. Figure 2 caption: “-log10 value” => “-log10 p-value”\nWe corrected the figure caption accordingly.\nReviewer #3 (Remarks on code availability):\nNo code is provided.\nWe apologize for the delay in providing the code. Setting up the Docker environment\ntook additional time due to conflicting package versions. In the meantime, we prepared\ncode to reproduce the main figures, as requested by the editors, since running the full\nanalysis pipeline on raw data is not feasible due to storage constraints and time\nlimitations. All code will be made publicly available upon publication."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors addressed all concerns. However, figure 2F-I is still unsatisfactory. Please show the immune cells using a different y axis scale. Since there are few immune cells, differences cannot be appreciated from the current figure", "reviewer_2": "Pauli et al. represent their investigation seeking to integrate multiple transcriptomic methodologies (bulk RNA seq, scRNAseq, and spatial transcriptomics) as a means to molecularly define carotid atherosclerotic disease. Using arterial and atheroma samples collected during CEA, the investigators attempt to elucidate a spatially-resolved molecular architecture underlying plaque development. The authors highlight unique plaque morphologies linked to different cell type compositions, disease heterogeneity (i.e. symptomatology), and clinical demographics (i.e. smoking, etc) with inference towards plaque development and stability.\nPrevious critiques raised questions regarding 1) associations with disease symptomatology, 2) Comparison of the presented results to previously published data (coronary and carotid), 3) Novelty and clinical utility of the results, and 4) clarification of methodologies and sample size. To this end, the investigators have provided a robust response that is highly responsive to the reviewer concerns. This includes adjustments to the data presentation, new analysis, and additional commentary throughout the manuscript.\nRegarding novelty, the team has provided sound arguments highlighting the relative impact their results may have on the molecular understanding of carotid artery disease and plaque progression. While the HMOX1 and TREM2 phenotypes are well established, the manuscript provides other novel findings including 1) spatial resolution of macrophages at the unstable plaque shoulder, 2) EC co-localization for neovascularization and immune cell entry into the necrotic core / plaque, 3) SMC transitions associated with plaque heterogeneity. While the clinical utility and implementation of many of these results remains to be determined, the work nonetheless provides a thorough examination of carotid atherosclerotic disease.\nOther Criticisms that have been adequately addressed:\n1. Association with symptomatology is uncertain. Concerns have been well addressed and edited within the manuscript.\n2. Surprising lack of DEG in bulk RNAseq. Concerns have been well addressed and edited within the manuscript.\n3. Underpowered study with significant demographic heterogeneity. Concerns have been well addressed and edited within the manuscript.\n4. Biased spatial transcriptomic approach based off Panels 1 and 2. Concerns have been well addressed and edited within the manuscript.", "reviewer_3": "All my comments are adequately addressed."}, "author_response": "Reviewer #1 (Remarks to the Author):\nThe authors addressed all concerns. However, figure 2F-I is still unsatisfactory. Please\nshow the immune cells using a different y axis scale. Since there are few immune cells,\ndifferences cannot be appreciated from the current figure\nAnswer: We appreciate the positive comments by this reviewer. We have made\nadaptations to Figures 2F-I in our revised manuscript to better address and\nvisulize the low amount of immune cells by choosing a varying scale for the y-\naxes.\nReviewer #2 (Remarks to the Author):\nPauli et al. represent their investigation seeking to integrate multiple transcriptomic\nmethodologies (bulk RNA seq, scRNAseq, and spatial transcriptomics) as a means to\nmolecularly define carotid atherosclerotic disease. Using arterial and atheroma samples\ncollected during CEA, the investigators attempt to elucidate a spatially-resolved\nmolecular architecture underlying plaque development. The authors highlight unique\nplaque morphologies linked to different cell type compositions, disease heterogeneity\n(i.e. symptomatology), and clinical demographics (i.e. smoking, etc) with inference\ntowards plaque development and stability.\nPrevious critiques raised questions regarding 1) associations with disease\nsymptomatology, 2) Comparison of the presented results to previously published data\n(coronary and carotid), 3) Novelty and clinical utility of the results, and 4) clarification of\nmethodologies and sample size. To this end, the investigators have provided a robust\nresponse that is highly responsive to the reviewer concerns. This includes adjustments\nto the data presentation, new analysis, and additional commentary throughout the\nmanuscript.\nRegarding novelty, the team has provided sound arguments highlighting the relative\nimpact their results may have on the molecular understanding of carotid artery disease\nand plaque progression. While the HMOX1 and TREM2 phenotypes are well\nestablished, the manuscript provides other novel findings including 1) spatial resolution\nof macrophages at the unstable plaque shoulder, 2) EC co-localization for\nneovascularization and immune cell entry into the necrotic core / plaque, 3) SMC\ntransitions associated with plaque heterogeneity. While the clinical utility and\nimplementation of many of these results remains to be determined, the work\nnonetheless provides a thorough examination of carotid atherosclerotic disease.\nOther Criticisms that have been adequately addressed:\n1. Association with symptomatology is uncertain. Concerns have been well addressed\nand edited within the manuscript.\n2. Surprising lack of DEG in bulk RNAseq. Concerns have been well addressed and\nedited within the manuscript.\n3. Underpowered study with significant demographic heterogeneity. Concerns have been\nwell addressed and edited within the manuscript.\n4. Biased spatial transcriptomic approach based off Panels 1 and 2. Concerns have\nbeen well addressed and edited within the manuscript.\nAnswer: We appreciate that the reviewer agrees with us on the clarifications\nmade during the revision process of our manuscript on novelty, clinical utility and\nselection of sample characterization utilized in our study.\nReviewer #3 (Remarks to the Author):\nAll my comments are adequately addressed.\nReviewer #3 (Remarks on code availability):\nThe code provided in code ocean are well organized. But the code ocean platform has\nsome issues running the code with the 'Reproducible Run'. Also, only R code shows up\nproperly. The ipython notebook must be downloaded and run locally, which is beyond my\nbandwidth.\nAnswer: We thank the reviewer for acknowledging that we have successfully\naddressed all of her/his previous concerns. We believe to have cleared all\nconfusion regarding code availability in Code Ocean."}]

Decision letter