Vascular smooth muscle cell state trajectories mediate molecular mechanisms of coronary disease risk

Domain: biology
Published: 2026-03-17
Editorial outcome: major_revision

Abstract
Vascular smooth muscle cells contribute to heritable coronary artery disease risk and undergo complex transitions to multiple disease-related phenotypes. To investigate the genetic basis of these trajectories, we develop a dense timecourse single-cell transcriptomic and epigenetic map of atherosclerosis in a murine disease model accompanied by high-plex in situ spatial data. Using temporal data and probabilistic fate modeling, we identify key transcription factors that drive cell state changes through a combination of network-based prioritization and in silico transcription factor perturbation. Parallel knockout studies of validated coronary artery disease gene Tcf21 uncover its molecular mechanisms in smooth muscle cell transition, due in part to a role regulating the transition of smooth muscle cells in the secondary heart field. Integrating the murine atlas with human coronary artery disease genetics pinpoint smooth muscle cell phenotypes that mediate disease risk, highlighting causal disease mechanisms. Together, these studies resolve atherosclerosis trajectories at single-cell resolution and identify genetic causal transcriptomic and epigenomic mechanisms of coronary artery disease risk.

Introduction
Cardiovascular diseases, principally coronary artery disease (CAD) and stroke, are the worldwide leading cause of global mortality1. Therapies that modify classical environmental and metabolic factors have ameliorated a portion of CAD risk2, but more than half of the risk can be attributed to common inherited genetic variation by affecting vessel wall pathways that mediate disease pathophysiology. These genetic factors remain unidentified and untreated3-6. While extensive studies have investigated the cellular and molecular features of atherosclerosis, they have been unable to establish causality and thus hinder the translation towards vascular wall directed therapies7.

The smooth muscle cell (SMC) lineage, which appears to make the largest contribution to disease risk8-10, undergoes extensive complex phenotypic transitions that have not been well characterized at the cellular and molecular level, making it difficult to assign causality to specific gene programs11-16. The few SMC GWAS genes studied in detail suggest that complex transcriptional programs can direct cellular trajectories that lead to fibroblast-like (fibromyocyte, FMC) or osteochondrogenic (chondromyocyte, CMC) phenotypes, with these cellular phenotypes being proposed to mediate opposing effects on disease risk12,13,17. However, this paradigm is in conflict with genomic analyses of limited single-cell data suggesting that all FMC transition to CMC over the course of atherosclerotic lesion development9,12,18. Molecular and cellular studies of SMC transitions to human genetic data linking locus associations to individual SMC gene causality and direction of effect, are needed to promote progress in the field.

Studies reported here are aimed at the comprehensive characterization of genes and gene programs that function in SMC to mediate phenotypic transitions and disease risk. We describe a mouse transcriptomic and epigenomic atlas of atherosclerosis with single-cell RNA and chromatin accessibility (ATAC) data collected over a dense timecourse in a well-accepted mouse model, and further correlate SMC cellular molecular phenotypes with their corresponding spatial niche using high throughput in situ RNA hybridization (Xenium). We map the SMC lineage cellular trajectories with real-time course gene expression, chromatin accessibility data and advanced trajectory inference methods such as the Waddington Optimal Transport algorithm, identifying the genes and gene regulatory networks that mediate the transitions to FMC and CMC cellular phenotypes. Further, we investigate how these trajectories are affected by the CAD protective Tcf21 gene, mapping the regulatory network downstream of this gene, further identifying collaborating transcription factors (TFs) that mediate genome-wide regulation of SMC phenotype transitions and disease risk.

Results
Single-cell atherosclerosis timecourse atlas construction
We performed single-cell RNA sequencing (scRNAseq) and single-nucleus assay of transposase-accessible sequencing (scATACseq) with aortic root tissue from ApoE-/- (ApoE KO) mice to characterize the complex genetic regulatory networks that control the developmental cascade of smooth muscle cell (SMC) phenotypic transitions in the context of atherosclerotic stress13. These mice, also expressing a tamoxifen-inducible Cre recombinase driven by the SMC-specific Myh11 promoter and a Cre-activatable tdTomato reporter gene, were fed a high fat diet (HFD) across 7 timepoints for scRNAseq and 6 timepoints for scATACseq assays (Fig. 1A). Aortic root tissues were collected, digested, and subjected to droplet-based cell capture, and independent RNA and transposed DNA sequencing13. The resulting data were subjected to dimensionality reduction, clustering and visualization with Seurat, providing individual cell clusters that we and others have identified previously (Fig. 1B)11-16. For this study, we subsetted lineage traced SMC clusters, transferred scRNASeq labels onto scATAC cells and co-embedded RNA and ATAC modalities, resulting in a total of 96,195 high quality cells across both modalities (Figs. 1C, D, Supplementary Figure. 1A-K, Supplementary Data 1, Methods)13,19.

Defining molecular and spatial patterns in the SMC lineage
We expanded upon previously identified SMC cell states11-16,19 and identified six distinct SMC lineage phenotypes supported by both RNA and chromatin accessibility patterns (Fig. 1C, D, Supplementary Figure. 1L, Supplementary Data 2). We observed high level expression of canonical SMC differentiation markers (Cnn1, Tagln, Acta2, Myh11) in three clusters (SMC-1, 2, 3). While SMC-1 primarily expressed these markers, SMC-2 demonstrated expression of markers such as Igfbp2, indicating an early phenotypic transition state20,21. A distinct SMC-3 population additionally enriched for myocardial markers Tnnt2 and Nkx2-5, identified the previously characterized population of SMC at the base of the aortic root arising from the secondary heart field(SHF)21,22. All of the phenotype clusters represented cells of the SMC lineage, as indicated by lineage tracing with the Myh11-Cre transgene (Supplementary Figure. 1H). It is important to note that recombination efficiency is quite high, ~95%, so cells expressing even low levels of Myh11 will activate the high-level expression of the tdTomato reporter, which does not really reflect endogenous Myh11 expression level that might reflect some heterogeneity in the SMC cell population.

All FMC were characterized by robust expression of fibroblast and epithelial mesenchymal transition (EMT) markers such as Vcam1, Lum, Bmp1 and Pdgfrb (Supplementary Figure. 1L). However, we identified two FMC sub-populations demonstrating discrete patterns of gene expression (Fig. 1E, F), chromatin accessibility (Fig. 1C, D, H, I), functional enrichment (Fig. 1K, L), and spatial localization (Fig. 1O-T, Supplementary Figure. 2A-G). Taken together, these data demonstrate the integrity of each FMC molecular subtype. One group, termed FMC-1, expressed specific markers including Tcf21, Fbln1, and targets of interferon gamma signaling including H2-Aa/Ab1/Eb1, Cd74 and Il33 (Supplementary Figure. 1L). Validation with in situ RNA hybridization of cluster-specific FMC-1 markers demonstrated expression largely encompassing the media and sparsely at the fibrous cap (Fig. 1R, Supplementary Figure. 2E, Supplementary Figure. 3A-C). Transcriptional functional enrichment with gene set enrichment analysis (GSEA) using GO biological process gene sets for FMC-1 cells identified processes such as immune activation and cytokine production (Fig. 1K, Supplementary Figure. 1L, Supplementary Data. 3). The second cluster of cells, termed FMC-2, were identified with specific markers Thbs1 and Notch3 localized to the fibrous cap (Supplementary Figure. 3 D, E) or spanning the intimal plaque (Col8a1, Fbln2, Ttc9) (Fig. 1S, Supplementary Figure. 2F, Supplementary Figure. 3F, G, H). Although FMC-2 shared a number of terms with FMC-1, expressed genes were uniquely associated with collagen fibril organization, regulation of lipid transport and wound healing (Fig. 1L, Supplementary Data. 3). Interestingly, a number of genes (Loxl1, Bmp1, Vcam1) were found to have vascular wall expression patterns encompassing those seen for both FMC-1 and FMC-2, possibly reflecting cells in transition (Supplementary Figure. 3I, J, K).

By comparison, the CMC cluster identified a distinct transition with a highly restricted chromatin accessibility pattern, and restricted expression of marker genes Col2a1 and Ibsp localized to the base of the plaque (Fig. 1G, J, Q, Supplementary Figure 1L, Supplementary Figure. 2G, Supplementary Figure. 3L, M) as we have shown previously13,17,19. CMC pathways were consistent with osteochondrogenic processes, as described by several groups11-16,19 (Fig. 1M, Supplementary Data. 3). To translate these spatial findings from mice to human, we further performed label transfer onto publicly available scRNA data from human donor transplant coronary arteries and observed a similar distribution of FMC-1, FMC-2 and CMC by single-cell features and spatial transcriptomics23 (Supplementary Figure. 2H-M).

SMC fate trajectories identify cell transition gene programs
We utilized force directed layouts to visualize transcriptional relationships, in real-time, across the SMC atherosclerosis trajectories (Fig. 2A). We first observed the appearance of FMC-1 cells at 5 weeks of HFD with a subsequent increase in FMC-2. This observation is consistent with a migratory path for SMC lineage cells from the media to the fibrous cap and subsequently down into the intimal plaque24-28 (Fig. 1O, R, S, Fig. 2A, B, C, Supplementary Figure. 3A-H). Consistent with this possibility, FMC-2 were also enriched for osteoblast progenitor gene signatures (Fig. 1M, Supplementary Figure. 1L, Supplementary Figure. 4A). Osteoblast progenitors have been characterized in bone development, and this module score29 determined with scRNAseq data characterizing these cells on their path toward osteochondroprogenitors. CMC abundance dramatically increased at 12 weeks of high fat diet, correlating with spatial localization at the plaque base (Supplementary Figure. 3L, M). In addition, transcriptomic patterns suggested increasing pathologic signatures including senescence, EMT, angiogenesis, apoptosis, and efferocytosis (Supplementary Figure. 4B-F, Supplementary Data. 4) as the cells transitioned towards the CMC cell state at the plaque base, which is consistent with the acellularity in this region13,25.

To examine the translocation of FMC from media to cap to plaque, we performed additional studies to map temporal gene expression of FMC-1 and FMC-2 marker genes Fbln1 and Fbln2 across serial sections to highlight their distinct expression patterns. RNAScope performed at baseline showed that FMC-1 marker Fbln1 expression was seen in the media while expression of FMC-2 maker Fbln2 was not (Supplementary Figure. 3N). Then, at 7 weeks of high fat diet (HFD), there was enrichment of Fbln1 in the media and at lower levels at the cap, while Fbln2 was more intimally restricted with lower medial expression compared with FMC-1. (Supplementary Figure. 3O). These findings aligned with our single-cell data which showed low level FMC-1 but no FMC-2 marker expression at baseline while FMC-1 and FMC-2 expression increased after week 7 of HFD (Fig. 2A, B, C).

We further investigated the lesion cellular anatomy at 16 weeks HFD with high quality Xenium 5k geneset data for the mouse aorta. With transfer of single-cell RNA sequencing labels to the spatial clusters, we demonstrated that FMC-1 are enriched in the media and also localized in the cap after 16 weeks HFD. Subsequently high levels of FMC-2 were identified in the cap and spanning the intima to the CMCs at the base of the plaque, further validating our inference based on selected RNAScope markers (Fig. 1N-S, Supplementary Figure. 2A).

To quantitatively model these observed cell state changes, we applied the Waddington Optimal Transport (WOT) algorithm to build a probabilistic model for cell state transitions30,31. WOT is a heuristic method that models growth rates based on cell cycle and apoptosis gene expression to perform developmental trajectory inference. To visualize the WOT predicted relative transition probabilities and emphasize the actual proportion of each cell state at each time point, we have visualized the transition probability results with a Sanke plot (Fig. 2D). The exact proportion of each cell state at each time point can be observed in Fig. 2C and Supplementary Data. 5. Applying this algorithm across timepoints, we found that while SMC-1/2/3 are present at baseline (Week 0), SMC-2 are the primary cell phenotype that transitions to FMC states (Fig. 2D, E, F). WOT modeling suggested that a small FMC-1 population existed at baseline, significantly expanded by 5 weeks, and directly contributes to both FMC-2 and CMC. FMC-2 appeared to primarily become CMC which were evident by 12 weeks of HFD. However, given the higher proportion of FMC-1 present relative to FMC-2, both FMC cell states likely contributed similarly to the total number of CMC.

WOT further captured TF enrichment across cell state transitions, using weeks 5 and 12 as reference points which encapsulated the majority of the cell state transitions. This approach allowed us to highlight the crucial TFs expressed in SMC that guide these cells towards their final fates (Supplementary Data 6). We then filtered these transition TFs with significantly enriched accessible transcription factor binding motifs identified from scATACseq data. This integrative approach enabled us to identify TFs that have previously been associated with SMC phenotypic transitions, and in some cases pointed to previously unappreciated functions for these genes (Fig. 2G). For example, we found Tcf21 expression to be the most enriched TF in cells fated to become FMC-1, consistent with its early role in phenotypic transitions, while WOT-based TF enrichment also suggested a prominent role in the FMC-CMC state transition 32,33. Across the enriched TFs found in our model, several patterns stand out. For example, TFs (Runx1/2, Zbtb7c, Zeb1/2, Smad3, Stat3, Hes1, Nfkb1, Cebpb) amongst this list have been readily associated in literature with AP-1, Klfs and the SMC regulatory Srf-myocardin complex. This suggests that pioneer factors and immediate-early AP-1 activation constitute a core early event that then allows for recruitment of transcription factors at specific enhancers to drive key cell state transition steps in conjunction with Srf-Myocardin. These core TFs then organize around central cellular processes such as TGFB signaling, hypoxia signaling, inflammatory response and osteochondrogenesis as summarized in Supplementary Data 734-38. Moreover, the top SMC-3 promoting TF was Isl1, reinforcing its SHF origins39.

The suggested trajectory paradigm is further supported by analyses with the RealTime kernel calculations in CellRank that computationally derive sustained states. This analysis recovered FMC-1, FMC-2, two CMC and the SMC-2 and SMC-3 as sustained states (Fig. 2H). Further, we visualized the WOT derived growth rates by cell state and found a greater proportion of cells with increased proliferation in the early transitioning SMC-2 cells, reduced growth rates in the FMC with FMC-2 having higher proportion of low growth rates and higher senescence score (Fig. 2I, Supplementary Figure. 4B), and lastly a higher growth rate in CMC, reflective of proliferative early-stage chondrocytes in developing bone40,41 (Fig. 2I). Given the high correlation between pseudotime and cell state stages confirmed by our real-time analysis (Pearson’s r = 0.57; spearman’s rho = 0.56, p < 2e-16) (Fig. 2J), we modeled the atherosclerosis trajectory with CellRank pseudotime kernel as an orthogonal approach to derive a probabilistic transition model.

Analyses using WOT and pseudotime suggested that both FMC-1 and FMC-2 phenotypes transition probabilities stabilize in mature lesions and are thus more likely to represent sustained phenotypes rather than a transition state. Complementary analyses with pre-existing 20 and 26 week high fat diet (ApoE and Ldlr-KO, respectively) atherosclerosis mouse scRNAseq datasets15 and brachiocephalic derived advanced lesions11 support the sustained existence of these phenotype cells into late stages of disease (Supplementary Figure. 4F-K). This is an attractive possibility, following the hypothesis that protective cells transition primarily to stable FMC and disease promoting cells undergo further transition to the CMC lineage. These findings are in contrast to previous pseudotime based analysis which have suggested that CMC are the common primary endpoint for all transitioning SMC in the plaque9,12.

The increased resolution provided by pseudotime allowed us to cluster gene expression across the different inferred trajectories, identify key TFs within trajectory gene clusters, and characterize pathways based on gene expression trends clustered over pseudotime. For example, we found that FMC-1 trajectory clustered genes demonstrate early activation of genes involved in inflammatory response and response to molecules of bacterial origin suggestive of a core set of genes involved in stress response including TFs such as Klf4, Cebpb, Runx1 and Nfkb1 (Supplementary Figure. 5A, cluster 3) followed by gene groups demonstrating epithelial cell migration and cell-substrate adhesion including TFs like Tcf21, Ar, Zeb2, and Twist1 (Fig 2K) and more unique processes such as antigen processing and presentation and cytokine mediated signaling that includes TFs such as Ahr, Gas7 and Stat1 (Supplementary Figure. 5A, cluster 2). For FMC-2, we observed an early cluster enriched for epithelium migration and regulation of Wnt signaling (Fig. 2L, Supplementary Figure. 5B, cluster 3), while later clusters showed enrichment for cell chemotaxis, angiogenesis and collagen metabolic process (Supplementary Figure. 5B, clusters 1, 2). Similarly, for CMC, there were robust signals for biomineralization and ossification terms (Supplementary Figure. 5C).

Accessible chromatin reveals cell state specific DNA motifs
To characterize the epigenetic processes that mediate the noted transcriptional effects, we investigated chromatin accessibility in the transitioning cells with scATACseq data (Fig. 1C, D). We observed high specificity of chromatin accessibility within each label transferred tdTomato lineage traced cluster and visualized top specific peaks along pseudotime (Fig. 3A). Using GREAT for functional enrichment42, we identified specific cellular processes including smooth muscle related processes in the SMC, cell migration, inflammation and response to platelet derived growth factor in the FMC, and ossification and chondrocyte development in the CMC (Fig. 3A, Supplementary Data 8).

To prioritize dynamic TF binding motifs, we applied ChromVar to calculate their motif binding accessibility probability distributions. We filtered these data using a core set of overrepresented TF motifs identified by Signac and HOMER and visualized motif deviations along pseudotime to capture the dynamic shifts of TF binding activity (Fig. 3B). We used the early pseudotime bins 3-24 (PseudoEarly) and late pseudotime bins 25-30 (PseudoLate) to approximate the SMC-FMC and FMC-CMC transitions (Fig. 3A). These analyses identified temporal trends for both known and previously unrecognized regulators of smooth muscle fate. For example, Srf, Mef2, and Zeb motif deviation was observed to be higher in SMC13,43. In addition, for the early transitions, we observed a number of TFs with the greatest motif accessibility across the SMC stages, including Tead, Zbtb7c, Meox2, Rfx, and Nfi factors. Factors showing greatest motif accessibility during the late transition stages included Tcf21, Smad3, Rbpj, Stat, Nfatc, Cebpb, and Runx factors, with many of these known to be involved in endochondral bone formation. KLF factors and numerous AP-1 factors showed a bimodal pattern with accessibility early in the quiescent SMC state and later, from FMC to pre-CMC bins, likely representing their pioneer factor functions at these differentiated cell states. These analyses summarized the chromatin landscape guiding TF pathways, extending previous studies using pseudotime data derived from baseline and sustained phenotypic cell states13 (Fig. 3B).

Network prioritization identifies top SMC-FMC transition TFs
Using our co-embedded scRNA and scATAC dataset, we leveraged complementary network inference methods, Pando44 and CellOracle45, to create a custom workflow to infer transcription factor-target interaction networks (gene regulatory networks, GRNs) that direct phenotypic transition and simulate cell identity changes with in silico TF perturbations (Fig. 3C, Supplementary Data 9,10, Methods). The dataset was divided by PseudoEarly and PseudoLate bins to infer GRNs for these analyses. We visualized a summary GRN colored by the average pseudotime-TF expression to reveal a cascade of network activation (Fig. 3D). SMC lineage phenotypes clustered separately, identifying numerous TFs that are likely to direct specific regulons that mediate cell states and transitions characteristic of the response to disease stresses (Supplementary Data 11). We computed TF activities for the PseudoEarly and PseudoLate SMC states, and a comparison of TF activity change from SMC to transitioning SMC cell state activity revealed patterns such as shifts in hypoxia inducible factor (HIF) activity highlighted by the known interaction between Hif1a and Epas146 and other factors such as Ahr that interact with HIF through the common heterodimer Arnt47 (Fig. 3E).

We then performed systematic in silico perturbations using our inferred TF-target links from SMC-FMC and FMC-CMC transitions to calculate perturbation scores measuring the potential of TFs to drive cell transition away from the originating cell state. In SMC, we identified enrichment of EMT factors such as Tcf21 and Zeb factors, and similar to the top TF activity changes, we identified enrichment for hypoxia inducible factors including Hif1a, and Epas1, reflecting their high network connectivity in early phenotypic transitions (Fig. 3F). For factors promoting change to CMC phenotypes, we identified known factors such as Klf411 and Runx1/248,49, in addition to unique factors such as Trps1, Sox6, Erg, Zbtb7c, Prrx1, Arntl2, and Snai1 that were nominated to drive the differentiation of cells from FMC to CMC (Fig. 3G).

We further created an aggregated TF transition ranking by combining normalized enrichment from WOT for core TFs promoting the FMC-1 and CMC fates and normalized SMC-FMC as well as FMC-CMC perturbation scores from CellOracle (unfiltered WOT and CellOracle heatmaps in Supplementary Figure. 6A-C). Through this combination of cell fate enrichment and network-based prioritization analyses, we found Ar, Zeb, Smad, Mecom, Prrx1, and Cebpb factors, as well as Tcf21, to be among central enriched TFs involved in the phenotypic transition from SMC to FMC (Fig. 3H, Supplementary Figure. 6B), and Runx, Klf, Egr2, Sox9, and Zbtb7c as top drivers for FMC to CMC transition (Fig. 3I, Supplementary Figure. 6C). Additionally, Runx and Zbtb7c factors were consistently highly ranked across both transitions.

To validate our analytical algorithms, we selected three highly ranked transcription factors predicted from our CellOracle analysis to affect the transition process, AR, EPAS1 and ZEB1and performed siRNA knockdown followed by qPCR for contractile markers including ACTA2 and TAGLN as well as FMC markers LUM and PDGFRB (Supplementary Figure. 7A-C). Knockdown of these genes provided robust evidence for increased contractile marker expression across these conditions. Furthermore, we also observed significant downregulation of LUM and PDGFRB for EPAS1 and ZEB1 knockdown while AR knockdown showed a smaller effect size but a trend towards decreased expression levels of these FMC markers. Together, this provides additional functional validation of our predictive modeling.

To provide functional validation of the integrity of the FMC-1 versus FMC-2 cellular phenotypes, we characterized cell expression signatures from primary human coronary artery smooth muscle cells (HCASMC) derived from different donors previously published in Liu, et al.50. Using top marker genes from our single-cell mouse timecourse, we performed non-negative least squares cellular deconvolution and identified cell lines 2105 and 1508 as exhibiting the top ‘FMC-1-like’ and ‘FMC-2-like’ cell states, respectively (Supplementary Figure. 7D-G). We then treated these primary cell lines with calcium phosphate media14 for seven days to simulate an osteochondrogenic CMC-like transition phenotype. Quantitative PCR identified a greater increase in CMC genes including RUNX2, SOX9, and HAPLN1 in the 1508 cell line, whereas there was minimal change in expression of these genes in the 2105 line (Supplementary Figure. 7H, I). This is in keeping with our observations that the FMC-2 population has a greater probability of transition to the CMC state. Furthermore, this important finding suggests that there are functional differences in the primary cell lines depending on their origin/disease state and they may harbor differential ability to respond to various simulated vascular stresses. Identification of these primary HCASMC lines with FMC-1 vs FMC-2 phenotype will aid future studies of molecular differences between these cell states.

Tcf21 loss significantly alters SMC transition probabilities
As a top TF predicted to direct SMC to FMC transition, we further examined the effect of exemplar CAD gene Tcf21 knockout on SMC trajectories with updated single-cell chemistries using the timecourse methods employed for the control dataset to elucidate Tcf21 SMC regulatory mechanisms that occur with vascular stress.

We first examined changes in SMC transition cell numbers that resulted from Tcf21-KO. In keeping with our prior work16, there was a significant decrease in the FMC populations at 12 and 16 weeks of HFD (Fig. 4A). Comparison of lineage-traced SMC proportions revealed a marked ~3-fold increase in the SMC-3 Tcf21-KO cells at 5 weeks that persisted to week 12. These Tnnt2-expressing cells have been lineage traced to the secondary heart field21,22, suggesting an early expansion of this aortic base contractile medial compartment in the context of Tcf21 loss. The possible involvement of Tcf21 in the regulation of these cells was supported by its expression in Nkx2-5 lineage traced cells as identified with scRNASeq (Supplementary Figure. 8A-E). At 12 weeks, there was a notable relative decrease in Tcf21-KO FMC-1, FMC-2 and CMC proportions that corresponded with a relative increase in SMC-2 and SMC-3 cluster proportions, suggesting that the Tcf21-KO cells were halted in disease associated transitions. At 16 weeks, there was a consistent decrease in FMC-1 without a decrease in FMC-2, suggesting potential compensation for the loss of Tcf21 expression in the SMC lineage cells over time (Fig. 4A). At both 12 and 16 weeks we also observed a decrease in CMC cell proportion in the Tcf21-KO relative to control, suggesting that Tcf21 directly promotes CMC development as suggested by our previous analyses (Figs. 2G, 3I).

WOT was used to investigate alterations in SMC trajectory transition probabilities in the context of Tcf21-KO (Fig. 4B). Focusing on cell state changes in the SMC lineage traced cells, there was a notable decrease in the overall transition probabilities across cell states for Tcf21-KO. The SMC-2 contribution to both FMC-1 and FMC-2 was decreased in KO diseased tissues, as was the FMC-1 contribution to the FMC-2 phenotype, and there was a low level of transition for FMC-1 and FMC-2 to CMC (Fig. 4B). Overall, there was evidence for a dramatic decrease in probability of transition to the CMC phenotype while transition to both the FMC-1 and FMC-2 phenotype was decreased but sustained. These decreased probabilities for SMC to FMC and FMC to CMC transitions in knockout mice accounted for the changes in SMC phenotype cell numbers but also highlighted the presence of compensatory processes which allowed continued phenotypic transition (Fig. 4A). As noted previously, at 16 weeks, the FMC-1 and FMC-2 states remained, suggesting that they represent a sustained phenotype rather than exclusively a transition to the CMC phenotype.

Tcf21-KO DEGs enrich for predicted network and CAD genes
Analysis of differentially expressed genes (DEGs) in the Tcf21-KO compared to control mice using DESeq2 provided insight into associated gene programs. DEGs for the FMC-CMC transition for Tcf21-KO versus control, identified 965 genes (Supplementary Data 12, 13). This list was enriched for TGFB family genes such as Ltbp1 and Tgfb2, and numerous CAD GWAS genes including Tgfb1, Zeb2, Lrp1, Palld, Col4a2, Lmod1, and Pdgfd (Fig. 4C). These DEGs were further analyzed with GSEA using GO-BP terms to gain insight into altered pathways and directionality of effect. We found Cellular response to Tgfb stimulus and actin filament bundle organization to have high positive DEG enrichment, consistent with Tcf21 promoted effects and increased representation of locomotion and wound healing possibly representing compensatory processes given the overall decreased phenotypic transition phenotype (Fig. 4D). Terms with a negative enrichment score (average downregulation with Tcf21-KO) were identified as suppressive for connective tissue development and endochondral bone morphogenesis, further suggesting that Tcf21 target genes likely have a role promoting the CMC phenotype.

SMC lineage cell state changes and related cellular trajectories can be modeled through gene-gene interactions in a GRN. To build a Tcf21 specific network we employed the regulatory network predicted from our PseudoLate control timecourse which spanned the majority of Tcf21 activation in transitioning SMC. Moreover, network targets by this method are more likely to represent direct interactions given the conditional need for epigenetically accessible Tcf21 binding sites to be linked to target genes. We compared Tcf21 target genes identified with this network analysis with Tcf21-KO DEGs and found significant overlap of DE genes and network genes, with 135 of 318 predicted network genes (42%) also showing differential expression with Tcf21-KO (Fishers exact test p<1e-4) (Fig. 4E). Further, we applied GSEA using the predicted Tcf21 network gene set ranked by TF-gene regulatory coefficient weights and observed a significant normalized enrichment score (NES = 1.24, p=0.031) of Tcf21-KO DEGs. These analyses showed congruence between distinct approaches and validated the utility of using predictive transcriptional modules from a comprehensive control dataset to infer perturbed pathways.

We leveraged the predictive pathway ability of inferred networks to augment the functional characterization of the Tcf21 signaling network. We selected differentially expressed TFs within the Tcf21 GRN and integrated these TF-centered GRNs to create a validated Tcf21-TF sub-network and find both repressed and activated TFs that are regulated by Tcf21 (Fig. 4F). We identified downstream TF networks with differential KO module scores and performed functional enrichment on this subset of TF networks to visualize the pathways affected. We found multiple TFs including Foxp2, Mecom, Zeb2, Meis2, Mecom and Tshz2 predicted to be involved in cell migration, angiogenesis and extracellular organization processes, while genes such as Meis1, Zeb2 and Foxp2 were also involved in contractile processes (Fig. 4G).

We computed TF fate correlations by comparing WOT transition probability with cell-level TF module scores and using TF-only expression as a control comparison (Fig. 4H). Interestingly, top modules that exhibit high correlation with FMC fates also have their central TF upregulated with Tcf21-KO, suggesting compensatory roles alongside Tcf21 given our observation of overall decreased FMC proportions. For example Zeb2, previously shown to drive phenotypic transition13, was increased in Tcf21-KO, and its later average TF-pseudotime suggested that its role is downstream of Tcf21 towards the FMC-2 fate. In contrast, TFs such as Meis1/2 have earlier average TF-pseudotime, suggesting compensatory upregulation possibly via feedback mechanisms. Conversely, top TF modules correlated with the CMC fate, including the known ossification regulator Sox9, showed decreased expression upon Tcf21-KO (Fig. 4I).

TCF21-TEAD epigenetic interactions modify CAD GWAS genes
Previous studies have shown that Tcf21 interacts with histone deacetylases to broadly alter the in vitro epigenetic landscape of human coronary artery SMC (HCASMC)51 but its effects in vivo have not been explored. We observed widespread motif accessibility deviations upon Tcf21-KO when visualizing differential ChromVar scores across pseudotime and generated differential ChromVar scores on a by cluster basis (Fig. 5A, Supplementary Data 14,15). We further used Jensen-Shannon divergence (JSD) scores to identify differentially deviated TF motifs across pseudotime between control and Tcf21-KO ChromVar distributions, finding significant differences in many core TF motifs including Zeb, Klf, Runx1/2, AP-1, Srf, and Ctcf, while Tcf21 was borderline significant (padj = 0.10) (Fig. 5B). TCF21 HCASMC ChIPseq was reprocessed from Zhao et al.52 and showed significant enrichment in TEAD and CEBP TF binding motifs within TCF21 peaks (Fig 5C, Supplementary Figure. 9A) in addition to TCF21 and AP-151. When these ChromVar scores were visualized along pseudotime in the mouse timecourse, there was enhancement of Tead1 motif accessibility with loss of Tcf21 while Cebpb shared a pattern similar to Tcf21, showing decreased accessibility with Tcf21-KO (Supplementary Figure. 9B). We performed additional ChIPseq in HCASMC for TEAD1 and confirmed significant overlap with TCF21 (Fig. 5C, Fisher’s exact test p<1e-4). Pathway analysis of shared peaks using GREAT predicted shared biological functions related to inflammation, apoptosis, TNF, and TGFB (Fig. 5D). Peaks shared by TCF21 and CEBPB were enriched for cell adhesion, ERK signaling, and cell motility terms (Supplementary Figure. 9C, D, Fishers exact test p<1e-4). Genome wide colocalization of TCF21 and CEBPB/TEAD1 motifs was identified by comparing the location of CEBPB/TEAD1 motifs in TCF21 ChIPseq peaks and TCF21 motifs in TEAD1 peaks (Supplementary Figure. 9D, E). Interestingly, performing enrichment analysis for GWAS SNPs at TCF21-TEAD1 co-bound loci using ChIPseq binding data and GWAS SNP localization with the GWASAnalytics package, we found that all TCF21 binding sites, including those with TEAD1 peaks, showed high level enrichment for CAD (-log p-value 10) (Fig. 5E, Supplementary Figure. 9F), but in the absence of TEAD1 peaks showed only low level enrichment for cancer and metabolism (Supplementary Figure. 9I). Also, TEAD1 peaks including TCF21 peaks showed strong enrichment for hypertension (-log p-value 10), while TEAD1 peaks only showed low level enrichment of SNPs for diabetes and hypertension (Supplementary Figure. 9G, H). Taken together, these data suggest that the interaction of TCF21 and TEAD1 plays a significant specific role toward CAD risk. This TCF21-TEAD1 relationship corresponded with murine trajectory analysis nominating Tead1 as an FMC-2 driver (Fig. 2L) and making it an intriguing Tcf21-interacting partner for further study. These results for TCF21-TEAD1 were contrasted with similar studies for CEBPB where inclusion of TCF21 loci detracted from the CEBPB metabolism signal (Supplementary Figure. 9J, K). Analyses for HNF1A peaks served as negative control for these studies (Supplementary Figure. 9L, M).

We further examined the shared genomic patterns between TCF21 and TEAD1 by partitioning their shared binding loci into TCF21+TEAD1 or TCF21 only loci. We observed greater enhancer profiles for TCF21+TEAD1 shared binding sites compared to TCF21 only, as indicated by overlap with H3K27ac (Supplementary Figure. 9N). TCF21+TEAD1 shared binding sites were located farther from the TSS regions (Supplementary Figure. 9O), consistent with enhancer co-localization. Further, pathway analysis of putative genes identified in TEAD1/TCF21/H3K27 regions by GREAT revealed enriched pathway keywords including differentiation, development and endopeptidase regulation while TCF21-only+H3K27 region genes showed enrichment for immune, viral and neutrophil related keywords (Supplementary Figure. 9P, Q).

We then investigated physical and functional interaction of these two TFs. Proximity ligation assays found that TCF21 and TEAD1 co-localized in the nucleus, suggesting direct protein-protein interaction (Fig. 5F). To detect a direct physical interaction between TCF21 and TEAD1, we performed Co-IP using a myc-tagged TCF21 transfected into HEK293 cells (Fig. 5G). We performed nuclear protein extraction followed by IP for the MYC-tag and western blot for TEAD1. These studies included IgG negative control and 5% input positive controls and provided evidence for TCF21-TEAD1 physical interaction. Further, immunohistochemistry in mice aortic root atherosclerosis sections also demonstrated intimal expression of Tead1 and Tcf21 (Supplementary Figure. 10A-C).

Given these findings, we investigated the expressed gene programs in SMC directed by TEAD1 by performing TEAD1 knockdown with siRNA in HCASMC along with bulk RNAseq. The transcriptomic changes with TEAD1 knockdown showed strong enrichment for processes such as cellular response to cAMP, extracellular matrix organization, connective tissue development, and chondrocyte development (Fig. 5H). These findings are in line with what we observed from our TEAD1 ChIPseq functional enrichment, further corroborating the hypothesis that TEAD1 plays a major role in the smooth muscle transition process to affect the development of phenotypically transitioning SMC.

Finally, to examine the functional interactions between TCF21 and TEAD1, we performed dual luciferase reporter gene transfection assays with A7r5 rat smooth muscle cells on a shared enhancer residing in an intron of the SRF gene, a master regulator of lineage contractile gene expression35, and two additional enhancers in CAD loci encoding ECM effectors of TGFB signaling, BMP1 and LOXL1. For the SRF enhancer, we showed that normal activation by SRF binding partner MYOCD was highly suppressed by TCF21 and to a greater degree by TEAD1 alone (Fig. 5I). There was an intermediate reporter activity when TCF21 and TEAD1 were both transfected, suggesting a competitive interaction between these TFs. Also, for the BMP1 and LOXL1 enhancers, both TFs showed repressor activity, but again, intermediate suppression when both were expressed in the same cells, suggesting competition (Fig. 5J, K). Taken together these data suggest that TCF21 and TEAD1 directly interact at shared regions across the genome to epigenetically regulate transcription.

Tcf21 mediates SMC CAD genetic risk via early rewiring
Single-cell methods have previously nominated genetic risk signals to have a unique high enrichment in SMC8. We investigated the relative disease related significance of our SMC cell states, using the scDRS algorithm5,53. At the single-cell level, we leveraged scDRS to integrate gene expression and GWAS gene z-score weights from MAGMA to generate disease relevance scores for each cell type. Because scDRS quantifies risk gene enrichment but not directionality, we also computed GWAS risk-weighted average expression by taking average individual gene expression multiplied by its scDRS gene weights and further inferred directionality using updated heritability adjusted S-PrediXcan modeling54 (Supplementary Data 16). Using this framework, we identified FMC-2 having overall statistical enrichment for excess CAD-associated risk gene expression by scDRS (Fig. 6A).

Interestingly, neither the transcriptionally similar FMC-1 nor the calcification-associated CMC exhibited significant scDRS enrichment. Next, from the aggregated GWAS risk-weighted expression analysis, of genes with available predicted risk directionality, we calculated greater averaged expression of CAD risk genes in CMC relative to FMC and SMCs. Consistent with this complex biology, FMC-2 cells express genes that both promote or suppress CAD risk with the net predicted risk direction more protective, prompting further investigation of distinguishing expression patterns of these transition phenotypes.

We used the TF networks identified from the control timecourse (Methods) to ask how individual TFs and their networks associate with CAD risk. First, we synthesized full TF network modules by assimilating all unique target connections from PseudoEarly and PseudoLate predicted regulatory networks for available TFs. Second, we extracted gene-level scDRS weights to derive normalized TF-scDRS correlations. Third, we averaged gene-scDRS correlation ranks of all genes within each TF network module to create a TF network average scDRS rank. At the network level, we observed TFs with higher average scDRS rank in the SMC to FMC pseudotime such as Tcf21, Nkfb1, Runx, Zeb, Hif and Smad factors (Fig. 6B, Supplementary Data 17). Many of these core TFs also overlap with in silico TF perturbation predictions (Fig. 3F, G). Moreover, this method nominated novel TFs with network enrichment of CAD risk genes that were postulated to drive the SMC-FMC transition, including Arntl, Prrx1, Tshz2, and Mecom (Supplementary Data 8) or the FMC-CMC transition such as Trps1, Zbtb7c, Snai1, and Sox5/6/9.

To further dissect the regulatory relationships of CAD GWAS genes in phenotypic transition, we examined the scDRS enrichment in Tcf21-KO and found an increase in SMC-3 scDRS z-score, meeting nominal significance (Fig. 6C). This shift in scDRS score implicated the ability of Tcf21 to coordinate CAD GWAS genes early in the phenotypic transition timeline. In addition, this observation was in agreement with Tcf21 expression showing SMC-3 enrichment (Fig. 1C) and demonstrating a basal level of accessible Tcf21 TF binding sites in early pseudotime (Fig. 3B). We then examined the DEGs in PseudoEarly bins and observed significant overlap with a curated set of putative GWAS genes (63/640, p < 0.0001) (Fig. 6D). GSEA enrichment of DEGs not only revealed similar increased contractile processes with Tcf21-KO but also highlighted a decreased response to stress signals such as ‘cytokine stimulus’ and ‘unfolded protein response’ (Fig. 6E). Integrating these genes with human genetic signals, we focused on the differential expression of genes that shared GWAS disease risk-eQTL correlation with Tcf21. For example, increases in proliferative factors such as PDGFD and SEMA3C or decreases in inflammatory transcription factor STAT3 were associated with increased CAD risk (Fig 6. F-H). Conversely, there was also an enrichment of genes that were correlated with decreased CAD risk such as LRP1 which has multifaceted coronary disease implications, COL4A2 which promotes basement membrane integrity or MYO9B that modulates vascular wound repair (Fig 6. I-K). Together, these relationships suggest a broad role for TCF21 in promoting risk through coordinated regulation of CAD GWAS genes in the phenotypic transition of disease SMC.

Methods
Mouse strains, induction of lineage marker, and sample collection
Our research complies with all relevant ethical regulations, with our animal study protocols (Protocol ID -10020, 10054) approved by the Institutional Animal Care and Use Committee at Stanford University. Control (final genotype - TgMyh11-CreERT2, Tcf21+/+, ROSAtdT/+, ApoE−/−) and Tcf21 (final genotype - TgMyh11-CreERT2, Tcf21ΔSMC/ΔSMC, ROSAtdT/+, ApoE−/−) knockout mice with floxed tdTomato fluorescent reporter to allow for SMC-specific lineage tracing were bred onto a C57BL/6 ApoE-/- background as previously described16. For all subsequent lineage tracing experiments, two doses of tamoxifen 48 hours apart via gavage was carried out when the mice reached 7.5 weeks of age. In the Tcf21-flox group, tamoxifen gavage also induces the Tcf21 knockout in addition to tdTomato lineage tracing. Following tamoxifen treatment, all mice were placed on a high-fat diet. At the designated timepoints, animals are euthanized after CO administration followed by cervical dislocation, consistent with the recommendations of the Panel on Euthanasia of the American Veterinary Medical Association (AVMA).The study involved two primary high fat diet experimental arms: one for single-cell RNA sequencing (scRNA) and the other for Assay for Transposase-Accessible Chromatin sequencing (ATAC).

Animals were housed in an AAALAC-accredited, specific-pathogen–free barrier facility in individually ventilated cages under a 12-hour light/12-hour dark cycle (lights on 0700–1900). Ambient temperature was maintained at 20-26°C with relative humidity at 30-70% per institutional standards. Sterile bedding was changed regularly, and mice had ad libitum access to standard rodent chow or high fat diet (21% anhydrous milk fat and 0.15% cholesterol (Dyets no. 101511; Dyets) and water. Environmental enrichment (e.g., nesting material and shelters) was provided.

For scRNA experiments, samples were collected from Control mice at multiple time points named by number of weeks on high fat diet: 0, 3, 5, 7, 9, 12, and 16 weeks. Samples from the Tcf21 knockout mice were collected at 5, 12, and 16 weeks. For ATAC experiments, samples were collected from the Control mice on the high-fat diet at 0, 5, 7, 9, 12, and 16 weeks and for Tcf21 knockout mice at 5, 12, and 16 weeks. 3-4 male mice are pooled for each single-cell collection as the Myh11CreERT2 transgene is inserted in the Y chromosome, limiting lineage tracing studies to the male sex.

Aortic digestion for 10x Genomics microfluidics
Samples from mouse aorta were dissociated into single cells for RNA and ATAC sequencing using the 10x Genomics Chromium platform. Euthanized mice are perfused with phosphate buffered saline and dissected to obtain the aortic root up to the level of the brachiocephalic artery. The tissue is washed with PBS and incubated in an enzymatic dissociation cocktail containing Liberase (5401127001; Sigma-Aldrich) and Elastase (LS002279; Worthington) in Hank’s Balanced Salt Solution with calcium (HBSS+ Ca2+) for 30 minutes. Mechanical dissociation is performed for 5 minutes followed by visualization under the microscope to ensure single-cell suspension. This suspension is strained through a 35µm nylon mesh snap cap into falcon test tubes (352235 BD) followed by FACS sorting (Sony SH800) for tdTomato positive and negative cells in parallel. For single-cell ATAC, nuclei was isolated per 10X recommended protocol and captured on the 10X scATAC platform. Each individually sorted cell suspension was loaded onto 10x GEM G orH chips with remainder per 10x protocols (Chromium Single Cell 3’ RNA V3.1, Chromium Single Cell ATAC V2). For both scRNA and scATAC runs, samples at some intermediate time points were processed as pooled samples with tdTomato and non-tdTomato loaded at 1:1 ratio after FACS and listed in Supplementary Data 1. Libraries were sequenced on the Illumina NovaSeq6000 platform with targeted depth of 40-50,000 reads per cell for RNA and 75,000 reads/cell for ATAC.

10x RNA and ATAC Data Preprocessing
RNA fastq files were processed using Cellranger v7.0.1 (10x Genomics) to obtain transcript count matrices and aligned to mouse transcriptome mm10-2020-A-2.0.0 (10x Genomics) with custom addition of lineage tracing tdTomato transcript. Samples at each timepoint were aggregated and analyzed with R package with Seurat (v4.3)68. Low-quality cells and mitochondrial-rich cells were filtered with parameters mitochondrial percentage < 6%, ribosomal percentage < 25%, and nFeature_RNA > 1250 and < 8000. Gene expression count matrices underwent log-transformation and library-specific scaling. Importantly, no additional batch correction for visualization was required in the pre-processing steps given the uniform processing of mice samples. Principal component analysis was used for dimensionality reduction followed by clustering using the Louvain algorithm.

For the full timecourse dataset, aligned scRNA files from Control and Tcf21-KO data were merged and processed with identical QC parameters as above. Clustering was then performed on this merged dataset in order to ensure comparable differential gene analysis and pseudotime comparisons. The processed dataset is then split into Control and Tcf21-KO objects for independent analysis. We further applied logistic regression to determine the optimal tdTomato expression cutoff and subset on SMC-derived lineage cells exclusive to tdTomato fluorescence-activated cell sorting (FACS).

For force directed layouts, we followed methods as described in Schiebinger et al. 50 dimensional diffusion components are calculated with SCANPY (1.9.3) with default parameters. For each cell, its 20 nearest neighbors was used to produce a nearest neighbor map. Then we applied leiden clustering at a resolution of 0.36 and generated cluster connectivities via PAGA and visualized the force-directed layout on the k-NN graph using ForceAtlas2.

scATAC Raw fastq files were uniformly processed using Cellranger-atac-2.1.0 (10x Genomics) and aligned to cellranger-arc-mm10-2020-A-2.0.0 reference genome (10x Genomics). The data were processed to remove low-quality cells and reads with low mapping quality with parameters peak_region_fragments > 2500, peak_region_fragments < 100000, pct_reads_in_peaks > 30, nucleosome_signal <= 4, TSS.enrichment > 2. The remaining reads were then used to construct a sparse binary matrix, representing chromatin accessibility states across individual cells and genomic regions. A unified set of peaks from all samples (Control and Tcf21-KO) were merged using overlapping and adjacent peaks using CellRanger-ATAC aggr protocol. This unified peak set was re-quantified by term frequency-inverse document frequency (TF-IDF) using Signac69 (1.10) RunTFIDF.

To preprocess scRNA and scATAC data for integration, gene activities are first calculated from chromatin accessibility data using GeneActivity() function from Signac with default parameters and log normalization. Datasets were integrated using canonical correlation analysis (CCA) with the Seurat RunCCA() function upon 2000 features using the SelectIntegrationFeatures() in Seurat. We then use the integrated dataset to perform pseudotime ordering using Slingshot70 (2.6) and split the dataset into 30 pseudotime bins for uniform downstream analyses.

scATAC Functional and Motif Analyses
To generate identify overrepresented TF binding motifs, we performed hypergeometric enrichment on the top 20,000 variable peaks with GC-matched controls using the FindMotifs function in Signac. We then further merged this motif set with HOMER’s (4.10) de novo motif enrichment method on a per-cluster basis taking the top overrepresented motifs along with all similar motifs above a similarity score of 0.7. This allowed identification of a comprehensive set of enriched transcription factor motifs to aid summarization of ATAC data. For mouse ATAC data, the background peak set utilizing the entire merged peak set generated by CellRanger-ATAC (version 2.0.1) Aggr.

We then calculate a cluster specificity score by dividing the detection percentage of accessible chromatin for each cluster by the detection percentage of all other clusters with a minimum 10% cutoff for within cluster accessibility. We filtered peaks with a cluster specificity of >1.5 for FMC-1, FMC-2, and CMC and specificity > 1.25 for SMCs. From these peaks, we selected the top 5,000 peaks with the highest specificity score in the FMC-1, FMC-2, and CMC clusters while for the aggregated SMC cluster (SMC-1, SMC-2, SMC-3), we use the top 2500 peaks. ChromVAR71 was used to calculate the transcription factor specific enrichment scores of accessible motifs across each pseudotime bin through the Signac wrapper RunChromVAR.

Trajectory Analysis with Waddington Optimal Transport and CellRank2
To model and infer cellular trajectories in our data, we employed the Waddington optimal transport algorithm30. This approach infers cell growth rates using single-cell gene expression to generate a transition probability distribution, enabling the identification of cellular transitions and differentiation paths.

We applied the default command line implementation of Waddington WOT as described in Schiebinger et al. with cell scores derived from updated KEGG cell cycle and Apoptosis gene signatures for estimation of cell growth and death. Growth rate tables are extracted and added to single-cell metadata for additional analyses. We generated a ‘fate correlation by comparing WOT predicted cell transition probability for the FMC-1, FMC-2 and CMC fate with each Tcf21 downstream TF network module score in the combined Tcf21-KO and control dataset. As negative control, we also correlated respective TF expression with fate probability to identify consistently more robust correlation with TF module scores compared to its central TF. Fate tables in Supplementary Data 6 are arranged in order of fraction expressed ratio.

Orthogonal trajectory analyses was performed using CellRank2. We applied the RealTimeKernel using default settings to identify sustained cell states and calculate driver gene trends plotted across our imported Slingshot derived pseudotime.

Data integration, gene regulatory network generation and in silico transcription factor perturbation
We utilize the Pando44 (1.0.1) to generate gene regulatory networks and CellOracle45 (0.18.0) for in silico TF perturbation experiments. By combining Pando’s network construction capabilities with CellOracle's in silico perturbation analysis, this approach enables a comprehensive exploration of gene regulatory dynamics in single-cell data. We divided the integrated dataset into PseudoEarly (pseudotime bins 3-24) and PseudoLate (pseudotime bins 25-30) for further analyses.

To generate GRNs from our scRNA and scATAC integrated dataset, we identified transcription factors and their respective binding sites in regions of accessible chromatin through motif scanning, and these transcription factor module candidate regions linked to genes by proximity through the PANDO framework. The bagging regression model was selected to infer the relationships between TF expression, binding-site accessibility and target gene expression to generate cell state specific networks.

We utilized the extended transcription factor motif database from the original Pando manuscript (including JASPAR202072, CIS-BTP73 as well as TFs without known motifs assigned by sequence homology) and further included all JASPAR2020 database mouse reference motifs for downstream motif scanning. Gene regulatory network inference was then performed using default PANDO parameters for selection of candidate regulatory regions from scATAC data, transcription factor motif scanning, selection of region-TF pairs, while the final regression model, we substitute the bagging ridge model for the default generalized linear model to match that of CellOracle. The output coefficient table was extracted and filtered for adjusted p-value < 0.05, R2 > 0.1, minimum number of variables (region-TF pairs) per model > 10, and minimum genes per module > 5 to generate the final regulatory graph that is provided as input for downstream analysis including:
1. Weighted pseudotime x TF module UMAP, is calculated by the get_network_graph() function in Pando to generate a subgraph with significant module TFs as features to generate a UMAP embedding which takes into account the ingoing connection of each node (TF) as well as coexpression of TFs. Additionally, the product of average pseudotime per cell by TF expression is colorized onto each TF module to provide a visual representation of the TF module relationships in the context of pseudotime.
2. TF activity score – calculated by the averaging the sum of TF expression x target gene coefficient for each TF module in the early (pseudotime bin 3-24) or late (pseudotime bin 25-30) SMC gene regulatory network.

For in silico perturbation, the merged control timecourse integrated dataset is downsampled to 40,000 cells to allow for optimal computation speed along with the top 3000 most variable genes calculated in Seurat and extended TF database from Pando converted into CellOracle compatible format to generate a baseline ‘Oracle’ file. The Pando derived TF-gene regulatory network table is converted into ‘links’ network file format as input for CellOracle’s link_data parameter for in silico systemic TF perturbation modeling with negative PS sums visualized as described in the CellOracle tutorial.

Xenium slide processing and analysis
Fresh frozen tissue sections were prepared following the 10x Genomics Xenium In Situ for Fresh Frozen Tissues Protocol (CG000579 Rev F). Briefly, tissues arranged in order of aortic root, ascending aorta, brachiocephalic artery with right subclavian branch, per-diaphragmatic descending, and abdominal descending aorta from a control (TgMyh11-CreERT2, Tcf21+/+, ROSAtdT/+, ApoE−/−) mouse were embedded in OCT compound immediately after dissection. Embedded tissue blocks were frozen on dry ice and immediately stored at -80C prior to sectioning. Sections of 10m thickness were cut using a Leica CM1860 cryostat. Sections were mounted within the sample area (10.45 x 22.45mm) of Xenium slides (PN-3000941) without overlap and stored at -80C until ready for use. All subsequent Xenium steps were processed at the Stanford Functional Genomics Core. Tissue morphology and quality was assessed using Hematoxylin and Eosin (H&E) and DAPI staining prior to downstream Xenium In Situ Gene Expression assays.

Processed Xenium datasets were imported into Seurat V5 for normalization and integration. Feature-barcode matrices were normalized using the SCTransform workflow with spatial coordinates retained for downstream deconvolution. To perform single-cell reference label transfer to the Xenium data, canonical correlation analysis (CCA) was implemented following strategies from the Seurat tutorial “Analysis, visualization, and integration of spatial datasets with Seurat” (https://satijalab.org/seurat/articles/spatial_vignette.html). RCTD (Robust Cell Type Decomposition) was applied to further deconvolve spot-level spatial data based on input cell reference.

Mouse to human scRNA label transfer and Spatial Slideseq analysis
Data is obtained and re-processed from the CZI Arterial Atlas courtesy of Zhao et al.23. We processed single-cell RNA data based on parameters published in Zhao et al. and subset the SMC clusters. We then transferred labels from our control smooth muscle control dataset and transferred labels using CCA with 1:30 dim and assigned SMC cell states to the human scRNA data. Slide-Seq object was aligned to hg38 by curioseeker pipeline and again underwent processing as described in Zhao et al.23. Robust Cell Type Decomposition was then applied to integrate this label transferred reference with spot-level data from the spatial transcriptomic sequencing using Seurat V574. Visualization was performed using Seurat built-in functions including SpatialDimPlot to visualize the label transferred cell states.

Processing of external datasets, Pan et al., Alencar et al., Cheng et al.
Data from Pan et al. (GSE155513)15, Alencar et al (GSE150644)11, Cheng et al. (PRJNA794806)12 were downloaded and reprocessed using standard 10x CellRanger pipeline as noted above to facilitate label transfer analysis with our timecourse dataset. All subsetting parameters and clustering parameters decisions were made based on available published methods, following prior settings as closely as possible. For the Pan et al. dataset, we performed further data pruning, subsetting based on fluorescent marker expression as previously described in Sharma et al.19. We then applied cutoffs of nFeature > 1900, percent.mt < 7.5 to match our control timecourse study. For clustering, we identified top 1500 variable genes and generated UMAP with 20 principal components (PCs) as previously described. For Alencar et al. we filtered on nFeature_RNA > 200 and < 5000, percent.mt < 10. 19 PCswere selected for clustering. The clustering resolution was not specified. This does not affect our downstream analysis, as the goal for this analysis is to understand the distribution of label-transfer clusters from our control timecourse study. For Cheng et al. we requested original processed script from the author to recreate the R-Smad clustering. After clustering of these datasets, we then performed CCA label transfer followed by generation of confusion matrices to compare and contrast the relationship between clusters and identify existing cell states.

ChIPseq analyses
Fastq files were mapped to hg38/GRCh37 genome with Bowtie2 (1.2.3). ChIP-seq peaks were then called with MACS2 (2.2.7.1) using default parameters. From this output, ‘robust’ peaks were selected by specifying a minimum fold-enrichment of 5.

GWAS trait SNP enrichment analyses
The intersection of GWAS loci and transcription factor binding (TEAD1 + TCF21) was defined as the SNPs located within any overlapping region with ChIP-seq peaks. Direct overlap of SNPs assimilated from GWAS Catalog + MVP CAD GWAS was performed with GWASanalytics script (https://github.com/zhaoshuoxp/GWASanalytics). The binomial overlap performed by this script has been described previously75.

scDRS (Single Cell Disease Relevance Score)
The scDRS53 algorithm calculates a disease relevance score for each cell by comparing its gene expression pattern with a reference disease gene signature obtained from the MVP coronary artery disease GWAS data. We utilized summary statistics obtained from the coronary artery disease GWAS meta-analysis from the Million Veterans Program (MVP) which also incorporates transancestry genetic data. MVP CAD GWAS was munged using (MAGMA) to obtain a weighted CAD-associated gene list that is applied towards the mouse timecourse data. Remainder scDRS analysis was performed as described in Zhang et al.53.

Bulk RNASeq Deconvolution
Bulk RNAseq fastq files were downloaded from GSE113348 and processed through a uniform pipeline with cutadapt (5.0), STAR (2.7.10b) and featureCounts (2.0.6) for trimming, alignment and generation of count matrices, respectively. Data was then converted to TPM with Kallisto (0.51.1) for comparison between cell lines. Following the prior publication, we removed nine cell lines from the original dataset due to poor data quality for a total of 52 cell lines used. We then applied non-negative least squares deconvolution using the python implementation of Nnnls from scipy to estimate the contributions of SMC-derived phenotypic clusters based on our tdTomato SMC subset reference. The top 100 marker genes per cluster (ranked by pct.1/pct.2 with at least 25% expression in the reference cluster) were used as features in the cell type signature matrix.

Curated CAD GWAS gene list
A curated list of nominated GWAS genes based on prior CAD GWAS meta-analysis was tabulated from Erdmann 2018, Koyama 2019, Matsunaga 2020, Tcheandjieu 2022, Aragam 202276 5 6 77,78.

Heritability adjusted PrediXcan, LocusZoom and eQTL correlation visualization
We estimated gene expression risk directionality inferred from the updated heritability adjusted S-PrediXcan modeling54. LocusZoomR79 (0.3.8) was used to visualize gene locus plots. eQTL colocalization was performed by extracting positional coordinates from MVP CAD GWAS and retrieving SNPs from dbGaP. SNP with lowest p-value near nominated GWAS gene was selected and all SNPs meeting GWAS significance 5e-8 and within 50kb up and downstream of lead SNP were selected. We then identify which SNPs had corresponding eQTLs using the V10 GTEx release. The beta for each SNP was correlated with the corresponding NES of each eQTL and plotted as a scatter graph with color representing the R2 (calculated with LDlinkR) with lead SNP.

In situ RNA hybridization (RNAScope)
Slides were processed according to the manufacturer’s instructions, using reagents from ACD Bio (ACD 322360-USM). Slides were washed in PBS, then immersed in 1× Target Retrieval reagent at 100°C for 5 minutes. After washing twice in deionized water, slides were immersed in 100% ethanol, air-dried, and sections were encircled with a liquid-blocking pen. Sections were incubated with Protease III reagent at 40°C for 30 minutes, then washed twice with deionized water. Sections were incubated with probes against Fbln1 (ACD 502881), Fbln2 (ACD 447931), Ibsp (ACD 415501), C3 (ACD 417841), Tcf21 (ACD 508661), Thbs1 (ACD 457891), Notch3(ACD 425171), Col8a1 (ACD 518071), Ttc9 (ACD 1113921-C1), Loxl1 (ACD 492531), Bmp1 (ACD 311151), Vcam1 (ACD 438641), Col2a1 (ACD 407221) or a negative control probe (ACD 310043) for 3 hours at 40°C. Multiplex fluorescence and colorimetric assays were performed per the manufacturer’s instructions.

HCASMC and A7r5 cell culture
Primary HCASMCs and HCASMC-hTERT80 were maintained in SMC basal medium (SmBM) supplemented with SmGM-2 SingleQuots kit (Lonza CC-3182) including human epidermal growth factor, insulin, human basic fibroblast growth factor and 5% FBS, according to the manufacturer’s instructions.

Rat aortic smooth muscle cells (A7r5) were purchased from ATCC and cultured in Dulbecco’s Modified Eagle Medium (DMEM) high glucose (Fisher Scientific, #MT10013CV) with 10% FBS at 37 °C and 5% CO2. A7r5 at passage 6-9 were used for experiments.

Calcification Assay
Primary cell lines (1508 and 2105) were split into 6 well plates at 60k cells per well. Cells were allowed to equilibrate over 12 hours prior to media change with SmBM with supplements as noted above. After 24 hours, once cells are noted to be confluent, we applied a calcification media cocktail consisting of SmBM basal media supplemented with 0.4% FBS, (Lonza CC-3182 without additives). Cells were incubated in calcification media for 3 days followed by media change with calcification media and allowed to grow to 7 days prior to RNA collection as previously described14.

siRNA Experiments
Vascular SMCs were transfected with siRNAs targeting TEAD1, AR, EPAS1, and ZEB1(Dharmacon ON-TARGETplus SMARTpool siRNA L-012603-00-0005, L-003400-00-0005, L-004814-00-0005, L-006564-01-0010; Horizon Discovery). Silencer Select negative control siRNA (ThermoFisher 4390844) was used and has been previously tested in our laboratory to ensure no cellular physiology changes. siRNA pool transfection was subsequently performed using Lipofectamine RNAiMax Transfection Reagent (ThermoFisher 13778150) according to manufacturer’s instructions and incubated for 12 hours post-transfection in serum free media followed by an additional 12 hours of recovery in SmBM supplemented media prior to RNA collection.

RNA Extraction and Quantitative PCR (qPCR) experiments
RNA extraction was performed using RNAEasy Mini kits (Qiagen 74106). 500ng of RNA per sample is then aliquoted for reverse transcription with High-Capacity RNA-to-cDNA Kit (Life Technologies 4388950). Quantitative PCR reactions were conducted with Taqman Universal Master Mix (44440048) and qPCR FAM probes for genes of interest (ThermoFisher) on a QuantStudio 6 Pro Real-Time PCR System. Relative transcript abundance was determined by the comparative Ct (Δ ΔCt) method, using the housekeeping gene UBC.

CEBPB, TEAD, H3K27ac ChIPseq
Approximately 1,000,000 HCASMCs were cross linked in 1% formaldehyde for 10 minutes and washed with PBS and replaced with hypotonic buffer (20 mM Hepes (pH 7.9), 10 mM KCl, 1 mM EDTA (pH 8) and 10% glycerol) and incubated on ice for 6 minutes. Cells were then sonicated using a Branson 250 Sonifier (using power setting 5, constant duty for 10 rounds of 30-second pulses) with confirmation of chromatin fragments at 250-400 base pairs. Lysates were then incubated overnight with 5 µg of anti-CEBPB (Santa Cruz sc-150), TEAD(Cell Signaling 12292S), or H3K27ac(Abcam 4729). Protein-DNA complexes were captured with Protein G agarose beads (Sigma 8104) and eluted in 1% SDS TE buffer at 65 °C. After reverse cross-linking followed by RNase A and proteinase K digestion, chromatin was purified using a QIAquick PCR purification kit (Qiagen 28706). ChIP DNA sequencing libraries were prepared using the Kappa HyperPrep (Roche 07962347001) and sequenced on a NovaSeq6000 with 150-base pair paired-end reads.

Proximity Ligation Assay
PLA was performed using manufacturer provided protocols for Sigma DuoLink In Situ Red Start Kit Mouse/Rabbit (DUO92101) with antibodies to TCF21 (Sigma AV33421, anti-rb), TEAD1 (sc-376113, anti-ms) and anti-mouse IgG (Vector I-2000-1).

Co-Immunoprecipitation (Co-IP)
Co-IP was performed using a myc-Tagged Tcf21 construct cloned into the PwPI vector. HEK293 cells were transfected with Tcf21-myc followed by protein extraction and nuclear isolation following manufacturer’s instructions (Active Motif Nuclear Complex Co-IP Kit 54001). 5% input control was kept as a positive control and samples were incubated with anti-rb TEAD1 (Cell Signaling 12292S) or anti-rb IgG (Abcam ab171870) as negative control. Co-IP was performed following instructions from the Active Motif Co-IP Kit, briefly, samples were incubated with 1:100 of TEAD1 antibody according to manufacturer recommendations and 5ug of anti-rb IgG (which is in excess compared to amount of TEAD1 antibody used). After overnight incubation, samples are washed and protein collected for western blotting. An anti-mouse HRP secondary antibody targeting the light chain was used for primary antibody detection as TEAD1 size is approximately 47kDA which would overlay with heavy chain antibody fragments.

Immunohistochemistry (IHC)
IHC was performed for TEAD1 and TCF21 in mouse (TgMyh11-CreERT2, ROSAtdT/+, ApoE−/−) aortic root sections fed 16 weeks of high fat diet to characterize their protein localization. Briefly, antigen retrieval was performed on sections after dissolving OCT in distilled water, and fixing in 4% paraformaldehyde for 5 minutes. After two PBS washes, antigen retrieval buffer (Biocare, DV2004) diluted in deionized water was preheated to 97–100°C in an Oster steamer. Sections were pre-treated with RNAscope hydrogen peroxide for 5 minutes twice, then transferred into antigen retrieval buffer for 6 minutes, ensuring the buffer remained at least 97°C. Slides were immediately moved to distilled water and washed twice by gentle dipping, air-dried, and a hydrophobic barrier was drawn around tissue with a hydrophobic pen. Section are incubated with Rodent Block M (Biocare, RBM961) for 30 minutes and washed in TBS. Primary anti-TEAD1 (Abcam 1335N33) and anti-TCF21 (Sigma, HPA013189) was diluted 1:100 in DaVinci Green Diluent (Biocare, PD900) and incubated overnight at 4°C. The next day, sections were washed in TBS, incubated with Rabbit-on-Rodent HRP Polymer for 30 minutes at room temperature, and washed in TBS. Vina Green Chromogen kit was prepared per manufacturer instructions, applied at 50 µL per section for 4 minutes, then slides were rinsed in deionized water twice, air-dried, and mounted with a xylene-based medium under coverslip for microscopy.

Dual luciferase assays
Enhancer/Promoter elements (SRF 2nd intron; chr15:73633709-73634109, BMP1 6th intron; chr8:22178981-22179700, Loxl1 5’ UTR; chr15:73926050-73926450) demonstrating overlapping between Tcf21 and Tead ChIP-seq sites were cloned into pWPI and transfected into A7r5 cells. A7r5 cells were seeded into 24 well plate (1.5 × 104 cells/well) in DMEM containing 10% FBS and incubated at 37 °C and 5% CO2 overnight. Cells were transfected with varying combinations of luciferase reporter plasmids (pLuc-MCS (empty), pLuc-enhancer, cDNAs (pWPI (empty), pWPI-TCF21, pWPI-TEAD1 and pWPI-MYOCD), and Renilla luciferase plasmid using Lipofectamine 3000 (Invitrogen, #L3000015). Six hours after transfection, the media was changed to fresh complete media. Relative luciferase activity (firefly/Renilla luciferase ratio) was measured by SpectraMax L luminometer (Molecular Devices) 24 h after transfection. All experiments were conducted in at least triplicate and normalized to the reporter plasmid after subtracting empty luciferase construct luminescence.

Statistics & Reproducibility
To identify differentially expressed genes in the scRNA-Seq data, we split our data by PseudoEarly (pseudotime bins 3-24) and PseudoLate (pseudotime bins 25-30) and employed the FindMarkers function with the DESeq2 wrapper and filter for genes with absolute value Log2FC > 0.15. For differential scATAC peak analysis, we applied the likelihood ratio test ‘LR’ from FindMarkers and used nCount_peaks as the latent.variable and min.pct as 0.001 on a per cell type cluster basis. For differential ChromVar analysis, we calculated row means and calculated the average difference in Z-score between conditions using the wilcox test from FindMarkers.

For Jensen-shannon divergence (JSD) calculation we compared the ChromVar motif deviation distribution for each motif across pseudotime by calculating the Kullback-Leibler (KL) divergence from each distribution to the midpoint distribution where:
(1) KL(p||m) = ∑ p·log₂(p/m) KL(q||m) = ∑ q·log₂(q/m)
Where p and q represent the two distributions being compared and m is the midpoint distribution calculated as the average of two normalized distributions
(2) m = (p + q)/2
Then the JSD was computed as the square root of the average of the two KL divergences.
(3) JSD(p||q) = √[(KL(p||m) + KL(q||m))/2]
We first randomly split the control dataset equally and calculated the JSD divergence to derive a control distribution for each motif. We then calculated JSD divergence between Tcf21-KO and control. To determine statistical significance, we constructed a null distribution using the control JSD measurements by performing 10,000 random sampling with replacement from the control values. One-sided empirical p-values were calculated for each TF using the formula p = (number of permuted values > observed value +1) / (total number of permutations + 1).

Fisher's exact test was used to assess the significance of overlaps between genomic regions. Adjusted P-values are corrected with FDR method of Benjamini-Hochberg and <0.05 were considered statistically significant. All data are presented as mean and error bars represent standard deviation (SD).

For dual luciferase experiments, the BMP1 and LOXL1 enhancers were tested in biological triplicates. For SRF 2nd intron, the test conditions were performed in quadruplicate, while control and MYOCD-only conditions were performed in duplicate. For SRF, the control and MYOCD-only conditions have been previously validated by Nagao et al.35 Each dual luciferase experiment was repeated three times independently with representative results shown.

No statistical method was used to predetermine sample sizes for in vitro and single-cell experiments. No data were excluded from the analyses. The experiments were not randomized. The Investigators were not blinded to allocation during experiments and outcome assessment.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In the paper “Vascular smooth muscle cell state trajectories mediate molecular mechanisms of coronary disease risk”, Li, Kundu et al. studied smooth muscle cell state transition measured by transcriptomics and chromatin accessibility in mouse aortic root across ApoE-/- induced atherosclerosis dense time points. The authors linked fibroblast and osteochondrogenic phenotype in smooth muscle cells to atherosclerotic states in transcriptomics and in special organization. The fibroblast-like cells (fibromyocytes, FMCs) had two subtypes, one of which was associated with immune activation and the other was associated with fibrous cap spanning the intimal plaque. The osteochondrogenic smooth muscle cells (chondrocytes, CMCs) had a distinct chromatin accessibility pattern. The authors further performed Waddington optimal transport to explore cell state transitions. There was an increase in FMC and CMC states at the atherosclerosis endpoint. The authors identified transcription factors enriched across cell state transitions. With aligned scRNA-seq and scATAC-seq data, the authors identified cell state specific transcription factor motifs from scATAC-seq data across pseudotime points. Gene regulatory network analysis in the co-embedded scRNA and scATAC dataset showed Tcf21 enrichment in SMC-FMC transition. The authors further performed knockout on Tcf21 and observed increase in SMC-3 cluster in early time points and decrease in FCM and CMC cell states at late atherosclerosis states. TCF21 was also found to interact with TEAD1 gene in vitro. An integrated analysis showed that FMC and CMC cell states have higher enrichment of CAD GWAS genes. Overall, the authors conducted a carefully designed time course experiment to identify critical cell state transitions, the underlying cellular process, and key genes that drive atherosclerosis in smooth muscle cells. The analysis is novel, and the underlying data will be a wonderful resource for the vascular biology field.\nQuestion/Suggestions:\n1. Was there heterogeneity in TdTomato expression in the cells analyzed? Would be good to show a feature plot for TdTomato expression to determine if there were cell-types that started with less expression of Myh11.\n2. The main figures have been rasterized, so many figures are unreadable. Figure 1 O-Q, for example, are completely unreadable.\n3. Discussion section would benefit from comparison between Waddington OT and other trajectory analysis methods? For example, CellRank is widely used, and would be helpful to know the nuances that motivated the use of Waddington OT.\n4. Current in situ hybridization is only show at 1 time point (Figure 1L,M,N). If this can be done at a few time points it would validate there is a migratory path.\nMinor comments:\n5. Is there spatial specific expression quantification for Fig 1L,M,N?\n6. The authors listed several transcription factors found through WOT analysis and TF enrichment in the paragraph from line 173 to 191. It will be good to have a summary describing how FMC and CMC regulate proliferation during atherosclerosis.\n7. How were the PseudoEarly and PseudoLate bins determined in scATAC-seq analysis? Since there is scATAC-seq data taken at each time point, could the authors perform enrichment analysis based on the cell state clusters within each time point? Would a random separation of pseudotime points show similar TF enrichment results? The random separation TF enrichment analysis could potentially show significance of the selected SMC-FMC and FMC-CMC transitions.\n8. Can the authors add Supp Fig 2 O slide-seq figure color bar?\n9. Is the cell population distribution in Figure 2J showing changes within the 2-to-3-week time? If so, the cell population fluctuates a lot within a 3-week time frame. Are the plots on the same scale across the weeks?\n10. How was the osteoblast progenitor score calculated?", "reviewer_2": "In this manuscript, the authors continue their longstanding strong work of using advanced multiome computational genomic tools with mouse and human genetics to elucidate the cell fate of SMC and their contribution to atherosclerosis progression. Here, the authors present a comprehensive single cell transcriptomic and epigenomic analysis of SMC populations and apply probabilistic Waddington Optimal transport (WOT) modeling across a timecourse of atherosclerosis progression (in a mouse ApoE KO model) to define dynamic SMC cell states and trajectories. By integrating scRNA-seq and scATAC-seq over time, with many more timepoints than most studies, they are able to infer a higher resolution TF regulon network to predict the top drivers of SMC phenotypic switching (SMC types, n=3) to fibromyocyte (FMC types, n=2) and to chondromyocytes (CMC). Furthermore, they use their previous model, Tcf21 KO, to gain further insight as to the changes in SMC transdifferentiation in this current work and to validate their in-silico TF perturbation computational predictions. Importantly, they integrated human CAD GWAS studies to associate changes in SMC subpopulations and cell states with disease risk in humans. Overall, in comparison to their previous work in Sharma, et al, ATVB, 2024, they now present a deeper genomic, transcriptional, and regulatory view rather than descriptive model of SMC phenotypic switching linked to human disease. This new work is technically rigorous and leverages multiome single cell sequencing data with advanced computational analysis to provide new regulatory insights into the relationship between SMC phenotypic switching and CAD.\nMajor Points:\n• The detailed time course, computational rigor and the use of Tcf21 KO data (an established CAD gene already well studied by this group) to probe the models are strengths. Yet the reliance alone on Tcf21 as the genetic probe has limitations and reduces the novelty given how much work this group has already done on Tcf21 KO mouse models and Tcf21 human genetics. And, despite the deep dive with multiple time points and integrated computational multi-omics approaches, the inferences are confusing in places. Based on human genetics and mouse genetic models prior work infers that TCF21 may serve a protective role in atherosclerosis. This serves as a grounding for interpretation of cell states, trajectories and their regulation and their function in disease. Yet, Tcf21 is found to regulate (increase) FMC cell state – inferred as protective - and also regulate (increase) CMCs - proposed to be harmful. Based on human genetics, the authors interpret the important disease related cell regulation as FMC increase and the inferred Tcf21-mediated CMC increase to be a downstream mass effect that is not important (not harmful) but there is no direct data to support that conclusion. Similarly, the authors provide evidence that FMC1 and FMC2 can be terminal cell states but that some FMCs can transition to more harmful CMC… yet it is completely unclear what might regulate the transition or lack of transition of FMC types to CMCs – which appears to be an important cellular functional event in disease. The authors should be more circumspect in their interpretation and interrogate more comprehensively these contradictions in their data. Ideally, they would interrogate separate data (even if limited and public) for an independent and opposing key regulator of SMC phenotype switching e.g., Klf4.\n• With respect to some subset of FMC1 and FMC2 populations being terminally differentiated based on data at 16 weeks of western diet (WD), more temporal data should be added especially as 16 weeks WD may not capture the complexity of later lesions and human disease. It is important to include a longer timepoint on WD, ~ 20 or 26 weeks, to determine if this holds true. Furthermore, this would also better elucidate the trajectory of CMC population with atherosclerosis progression and more advanced disease.\n• The authors use WOT algorithm that was initially developed by Eric Lander’s and Aviv Regev’s labs to improve upon already available methods for cell trajectory or lineage inferring across developing tissues. However, all these methods, including the one in this manuscript, lacks the cellular resolution to accurately infer lineage tracing, specially for rarer populations. This seems evident by the apparent lack of origin of FMC2 subcluster at the earlier timepoints (~ 5 Wk). While the Sankey representation of inferred transition of state in Fig 3D clearly shows the transition SMC2 -> FMC1 -> FMC2 or SMC2 -> FMC2 in the later timepoints, it seems to have a problem inferring transition state changes for FMC2 at 5 wk. A similar issue is encountered with CMC population in the Tcf21 KO at 12 and 16 wks. The authors must comment about the temporal and cellular resolution limits of this approach and strategies to possibly improve it (i.e. adding timepoints in between?).\n• In this context, it is also notable that the authors appear to use a more stringent transcript threshold than their prior work (Cheng et al., Nat Card Res. 2022 and Sharma et al., 2024) and this may miss rarer smaller clusters down the lineage lines. Indeed, In Cheng et al., Nat Card Res. 2022, the authors previously identified, albeit at a lower resolution, a smaller R-SMC cluster located in the basal plaque and primarily identified by their high Mmp3 expression, but distinct from SMC, FMC and CMC, and they hypothesized that these derived from FMC by pseudotime analysis with Slingshot. What happened to this SMC subcluster? Was it integrated into the new CMC cluster or was it not detected on this dataset?\n• The authors should show the actual percentage change in FMC2 population between control and Tcf21-KO at 16 weeks? From figure 4A, the overall change seems negligible when compared to the decrease in FMC1 and CMC. Is the overall increase in atherosclerosis in Tcf21-KO attributed mainly to the expansion of SMC3, given the lack of changes FMC2 population? Another mechanism that might play a role beside the changes in phenotypic switching or percentages of modulated SMC population in Tcf21-KO is their localization. Is there a significant change in the localization of FMC2, FMC1, CMC and particularly SMC3 population in control vs Tcf21 KO?\nMinor Points\n• Does the osteoblast, senescence, EMT, apoptosis, angiogenesis or efferocytosis cell state of SMC3 changes between control and Tcf21?\n• Can the authors provide the actual percentages of the 6 different SMC subpopulations at all timepoints in Fig 2C?\n• Given that control TCF21 + IgG showed some unspecific background PLA signal, it would be worth it to show the second control (TEAD1 + IgG) and quantifications of the background signal in the negative controls as well as in TCF21 + TEAD1. It seems that there is also significant interaction between Tcf21p and Tead1p in the cytosol as well. Given this, it would be better to include more information regarding the images. Are they epifluorescence or laser confocal images? If epifluorescence, are they deconvolved and how? Are they max projection or single slices (mentioning objective + NA and pinhole size, or optical slide depth)?\n• PLA assays detects proximity of proteins up to max of 35~ 40 nm, well above the Rmin of the predicted dimer for TcF21-Tead1 interaction, thus it is still not direct proof of binding/interaction between these protein. Can/did the authors try coIP? (acknowledging coIP from nuclear lysates can be difficult). This would provide better evidence of physical interaction between these two protein or coIP their intermediate if the interaction is not direct.\n• It seems that other TF besides the KLF family and AP1 factors, such as Smad2 and Smad3, also show a bimodal open motif distribution. Could they be playing a dual role in early SMC to FMC transition as well?\n• Figure 4B: label the top Sankey plot as “Control” and bottom as “Tcf21-KO” in the graph, not only mentioning it in the figure legend", "reviewer_3": "This is a detailed and comprehensive single-cell and epigenomic study on VSMC transitions in atherosclerosis in the context of CAD. It aimed to study trajectories of SMC state transitions, identify novel TFs and gene regulatory networks contributing to these transitions, and how the Tcf21 gene regulates these processes. Key findings included the identification of six SMC states, with the FMC-1 and FMC-2 states contributing to CMC. Tcf21 emerged as the top regulator for the SMC to FMC transition, along with other TFs, ie. Runx, Klf4, Smad3, HIF etc. Tcf21 KO resulted in reduced FMC and CMCs and increased SMCs, and TEAD1 was found to interact with Tcf21 to modulate TGFB and ECM-related profibrotic pathways.\nStrengths:\nThe authors utilized advanced multi-omic tools such as single-cell multi-omics (scRNA+scATAC) and technical computational modeling (ie. WOT, CellOracle, Pando..etc). They conducted lineage traced atherosclerosis mouse model experiments to assess the transition of SMC switching across 4 time points and the role of Tcf21. They assessed translational relevance by integrating with human genetic data (GWAS, scDRS, S-PrediXcan). They discovered a novel mechanism in the Tcf21-TEAD1 interaction and importantly identified FMC2 as a key CAD-risk state.\nWeaknesses:\n1. Three key findings were that 1) Tcf21 plays a strong role as a driver of contractile SMCs transitioning to FMC-1 and that 2) FMC-2 was identified as the key FMC cluster associated with CAD risk. 3) Knocking out Tcf21 led to decreased FMC1 and FMC2 and CMC proportions and increased SMC2 and SMC3 at 12 weeks, and further decrease of FMC1 but not FMC2 at 16 weeks. Since this was the case, did the authors consider another major candidate regulator/TF associated with driving SMCs to FMC2 (ie. TEAD1?) since this seems to be a very important natural point to the study’s narrative?\n2. The authors provide important in situ transcript validation using RNAScope. However, additional protein-level validation of FMC markers would strengthen the evidence beyond the RNA landscape towards functional cell phenotypes.\n3. The evidence of the role of TEAD1 from ChIP-seq and reporter assays support TEAD1 as a TCF21 partner. Functional perturbation of TEAD1 would significantly strengthen the mechanistic conclusions.\n4. If the conclusion is that Tcf21 promotes the protective FMC1 and competes with TEAD1 which promotes the CAD-risk-enriched state FMC2, it may be be important to know whether TEAD1 acts as an independent disease promoting factor or modulates risk through its interaction with TCF21.\n5. The use of CellOracle provided a rich set of candidate TFs for driving SMC phenotypic switching. With the exception of Tcf21, some (even limited) validation of additional TFs could strengthen conclusions made.\n6. Some functional validation in primary SMCs to identify and distinguish FMC1 and FMC2 would strengthen the story inferred from the multi-omics data.\nMinor corrections:\nFigure 6 legend: repeat of “and” in title."}, "author_response": "Response to Reviewer Comments:\nReviewer #1 (Remarks to the Author)\nIn the paper “Vascular smooth muscle cell state trajectories mediate molecular mechanisms of coronary disease risk”, Li, Kundu et al. studied smooth muscle cell state transition measured by transcriptomics and chromatin accessibility in mouse aortic root across ApoE-/- induced atherosclerosis dense time points. The authors linked fibroblast and osteochondrogenic phenotype in smooth muscle cells to atherosclerotic states in transcriptomics and in special organization. The fibroblast-like cells (fibromyocytes, FMCs) had two subtypes, one of which was associated with immune activation and the other was associated with fibrous cap spanning the intimal plaque. The osteochondrogenic smooth muscle cells (chondrocytes, CMCs) had a distinct chromatin accessibility pattern. The authors further performed Waddington optimal transport to explore cell state transitions. There was an increase in FMC and CMC states at the atherosclerosis endpoint. The authors identified transcription factors enriched across cell state transitions. With aligned scRNA-seq and scATAC-seq data, the authors identified cell state specific transcription factor motifs from scATAC-seq data across pseudotime points. Gene regulatory network analysis in the co-embedded scRNA and scATAC dataset showed Tcf21 enrichment in SMC-FMC transition. The authors further performed knockout on Tcf21 and observed increase in SMC-3 cluster in early time points and decrease in FCM and CMC cell states at late atherosclerosis states. TCF21 was also found to interact with TEAD1 gene in vitro. An integrated analysis showed that FMC and CMC cell states have higher enrichment of CAD GWAS genes. Overall, the authors conducted a carefully designed time course experiment to identify critical cell state transitions, the underlying cellular process, and key genes that drive atherosclerosis in smooth muscle cells. The analysis is novel, and the underlying data will be a wonderful resource for the vascular biology field.\nQuestion/Suggestions:\n1. Was there heterogeneity in TdTomato expression in the cells analyzed? Would be good to show a feature plot for TdTomato expression to determine if there were cell-types that started with less expression of Myh11.\nThank you, this is an important question. We show in Supp Fig 1 that there is a bimodal distribution of tdTomato expression which likely represents ambient RNA or doublets which incorporated tdTomato expressing cells (as discussed in Sharma et al.1). To address your question more directly, we have now included an additional panel (Supp Fig 1H) to show the distribution of tdTomato at baseline to demonstrate the specificity of Myh11 driven tdTomato distribution. It is important to remember that recombination efficiency is quite high, ~95%, so cells expressing even low levels of Myh11 will activate the high-level expression of the tdTomato reporter, which does not really reflect endogenous Myh11 expression level.\nResults, line 119:\nAll of the phenotype clusters represented cells of the SMC lineage, as indicated by lineage tracing with the Myh11-Cre transgene (Fig. 1H). It is important to note that recombination efficiency is quite high, ~95%, so cells expressing even low levels of Myh11 will activate the high-level expression of the tdTomato reporter, which does not really reflect endogenous Myh11 expression level that might reflect some heterogeneity in the SMC cell population.\n2. The main figures have been rasterized, so many figures are unreadable. Figure 1 O-Q, for example, are completely unreadable.\nThank you for pointing this out. We apologize for this mistake and have re-rendered all of the images to fit the quality standards for the submission.\n3. Discussion section would benefit from comparison between Waddington OT and other trajectory analysis methods? For example, CellRank is widely used, and would be helpful to know the nuances that motivated the use of Waddington OT.\nAt the time of our data development and temporal analyses, Waddington OT was a well-utilized trajectory method that incorporated experimental time as a key feature in addition to modeling of variable cell growth rate/death (cell cycle and apoptosis) signatures derived from transcriptomic profiles rather than just comparing cell-cell expression similarity. We further selected time points based on clear real time phenotypic stages that occur in the data, which is another key assumption in this method. This allowed for a unified framework for further time-course based analyses.\nSimilar real-time based methods such as those incorporated in CellRank22 had been just described in preprint and so we again, decided to apply well-utilized orthogonal methods from the original CellRank3, such as its pseudotime kernel, to generate a probabilistic model, which appears to follow similar patterns as our optimal transport approach. This provided the added advantage of placing critical TFs across time in a sequential pattern to follow the cascade of transcription factor activation.\nTaken together, these methods expand upon initial studies which only used pseudotime to infer disease trajectories in this murine atherosclerosis model4,5. We have now incorporated these points into the discussion section with the following text.\nDiscussion, line 600:\nComputational algorithms used in these studies were identified to be those most aligned and broadly tested with the genomic and genetic data types provided in our timecourse study. In particular, novel trajectory methods that become available may continue to offer useful and more nuanced details regarding the complex cell state changes that the SMC undergo during the disease process and we make our data available for such analyses. While the WOT approach was specifically designed for timecourse scRNAseq data, the CellRank22 algorithm is now considered to potentially yield more accurate and robust results by integrating additional biological information, such as metabolic labeling, and can use optimal transport as one component of its broader multivariate analysis. We did in fact utilize orthogonal methods from the original CellRank3, such as its pseudotime kernel, to generate a probabilistic model, which appeared to follow similar patterns as our WOT approach to generate interpretable summary models.\n4. Current in situ hybridization is only show at 1 time point (Figure 1L,M,N). If this can be done at a few time points it would validate there is a migratory path.\nWe agree that this is an important point and performed additional RNAScope to follow-up on this recommendation. We performed RNAScope at baseline and 7 weeks of high fat diet (HFD) to resolve the localization of top FMC-1 and FMC-2 markers (Fbln1, Fbln2, respectively) to show that at baseline, Fbln1 is seen at low level in the media and Fbln2 is not. Then, at 7 weeks of HFD, there is enrichment of Fbln1 in the media and at the cap (at lower levels) while Fbln2 is more restricted to the intima when compared with FMC-1. It is important to note that the comparisons of Fbln1 and 2 in these slides represent serial sections, therefore, highlighting the beginning differences in expression pattern between these marker genes. This aligns closely with our single cell data which shows low level FMC-1 but no FMC-2 marker expression at baseline while FMC-1 and FMC-2 expression ramps up by week 7 of HFD. These results are shown in the revised manuscript (Fig. 2B, C, S-Fig. 3N, O).\nFurther, we expanded our cell transition localization analysis with high quality Xenium 5k geneset data for the mouse aorta and its branches. With transfer of single cell RNA sequencing labels to the spatial clusters, we demonstrated that FMC-1 signals are enriched in the media and also sporadically localized in the cap by 16 weeks of HFD. Subsequently high levels of FMC-2 were identified in the cap and spanning the intima to the CMCs at the base of the plaque, further validating our inference based on selected RNAScope markers. The additional Xenium data are provided (Fig. 1P, Q, S, S-Fig. 2A) and related text is provided in the Results section of the manuscript.\nResults, line 167:\nTo examine the translocation of FMC from media to cap to plaque, we performed additional studies to map temporal gene expression of FMC-1 and FMC-2 marker genes Fbln1 and Fbln2 across serial sections to highlight their distinct expression patterns. RNAScope performed at baseline showed that FMC-1 marker Fbln1 expression was seen in the media while expression of FMC-2 maker Fbln2 was not (S-Fig. 3N). Then, at 7 weeks of high fat diet (HFD), there was enrichment of Fbln1 in the media and at lower levels at the cap, while Fbln2 was more intimally restricted with lower medial expression compared with FMC-1. (S-Fig. 3N, O). These findings aligned with our single cell data which showed low level FMC-1 but no FMC-2 marker expression at baseline while FMC-1 and FMC-2 expression increased after week 7 of HFD (Fig. 2A, B, C). We further investigated the lesion cellular anatomy at 16 weeks HFD with high quality Xenium 5k geneset data for the mouse aorta. With transfer of single cell RNA sequencing labels to the spatial clusters, we demonstrated that FMC-1 are enriched in the media and also localized in the cap after 16 weeks HFD. Subsequently high levels of FMC-2 were identified in the cap and spanning the intima to the CMCs at the base of the plaque, further validating our inference based on selected RNAScope markers (Fig. 1O-T, S-Fig. 2A).\nFor our spatial analysis, we have incorporated the new methods in the Materials and Methods section of the manuscript.\nMaterials and Methods, line 811:\nXenium Slide Processing and analysis\nFresh frozen tissue sections were prepared following the 10x Genomics Xenium In Situ for Fresh Frozen Tissues Protocol (CG000579 Rev F). Briefly, tissues arranged in order of aortic root, ascending aorta, brachiocephalic artery with right subclavian branch, per-diaphragmatic descending, and abdominal descending aorta from a control (TgMyh11-CreERT2, Tcf21+/+, ROSAtdT/+, ApoE−/−) mouse were embedded in OCT compound immediately after dissection. Embedded tissue blocks were frozen on dry ice and immediately stored at -80C prior to sectioning.\nSections of 10m thickness were cut using a Leica CM1860 cryostat. Sections were mounted within the sample area (10.45 x 22.45mm) of Xenium slides (PN-3000941) without overlap and stored at -80C until ready for use. All subsequent Xenium steps were processed at the Stanford Functional Genomics Core. Tissue morphology and quality was assessed using Hematoxylin and Eosin (H&E) and DAPI staining prior to downstream Xenium In Situ Gene Expression assays.\nProcessed Xenium datasets were imported into Seurat V5 for normalization and integration. Feature-barcode matrices were normalized using the SCTransform workflow with spatial coordinates retained for downstream deconvolution. To perform single-cell reference label transfer to the Xenium data, canonical correlation analysis (CCA) was implemented following strategies from the Seurat tutorial “Analysis, visualization, and integration of spatial datasets with Seurat” (https://satijalab.org/seurat/articles/spatial_vignette.html). RCTD (Robust Cell Type Decomposition) was applied to further deconvolve spot-level spatial data based on input cell reference.\nMinor comments:\n5. Is there spatial specific expression quantification for Fig 1L,M,N?\nGiven the specificity of the cellular localization and the low background with the RNAscope hybridization, and the strong correlation with the single cell data we did not feel that quantitation of the RNAscope was necessary. We have since also added Xenium data to further validate our observations with selected RNAScope markers and also as a public resource for future studies. We feel that quantitation of cell numbers for the different transition phenotypes is best reflected by the cluster cell numbers in the scRNAseq data.\n6. The authors listed several transcription factors found through WOT analysis and TF enrichment in the paragraph from line 173 to 191. It will be good to have a summary describing how FMC and CMC regulate proliferation during atherosclerosis.\nWe thank the Reviewer for this suggestion. We have now included an additional supplementary table (S-Table 7 ) to describe the function of TFs related to optimal transport/cell oracle nominated TFs, including regarding proliferation and other cell state changes.\nResults, line 204:\nAcross the enriched TFs found in our model, several patterns stand out. For example, TFs (Runx1/2, Zbtb7c, Zeb1/2, Smad3, Stat3, Hes1, Nfkb1, Cebpb) amongst this list have been readily associated in literature with AP-1, Klfs and the SMC regulatory Srf-myocardin complex. This suggests that pioneer factors and immediate-early AP-1 activation constitute a core early event that then allows for recruitment of transcription factors at specific enhancers to drive key cell state transition steps in conjunction with Srf-myocardin. These core TFs then organize around central cellular processes such as TGFB signaling, hypoxia signaling, inflammatory response and osteochondrogenesis as summarized in S-Table 76-10.\n7. How were the PseudoEarly and PseudoLate bins determined in scATAC-seq analysis? Since there is scATAC-seq data taken at each time point, could the authors perform enrichment analysis based on the cell state clusters within each time point? Would a random separation of pseudotime points show similar TF enrichment results? The random separation TF enrichment analysis could potentially show significance of the selected SMC-FMC and FMC-CMC transitions.\nThank you for this comment. The PseudoEarly and PseudoLate bins are based on the distribution of cell state clusters where the majority of phenotypically transitioning cluster cells are covered by the PseudoLate bins (bins 25-30). This was demonstrated by the distribution of cell clusters in Figure 3A. However, this was shown only by a small graphic above the heatmap, and we now have highlighted the timepoint and cluster designations with expanded bar graphs above the ATAC pathway clustering in Figure 3A. Further, we have made additional clarifications in the text and figure legends to highlight the importance of these expanded graphics Although there are multiple timepoints for the scATAC data, because all clusters are represented in the later timepoints, TF enrichment analyses on a timepoint basis may not accurately delineate transitioning TF motif enrichment. Therefore, we took a pseudotime approach, which not only showed a strong correlation with the appearance of cell state clusters over time, but also helped unify our analysis across all the timepoints.\nWe further investigated the question regarding random separation of pseudotime bins to assess how TF enrichment may shift. We utilized the pseudotime aggregated ChromVar scores and randomly shuffled the 30 pseudotime bins 1000 times by bootstrapping and created an average pseudotime bin distribution. Here, we clearly observed a significant visual shift in the distribution of TF motif deviations from a discrete to a uniform pattern upon shuffling of pseudotime, see embedded Reviewer figure below (Reviewer Figure 1).\nResults, line 248:\nTo characterize the epigenetic processes that mediate the noted transcriptional effects, we investigated chromatin accessibility in the transitioning cells with scATACseq data (Fig. 1C, D). We observed high specificity of chromatin accessibility within each label transferred tdTomato lineage traced cluster and visualized top specific peaks along pseudotime (Fig. 3A). Using GREAT for functional enrichment11, we identified specific cellular processes including smooth muscle related processes in the SMC, “cell migration”, “inflammation” and “response to platelet derived growth factor” in the FMC, and “ossification” and “chondrocyte development” in the CMC (Fig. 3A, S-Table 8).\nReviewer Figure 1. (Left) The left heatmap visualizes aggregated ChromVar scores randomly shuffled with bootstrapping. (Right) Heatmap of organized motif deviations (ChromVar scores) across pseudotime.\n8. Can the authors add Supp Fig 2 O slide-seq figure color bar?\nYes, we have now included the color legend for slide-seq, now in updated S-Fig. 2L.\n9. Is the cell population distribution in Figure 2J showing changes within the 2-to-3-week time? If so, the cell population fluctuates a lot within a 3-week time frame. Are the plots on the same scale across the weeks?\nThe cell populations shown in the panels in Figure 2J represent the pseudotime distribution of cells at each timepoint along the same pseudotime scale rather than specifically the changes in cell population distribution. It shows that, at each timepoint, there is a distribution rather than discrete set of cells which only exist at that timepoint, thus representing our nominated cell states along the pseudotime continuum. Furthermore, there is a trend towards greater proportion of phenotypically transitioning cells with later pseudotime designation as we progress with weeks of high fat diet. These panels depict the expected distribution for all cells along pseudotime on the y-axes, with lower distributions that gradually trend toward greater pseudotime with high fat diet.\n10. How was the osteoblast progenitor score calculated?\nThe osteoblast progenitor score is obtained from https://doi.org/10.1073/pnas.2019152118 Haseeb et al. referenced as citation 29. We have now highlighted this point more clearly in the text as follows.\nResults, line 158:\nConsistent with this possibility, FMC-2 were also enriched for osteoblast progenitor gene signatures (Fig. 1M, S-Fig. 1L, S-Fig. 4A). Osteoblast progenitors have been characterized in bone development, and this module score12 determined with scRNAseq data characterizing these cells on their path toward osteochondroprogenitors.\nReviewer #2 (Remarks to the Author):\nIn this manuscript, the authors continue their longstanding strong work of using advanced multiome computational genomic tools with mouse and human genetics to elucidate the cell fate of SMC and their contribution to atherosclerosis progression. Here, the authors present a comprehensive single cell transcriptomic and epigenomic analysis of SMC populations and apply probabilistic Waddington Optimal transport (WOT) modeling across a timecourse of atherosclerosis progression (in a mouse ApoE KO model) to define dynamic SMC cell states and trajectories. By integrating scRNA-seq and scATAC-seq over time, with many more timepoints than most studies, they are able to infer a higher resolution TF regulon network to predict the top drivers of SMC phenotypic switching (SMC types, n=3) to fibromyocyte (FMC types, n=2) and to chondromyocytes (CMC). Furthermore, they use their previous model, Tcf21 KO, to gain further insight as to the changes in SMC transdifferentiation in this current work and to validate their in-silico TF perturbation computational predictions. Importantly, they integrated human CAD GWAS studies to associate changes in SMC subpopulations and cell states with disease risk in humans. Overall, in comparison to their previous work in Sharma, et al, ATVB, 2024, they now present a deeper genomic, transcriptional, and regulatory view rather than descriptive model of SMC phenotypic switching linked to human disease. This new work is technically rigorous and leverages multiome single cell sequencing data with advanced computational analysis to provide new regulatory insights into the relationship between SMC phenotypic switching and CAD.\nWe thank the Reviewer for these laudatory comments.\nMajor Points:\n• The detailed time course, computational rigor and the use of Tcf21 KO data (an established CAD gene already well studied by this group) to probe the models are strengths. Yet the reliance alone on Tcf21 as the genetic probe has limitations and reduces the novelty given how much work this group has already done on Tcf21 KO mouse models and Tcf21 human genetics. And, despite the deep dive with multiple time points and integrated computational multi-omics approaches, the inferences are confusing in places.\nBased on human genetics and mouse genetic models prior work infers that TCF21 may serve a protective role in atherosclerosis. This serves as a grounding for interpretation of cell states, trajectories and their regulation and their function in disease. Yet, Tcf21 is found to regulate (increase) FMC cell state – inferred as protective - and also regulate (increase) CMCs - proposed to be harmful. Based on human genetics, the authors interpret the important disease related cell regulation as FMC increase and the inferred Tcf21-mediated CMC increase to be a downstream mass effect that is not important (not harmful) but there is no direct data to support that conclusion.\nWe believe that these studies have provided considerable novel information regarding the specific role of TCF21 toward regulation of SMC phenotype transitions and disease risk. The Reviewer is correct that TCF21 is now a highly validated CAD GWAS causal gene and represents an archetypal contributor to SMC phenotypic transitions that by inference from human genetics represent a protective mechanism toward disease risk. However, our original single cell RNA sequencing studies published in Wirka et al.13 employed early (v2) 10x chemistry and thus provided a limited characterization of the cellular and molecular effects of Tcf21 in the disease setting, and was not paired with single cell epigenetic profiling. Importantly, we did not have sufficient cells and identified genes in that study to fully characterize the chondromyocyte (CMC) phenotype. We thus included Tcf21 knockout in these comprehensive temporal transcriptomic and epigenomic studies to place it in the SMC transitional regulatory framework, and link TCF21 to other CAD associated factors that promote or inhibit risk in this context. Through these studies we have shown that Tcf21 inhibits transition of secondary heart field cells to a differentiated SMC phenotype we have labeled SMC-3, placed Tcf21 function in the FMC1 subpopulation of the fibromyocyte cell phenotype, and characterized the disease regulon of Tcf21 target genes. These findings have significantly expanded our understanding of its cellular and molecular functions, but as discussed below, also identified some confusing features regarding its role in the transition from FMC to CMC transition phenotype.\nThe comments by the Reviewer have also prompted us to reconsider our previous studies, and to incorporate the Klf4 KO analyses into our considerations regarding SMC disease causality and the FMC2 and CMC transition phenotypes. Previous studies in our initially suggested that transition to the CMC phenotype might be the mechanism of disease risk for SMC lineage cells. Mouse knockout studies for the CAD protective orthologous human SMAD3 gene revealed increased CMC transition, vascular calcification, and disease burden4. Knockout of the CAD disease promoting PDGFD mouse ortholog revealed decreased CMC numbers and decreased vascular calcification, and while there was no change in plaque burden there were a number of other cellular features consistent with plaque destabilization, including increased inflammatory profile and macrophage recruitment. Knockout of the mouse Ahr gene further substantiated the correlation of CMC transition phenotype to vascular calcification and disease burden. These findings, in conjunction with extensive clinical and epidemiology data suggesting that vascular calcification is a strong predictor of CAD risk, and likely causal contributor to disease pathology14-16, our initial hypothesis was that the CMC phenotype and related vascular calcification is the mechanism by which CAD associated primarily SMC lineage genes promote disease risk. It was unclear at that time whether FMC represented a transition phenotype or possibly a stable phenotype that might have a protective effect, possibly retarding transition to the disease promoting CMC transition. Thus our expectation going into these studies was that knockout of the protective Tcf21 gene would lead to an increase in CMC and that this cellular phenotype would show evidence of increased risk.\nIn an effort to assign genetic disease risk to SMC state transitions we conducted studies reported here employing single cell disease relevance score (scDRS) analyses to gain insights linking SMC transition phenotype and genetic disease risk and these analyses pointed to the FMC-2 phenotype. It was a surprising but also important finding that these methods did not assign genetic risk to the CMC, but this observation in conjunction with the strong evidence that protective gene Tcf21 promotes transition to the CMC phenotype led us to propose a unifying model that genetic risk predominantly resides in the FMC-2 phenotype, and not the CMC phenotype, with the corollary that Tcf21 promotes an increase in non-disease promoting CMC through a simple “mass action” effect.\nThese Reviewer comments have led us to re-examine these data and to try and better integrate the directionality from genetic risk assessment with disease lesion features and risk prediction as linked to single cell gene expression and phenotype.\nAs noted by the reviewer we do not have evidence that convincingly eliminates the possibility that CMC contribute CAD risk, and this observation has led us to examine the correlation between SMC phenotype and the contribution of the extensively validated gene Klf4 to disease risk. To accomplish this, we did label transfer of our atherosclerosis time-course trajectory to the Klf4-KO data17, and found that knockout of this gene produced a prominent decrease in CMC number (20% vs 52%) that correlated with the previously described decrease in disease burden, fibrous cap contribution, and the identified osteochondrogenic cluster 7 transcriptional phenotype. These data in conjunction with our previous correlation of disease features with causal gene functions for Smad3, Pdgfd and Ahr has provided a compelling case for us to re-consider our model, focusing on the statistical and bioinformatic analyses that are in line with disease lesion features, and suggests a role for the CMC phenotype toward clinical disease risk.\nWe consider our finding of FMC-2 disease causality to be robust, and the analysis is supported by work published by Ord et al.5 using linkage disequilibrium analysis to show cells expressing Vcam1, essentially FMC-1 and FMC-2 combined, are the primary contributors to disease risk for the SMC lineage. Further, our scDRS analysis of data from Pan et al.18 also finds similar FMC-2 enriched genetic risk (S-Fig. 4A-C). Importantly, there are a number of reasons why CMC might not show risk with the scDRS algorithm. An important point is that there are actually few osteochondrogenic genes that reside in CAD risk loci, and this is surprising given the extensive list of chondrogenic development genes that are expressed in these cells. Perhaps allelic variation arising in these loci is lost through purifying selection. Interestingly, we did find greater differences in gene expression for the individual CMC genes that were included in the scDRS analysis, but the number of genes was not sufficient to provide a statistical association. For these reasons, we now take a broader view of the relationship of SMC transitions to disease risk, with both the FMC-2 and CMC considered as likely contributors to risk.\nCombining all available data, we now propose the hypothesis that FMC-2 provide a primarily protective role while CMC likely function to promote vascular calcification and disease risk.\nThus in the revised manuscript we continue to support the risk analyses that identify causal evidence for a risk-determining protective effect of FMC-2 phenotype cells. New to the revised manuscript, we propose a heuristic model that assigns CAD risk to the pro-calcification CMC phenotype cells, based on correlation of this cell state with disease burden and direction of human genetic risk. We point to a need for future research that identifies a greater number of CAD causal genes, especially those that are in keeping with the osteochondrogenic CMC molecular processes, that will allow more sophisticated assessment of risk ascertainment as focused on the transition to this SMC lineage phenotype. Further, we no longer present the mass action hypothesis to explain the TCF21 generation of CMC that do not contribute to disease risk but highlight the need to study the functional role of the CMC phenotype in general, and especially Tcf21-KO lineage CMC, i.e., how they may differ from non-Tcf21-KO lineage CMC in the context of disease risk. These points are now articulated in the Discussion with the following text:\nDiscussion, line 570:\nWe and others have used CAD GWAS findings along with gene expression or chromatin accessibility data to show that much of the risk for CAD resides in the coronary vascular SMC lineage5,19-21. Our high-resolution dataset extends upon previous observations and identifies FMC-2 as the cell state harboring significant CAD risk, based on disease gene expression integration with scDRS and PrediXcan, and suggests that the FMC-2 cell state is a protective mechanism toward human CAD risk (Fig. 6A). In situ RNA staining localizes the FMC-2 population to both the fibrous cap (Notch3, Thbs1) and the intimal plaque (Fbln2, Ttc9) confirmed with label transfer using Xenium spatial transcriptomics. Further molecular analysis finds FMC-2 at the juncture of critical gene modules for senescence and apoptosis while expressing numerous CAD associated genes that we have linked to atherosclerosis, including ZEB2 and SMAD34,22.Using this scDRS score leverages the power of GWAS genes to identify likely causal cell state determining genes whose functions are critical to prime cells towards phenotypic transition.\nThe scDRS analysis did not identify CAD risk in the CMC phenotype cells. This is surprising, given that mouse knockout disease models of orthologs of human CAD associated genes SMAD3 and PDGFD4,23, and other genes not yet linked to CAD by GWAS, KLF4 and AHR17,24, have demonstrated a significant correlation between CMC number, plaque burden, and vascular calcification. This is true for both disease promoting and disease protective gene functions, as determined by effect allele identification and genetics of gene expression data. It is possible that there is a dearth of informative allelic variation in CAD loci that determine the CMC phenotype. We did observe that aggregate expression of CAD disease risk genes in the CMC cell state showed a higher relative average expression of CAD GWAS risk genes, but this did not provide a statistically significant result. This observation suggests that a bias against regulatory variation at the CMC determining gene loci resulting in a low number of informative genes may explain the lack of an scDRS finding. Also, an important consideration is that FMC2 are the most similar to CMC in terms of gene expression phenotype, and the risk identified in these cells may not be protective but may in fact promote disease risk through directing differentiation to the CMC phenotype. eQTL based directionality inference is limited by the number of genes that have statistically significant eQTL links and this can bias results. Finally, the observation that Tcf21 promotes transition to the CMC phenotype is difficult to reconcile with these considerations but may simply reflect that the protective effect of increased FMC number outweighs the disease promoting effect of increased CMC.\nDiscussion, line 626:\nFuture studies are needed to confirm the functions of nominated novel genes implicated in the cell state changes characterized through these studies. Only through study of the larger causal gene regulatory network will we be able to understand which aspects of the complex cellular phenotypic changes, migratory behaviors and cell-cell interactions that are responsible for the risk that is modulated by the SMC lineage. Specifically, studies are needed to characterize in greater detail the FMC-2 phenotype genes that regulate CAD risk, what is the nature of CMC disease risk, and how can disease protective genes like TCF21 promote CMC formation without increasing risk.\nThe required expansive causal CAD gene characterization can only come from large-scale in vitro CRISPR screening with highly relevant disease cellular models, and similar screens conducted in disease model mice which will provide in vivo disease transcriptomic phenotype information regarding the function of SMC genes that both regulate SMC phenotype and CAD risk.\nSimilarly, the authors provide evidence that FMC1 and FMC2 can be terminal cell states but that some FMCs can transition to more harmful CMC… yet it is completely unclear what might regulate the transition or lack of transition of FMC types to CMCs – which appears to be an important cellular functional event in disease. The authors should be more circumspect in their interpretation and interrogate more comprehensively these contradictions in their data. Ideally, they would interrogate separate data (even if limited and public) for an independent and opposing key regulator of SMC phenotype switching e.g., Klf4.\nWe agree with the Reviewer that it is critical to identify genes that promote and inhibit the SMC-FMC and FMC-CMC transitions, and the current study was performed to further address this specific need. With work reported in this manuscript we have taken the first step, to identify and characterize in detail the cellular transition phenotypes and the genes and regulons that are predicted to regulate the critical cell state changes. Perhaps most important is the gene regulatory network that is functional in the FMC2 population of cells, which are the closest in phenotype to the CMC and likely are the “master regulators” of whether these cells transition to disease risk CMC or remain in a protective FMC2 phenotype. We have already shown that CAD protective gene Zeb2 promotes transition to the FMC phenotype, and possibly more importantly disrupts transition to the CMC phenotype22. Other genes that are identified with WOT analyses described in this manuscript, in addition to TCF21, include Ahr, Runx2, and Sox9, Gas7 and an inflammatory module that includes Irf5 and Nfkb, all of which are also likely regulators of the FMC2-CMC transition. Further studies are needed to assess these factors as modulators of the FMC2-CMC transition. We have now established methodology for in vivo Perturb-seq studies and will be able to investigate changes in transcriptional phenotype of SMC lineage cells in the disease setting with CRISPRi mediated gene knockdown and individual cell whole genome scRNAseq. These studies should allow us to identify which of these genes promote the FMC-CMC transition.\nAlso, as noted above we integrated scRNAseq data from the Owens lab characterizing the cellular and molecular features associated with Klf4 expression in the SMC lineage in the disease setting. Comparison of wildtype and Klf4 KO clusters showed a significant decrease in FMC-2 and CMC cell numbers with Klf4 deletion. This group of cells corresponds with the author’s Phactr1 staining which showed SMC media and fibrous cap enrichment, where we have localized this cell phenotype. Also, these findings are in line with our CellOracle predictions (Fig 3I) where Klf4 was predicted to be a lead gene promoting FMC to CMC transition. They also demonstrate the importance of the FMC2 population whereby decreased cells of this phenotype then results in decreased CMC, especially in the context of Klf4 knockout.\nWith respect to some subset of FMC1 and FMC2 populations being terminally differentiated based on data at 16 weeks of western diet (WD), more temporal data should be added especially as 16 weeks WD may not capture the complexity of later lesions and human disease. It is important to include a longer timepoint on WD, ~ 20 or 26 weeks, to determine if this holds true. Furthermore, this would also better elucidate the trajectory of CMC population with atherosclerosis progression and more advanced disease.\nThank you for this suggestion. We have used this opportunity to analyze prior data from Pan et al18 which utilized 22 week and 26 week high fat diet (HFD) exposure on an ApoE KO and Ldlr KO background, respectively. The sequencing chemistry from this legacy 10x data set is V2, but we were able to robustly identify the FMC-1 and FMC-2 populations from this dataset with label transfer, with cell state overlap visualized by a confusion matrix (Reviewer Fig. 2, new S-Fig. 4A-F, see below). Because this dataset excluded the aortic root, we did not observe the SMC-3 population in this analysis. Similar to the Pan et al. data, we further performed label transfer of our cell states onto the lineage traced Klf4 dataset from Alencar et al.17 which studied cells harvested from 18 weeks HFD. We similarly identified distinct FMC-1, FMC-2, and CMC populations with good correlation to our time-course data.\nWe apologize for using the presumptuous label “terminal phenotype” to describe the FMC-1 and FMC-2 clusters. In the revised manuscript we now use the term “sustained” to better reflect that these are cellular states that arise early and exist across extended timepoints for cells exposed to high fat diet, and also acknowledge that we have not proven that they are truly a terminal state that is resistant to further phenotypic changes. However, as recommended by the reviewer, we were able to consistently observe the presence of these cell states even with longer duration of high fat diet such as 18 and 22 weeks for ApoE and 26 weeks for Ldlr KO atherosclerosis models.\nReviewer Fig. 2 (S-Fig 4. F-K). F. scRNA UMAP demonstrating SMC cell states reanalyzed from Pan et al. ascending aorta split between 26 weeks of high fat diet on LDLR-KO background and 22 weeks of high fat diet on ApoE-KO background with label transfer from control timecourse scRNA. G. Confusion matrix showing Pan et al. seurat clusters and timecourse cellular profile similarity. H. FeaturePlots showing localization of Fbln1, Fbln2, and Ibsp as representative markers for FMC-1, FMC-2 and CMC. I. scRNA UMAP demonstrating SMC cell states reanalyzed from Alencar et al. Klf4-KO brachiocephalic artery data after 18 weeks of high fat diet after label transfer from control timecourse scRNA. J. Alluvial plot showing change in cell state proportions using transferred labels. CMC decrease from 52.4% to 20.5%. K. Confusion matrix showing Alencar et al. seurat clusters and timecourse cellular profile similarity.\nResults, line 534:\nOur comprehensive multi-omic dataset has enabled the profiling of TF motif accessibility gradients across time and leveraged these data to generate regulatory networks, prioritize key driver genes and evaluate their functional pathways. Analyses using both WOT and pseudotime suggested that both FMC-1 and FMC-2 phenotypes transition probabilities stabilize in mature lesions and are thus more likely to represent sustained phenotypes rather than a transition state. Complementary analyses with pre-existing 20 and 26 week high fat diet (ApoE and Ldlr-KO, respectively) atherosclerosis mouse scRNAseq datasets18 and brachiocephalic derived advanced lesions 17 support the sustained existence of these phenotype cells to late stages of disease. This is an attractive possibility, following the hypothesis that protective cells transition primarily to stable FMC and disease promoting cells undergo further transition to the CMC lineage.\nWe incorporate new methods for analysis of these publicly available data in the Materials and Methods section of the manuscript:\nMaterials and Methods, Line 842:\nProcessing of external datasets Pan et al., Alencar et al., Cheng et al.\nData from Pan et al. (GSE155513)18, Alencar et al (GSE150644)17, Cheng et al. (PRJNA794806)4 were downloaded and reprocessed using standard 10x CellRanger pipeline as noted above to facilitate label transfer analysis with our timecourse dataset. All subsetting parameters and clustering parameters decisions were made based on available published methods, following prior settings as closely as possible. For the Pan et al. dataset, we performed further data pruning, subsetting based on fluorescent marker expression as previously described in Sharma et al.1. We then applied cutoffs of nFeature > 1900, percent.mt < 7.5 to match our control timecourse study. For clustering, we identified top 1500 variable genes and generated UMAP with 20 principal components (PCs) as previously described. For Alencar et al. we filtered on nFeature_RNA > 200 and < 5000, percent.mt < 10. 19 PCs were selected for clustering. The clustering resolution was not specified. This does not affect our downstream analysis, as the goal for this analysis is to understand the distribution of label-transfer clusters from our control timecourse study. For Cheng et al. we requested original processed script from the author to recreate the R-Smad clustering. After clustering of these datasets, we then performed CCA label transfer followed by generation of confusion matrices to compare and contrast the relationship between clusters and identify existing cell states.\nThe authors use WOT algorithm that was initially developed by Eric Lander’s and Aviv Regev’s labs to improve upon already available methods for cell trajectory or lineage inferring across developing tissues. However, all these methods, including the one in this manuscript, lacks the cellular resolution to accurately infer lineage tracing, specially for rarer populations. This seems evident by the apparent lack of origin of FMC2 subcluster at the earlier timepoints (~ 5 Wk). While the Sankey representation of inferred transition of state in Fig 3D clearly shows the transition SMC2 -> FMC1 -> FMC2 or SMC2 -> FMC2 in the later timepoints, it seems to have a problem inferring transition state changes for FMC2 at 5 wk. A similar issue is encountered with CMC population in the Tcf21 KO at 12 and 16 wks. The authors must comment about the temporal and cellular resolution limits of this approach and strategies to possibly improve it (i.e. adding timepoints in between?).\nWe initially attempted to create a unifying transition probability figure for 2C and 3D that that not only showed representative cell proportions at each timepoint but also the transition probability with size of transition bands relative to the starting band width at each time point. Unfortunately, due to problems with rasterization of the figure images, the lines in the Sankey plots connecting FMC-1 to FMC-2 at the 5-week time-point may not be readily visualized. This is due to a small, but present, population of FMC-1 at baseline, reducing the probability of transition of FMC-1 to FMC-2. With this caveat in mind, we have recreated the visualization to just show relative proportion of change in Fig. 2D and emphasize in the text that the actual proportion of each cell state at each time point can be observed in Fig. 2C. With our current visualization, we try to balance relative proportional changes as well as the initial proportion of each cluster. We did not pursue this visualization initially, because we did not want to create the false sense that cell proportions are equal at each timepoint but unfortunately clearly introduced other unintended confusion with this plot. These new plots have now replaced the initial Fig 2C and 3D within the Results section of the manuscript. With this change, the relative proportion changes are more apparent for the FMC1-FMC2 transition in 2D with ~20% probability and FMC-CMC probability in 3C with low probability from both FMC-1 and FMC-2.\nOverall, as the reviewer correctly pointed out, this reflects the important observation that current methods may lack the fidelity to infer transition for more rare populations, which is why we utilized more conservative cutoffs throughout our study to develop a reproducible trajectory model for SMC transitions. Furthermore, this highlights the limitation of available methods where adequate representation of cell states, is necessary due to the inability to perform accurate zero-shot inference while simply adding additional datapoints also makes trajectory inference more difficult due to incremental changes in cellular proportion. Our strategy with available methods straddles the middle ground where we generate multiple timepoints to confirm trends in cellular proportion change and carefully select representative timepoints that show appropriate cell state proportion change to facilitate trajectory analysis.\nResults, line 185:\nTo visualize the WOT predicted relative transition probabilities and emphasize the actual proportion of each cell state at each time point, we have visualized the transition probability results with a Sanke plot (Fig. 2D). The actual proportion of each cell state at each time point can be observed in Fig. 2C and S-Table. 5.\nIn this context, it is also notable that the authors appear to use a more stringent transcript threshold than their prior work (Cheng et al., Nat Card Res. 2022 and Sharma et al., 2024) and this may miss rarer smaller clusters down the lineage lines. Indeed, In Cheng et al., Nat Card Res. 2022, the authors previously identified, albeit at a lower resolution, a smaller R-SMC cluster located in the basal plaque and primarily identified by their high Mmp3 expression, but distinct from SMC, FMC and CMC, and they hypothesized that these derived from FMC by pseudotime analysis with Slingshot. What happened to this SMC subcluster? Was it integrated into the new CMC cluster or was it not detected on this dataset?\nWe thank the reviewer for this helpful comment. We re-examined the dataset from Cheng et al.4, to re-derive the R-SMC cluster. After label transfer from our time-course data, we discovered that the R-SMC population is a subset of the FMC-1 cell state identified in our current study, with the prior Mmp3 RNAScope data showing similar localization consistent with our FMC-1 population (Reviewer Fig. 3). These cells primarily span the SMC media, and at a low level also the fibrous cap of the plaque. When observed on a time-course scale, we see early expression of Mmp3 at 5 weeks, although its expression can also be noted at a low level in the FMC-2. After label transfer, there is still a significantly higher proportion of CMCs in the Smad3 knockout data, demonstrating both the robust translatability as well as unifying ability of our time-course data. We note that there was considerable heterogeneity of cluster definitions in the earlier vascular smooth muscle cell studies and one goal of this time-course study was to provide a more robust and consistent cell state definition for future studies in the field. Therefore, we utilized stringent filtering thresholds to promote reproducibility across future studies.\nReviewer Figure 2. (Left) Label transfer of timecourse data onto SMC subset from Cheng et al shows that R-SMC (Cluster 7) reside in the timecourse FMC-1 cluster. (Right) Confusion matrix showing that cluster 7 expression from Cheng et al. is correlated with the FMC population, and primarily with FMC-1.\nThe authors should show the actual percentage change in FMC2 population between control and Tcf21-KO at 16 weeks? From figure 4A, the overall (FMC-2) change seems negligible when compared to the decrease in FMC1 and CMC. Is the overall increase in atherosclerosis in Tcf21-KO attributed mainly to the expansion of SMC3, given the lack of changes FMC2 population? Another mechanism that might play a role beside the changes in phenotypic switching or percentages of modulated SMC population in Tcf21-KO is their localization. Is there a significant change in the localization of FMC2, FMC1, CMC and particularly SMC3 population in control vs Tcf21 KO?\nThese are great questions, and we are working to address them. The specific roles that these cells play in the development of the lesion anatomy will provide insights into how SMC function to create the disease plaque and fibrous cap. However, we feel that we will need to be able to lineage trace these cells to adequately address the phenotypic changes these cells undergo and contribute to the disease process. We are currently sourcing lineage tracing mice to further investigate the differential contribution of these cells to the disease plaque in mouse atherosclerosis.\nMinor Points\nDoes the osteoblast, senescence, EMT, apoptosis, angiogenesis or efferocytosis cell state of SMC3 changes between control and Tcf21?\nWe compared the module scores in SMC-3 between Tcf21-KO and control using Wilcoxon rank sum test and made the following observations using a Bonferroni corrected p=0.008 as the significance cutoff threshold (for 6 comparisons in SMC-3). These analyses identified a significant decreased senescence module score (p=0.0002) , increased EMT score (p < 2e-16), angiogenesis (p = 3.1e-6), and osteoblast scores (p=1.6e-6), while there was no change in apoptosis (p=0.66) and efferocytosis scores (p=0.015). These data suggest that Tcf21-KO has a direct impact on the expression profile of progenitor cells within the SMC-3 population. However, the direction of module scores may not indicate the directionality of the functional profile of the cell and would require additional study to identify how these populations function as a progenitor population.\nCan the authors provide the actual percentages of the 6 different SMC subpopulations at all timepoints in Fig 2C?\nGiven the number of panels required to show the results for Figure 2, adding percentages within the figure was difficult, so we have added an additional supplemental table (S-Table 5) with cell proportion percentages for each .\nGiven that control TCF21 + IgG showed some unspecific background PLA signal, it would be worth it to show the second control (TEAD1 + IgG) and quantifications of the background signal in the negative controls as well as in TCF21 + TEAD1. It seems that there is also significant interaction between Tcf21p and Tead1p in the cytosol as well. Given this, it would be better to include more information regarding the images. Are they epifluorescence or laser confocal images? If epifluorescence, are they deconvolved and how? Are they max projection or single slices (mentioning objective + NA and pinhole size, or optical slide depth)?\nOur epifluorescence microscopy was performed using the ECHO Revolve with a 20x Olympus Plan Achromat (NA 0.8, WD 10mm), taking wide-field single slice images and noise reduction performed by proprietary Digital Haze Reduction (DHR) software (ECHO).\nPLA assays detects proximity of proteins up to max of 35~ 40 nm, well above the Rmin of the predicted dimer for TcF21-Tead1 interaction, thus it is still not direct proof of binding/interaction between these protein. Can/did the authors try coIP? (acknowledging coIP from nuclear lysates can be difficult). This would provide better evidence of physical interaction between these two protein or coIP their intermediate if the interaction is not direct.\nIn the context of our study, we initially hypothesized that TCF21 and TEAD1 act at similar genomic regions and thus performed the PLA assays. To follow-up on this mechanistic question, we performed Co-IP using a myc-tagged TCF21 transfected in HEK293 cells. We performed nuclear protein extraction followed by IP for the MYC-tag (mouse ab, abcam) and performed western blot for TEAD1 (rabbit ab, abcam) to find evidence for TCF21-TEAD1 physical interaction. These studies also included IgG negative control and 5% input positive controls. These data are now described in the manuscript and presented in the Results section of the manuscript (Fig. 5G).\nResults, line 452:\nTo detect a direct physical interaction between TCF21 and TEAD1, we performed Co-IP using a Myc-tagged TCF21 transfected into HEK293 cells (Fig. 5G). We performed nuclear protein extraction followed by IP for the Myc-tag and western blot for TEAD1. These studies included IgG negative control and 5% input positive controls and provided evidence for TCF21-TEAD1 physical interaction.\nWe incorporated new methods for our Co-IP in the Materials and Methods section of the manuscript:\nMaterials and Methods, Line 972:\nCo-Immunoprecipitation\nCo-IP was performed using a myc-Tagged Tcf21 construct cloned into the PwPI vector. HEK293 cells were transfected with Tcf21-myc followed by protein extraction and nuclear isolation following manufacturer’s instructions (Active Motif Nuclear Complex Co-IP Kit 54001). 5% input control was kept as a positive control and samples were incubated with anti-ms Myc-tag (Cell Signaling 2276S) or anti-ms IgG (Vector I-2000-1) as negative control. Co-IP was performed following instructions from the Active Motif Co-IP Kit, briefly, samples were incubated with 1:50 of Myc-tag antibody according to manufacturer recommendations and 1ug of anti-rb IgG (which is in excess compared to amount of Myc antibody used). After overnight incubation, samples are washed and protein collected for western blotting with anti-rb TEAD1 (Cell Signaling 12292S). An anti-rabbit HRP secondary antibody targeting the light chain was used for primary antibody detection as TEAD1 size is approximately 47kDA which would overlay with heavy chain antibody fragments.\nIt seems that other TF besides the KLF family and AP1 factors, such as Smad2 and Smad3, also show a bimodal open motif distribution. Could they be playing a dual role in early SMC to FMC transition as well?\nThank you for this comment. We utilized for our studies in part the JASPAR2020 motif database that was described in the Pando paper25. After reviewing the motif sets from that version of the database, we find that one instance of the 3 Smad2/Smad3 motifs actually correspond to AP-1 binding sequences (MA1622.1). This error has subsequently been corrected in later versions of JASPAR, and we have corrected our heatmap accordingly, and the heatmap in Fig. 3B has been edited to remove MA1622.1.\n• Figure 4B: label the top Sankey plot as “Control” and bottom as “Tcf21-KO” in the graph, not only mentioning it in the figure legend\nThank you for pointing this out, we have now corrected the labeling in Fig. 4B as noted.\nReviewer #3 (Remarks to the Author)\nThis is a detailed and comprehensive single-cell and epigenomic study on VSMC transitions in atherosclerosis in the context of CAD. It aimed to study trajectories of SMC state transitions, identify novel TFs and gene regulatory networks contributing to these transitions, and how the Tcf21 gene regulates these processes. Key findings included the identification of six SMC states, with the FMC-1 and FMC-2 states contributing to CMC. Tcf21 emerged as the top regulator for the SMC to FMC transition, along with other TFs, ie. Runx, Klf4, Smad3, HIF etc. Tcf21 KO resulted in reduced FMC and CMCs and increased SMCs, and TEAD1 was found to interact with Tcf21 to modulate TGFB and ECM-related profibrotic pathways.\nStrengths:\nThe authors utilized advanced multi-omic tools such as single-cell multi-omics (scRNA+scATAC) and technical computational modeling (ie. WOT, CellOracle, Pando..etc). They conducted lineage traced atherosclerosis mouse model experiments to assess the transition of SMC switching across 4 time points and the role of Tcf21. They assessed translational relevance by integrating with human genetic data (GWAS, scDRS, S-PrediXcan). They discovered a novel mechanism in the Tcf21-TEAD1 interaction and importantly identified FMC2 as a key CAD-risk state.\nWeaknesses:\n1. Three key findings were that 1) Tcf21 plays a strong role as a driver of contractile SMCs transitioning to FMC-1 and that 2) FMC-2 was identified as the key FMC cluster associated with CAD risk. 3) Knocking out Tcf21 led to decreased FMC1 and FMC2 and CMC proportions and increased SMC2 and SMC3 at 12 weeks, and further decrease of FMC1 but not FMC2 at 16 weeks. Since this was the case, did the authors consider another major candidate regulator/TF associated with driving SMCs to FMC2 (ie. TEAD1?) since this seems to be a very important natural point to the study’s narrative?\nWe thank the Reviewer for this important comment. Based on our findings, we interpret the cell state changes from Tcf21 knockout as a direct effect on FMC-1 and an indirect effect on FMC-2 with compensation by week 16. Given our prediction from fate correlation analysis that Tcf21 acts upstream of Zeb2, and experimental evidence confirming that TCF21 binds and regulates ZEB2 expression, we re-analyzed our prior Zeb2 mouse atherosclerosis knockout data22. We performed label transfer and found that there was also a significant decrease in FMC-1 but not FMC-2. We interrogated differentially expressed genes with DESeq2 and identified increases in Tead1 as well as co-activator Yap1 within the FMC-2 cluster, suggesting that these factors are compensating for Tcf21 loss and maintaining FMC-2 cell numbers.\nIn addition, we also further interrogated our control mice gene regulatory networks and were intrigued to identify Tead1 as a downstream target of Zeb2, and Mecom – the two top transcription factors that are predicted to be Tcf21 targets. Both Zeb2 and Mecom were predicted to suppress Tead1 expression in this context which is validated with Tead1 increase in the Zeb2 knockout study.\nMoreover, our repeat analysis of the Klf4 KO data has shown that this factor, which promotes FMC-CMC transition, regulates Tead1 expression. Other transcription factors predicted to be involved with Tead1 include AP-1 factors such as Fos and Atf3, which are core TFs that have a role in regulating smooth muscle identity.\nAs outlined below, we have performed extensive in vitro studies to link Tead1 and the Yap pathway with Tcf21 interactions in the context of SMC transitions, DNA binding and transcription, and collaboration with Tcf21 to regulate disease risk mechanisms. We fully agree with the Reviewer’s hypothesis regarding the SMC transition role for Tead1, however, to extend our understanding of Tead1 in greater detail will require more extensive studies with a direct focus on this factor, employing a cell-specific model and extensive study of TEAD1 contribution to CAD risk, and we feel that these studies are beyond the scope of the current work.\n2. The authors provide important in situ transcript validation using RNAScope. However, additional protein-level validation of FMC markers would strengthen the evidence beyond the RNA landscape towards functional cell phenotypes.\nWe appreciate the reviewer’s suggestion. While we agree that protein-level confirmation would be valuable, we felt that in our hands individual antibody studies are extremely time-consuming and difficult to scale at this point. We are looking at spatial proteomics as a method that will provide the type of data that will be truly transformative for studies such as those noted by the reviewer, and we hope to implement this work in the coming months. Other protein-based imaging (e.g., high throughput immunofluorescence or imaging mass cytometry) will also further strengthen the single cell transcriptomic findings as high‑quality reagents can be identified.\nHowever, to provide robust, orthogonal validation that preserves spatial context, we performed in situ spatial transcriptomics using the Xenium 5k panel, a state‑of‑the‑art platform that quantifies thousands of transcripts with single‑cell/subcellular resolution across multiple segments of the mouse aorta to further validate our nominated cell states (Fig. 1P-T, S-Fig. 2A). To the best of our knowledge, this will be the first such single cell resolution dataset to be generated for mouse aortic atherosclerosis and will significantly aid the research of groups in the field. As another tool that expands the utility of single cell transcriptomics, we leveraged label transfer from our time-course to assign cell‑state identities using the robust cell type decomposition (RCTD) algorithm. These analyses demonstrate strong concordance between the cell states nominated in our discovery data and those observed in situ, further validating the anatomic niches and co-localization patterns described in our study.\nThe evidence of the role of TEAD1 from ChIP-seq and reporter assays support TEAD1 as a TCF21 partner. Functional perturbation of TEAD1 would significantly strengthen the mechanistic conclusions.\nWe performed TEAD1 knockdown with siRNA and evaluated bulk RNAseq in HCASMC. The transcriptomic changes with TEAD1 knockdown showed strong enrichment for processes such as ‘small GTPase mediated signal transduction’, ‘extracellular matrix organization’, ‘wound healing’, ‘cell-substrate adhesion’, and ‘regulation of actin filament-based processes.’ These findings are in line with what we observed from our TEAD1 ChIP-seq functional enrichment, further corroborating the hypothesis that TEAD1 plays a major role in the smooth muscle transition process to affect the development of phenotypically transitioning SMCs. We have now summarized these data in the Results and Discussion sections of the manuscript.\nResults, line 457:\nGiven these findings, we investigated the expressed gene programs in SMC directed by TEAD1 by performing TEAD1 knockdown with siRNA in HCASMC along with bulk RNAseq. The transcriptomic changes with TEAD1 knockdown showed strong enrichment for processes such as ‘cellular response to cAMP, ‘extracellular matrix organization’, ‘connective tissue development”, and ‘chondrocyte development’ (Fig. 5H). These findings are in line with what we observed from our TEAD1 ChIPseq functional enrichment, further corroborating the hypothesis that TEAD1 plays a major role in the smooth muscle transition process to affect the development of phenotypically transitioning SMC.\nIf the conclusion is that Tcf21 promotes the protective FMC1 and competes with TEAD1 which promotes the CAD-risk-enriched state FMC2, it may be important to know whether TEAD1 acts as an independent disease promoting factor or modulates risk through its interaction with TCF21.\nThis is an important question, given that the distribution of TEAD1 expression tracks with the FMC-2 population. First, we investigated summary level Million Veteran Program CAD GWAS data at the TEAD1 locus and noted there is a highly significant peak ~300,000 basepairs from the TEAD1 gene. There are no vascular tissue eQTLs among the CAD associated SNPs, but modest eQTLs identified in testis and blood cells point to the closer ARNTL gene. As we have identified numerous examples where a single disease associated enhancer can regulate expression of multiple genes, we cannot rule out the possibility that variation at this locus regulates TEAD1 and has a direct causal effect on CAD risk. This is obviously an important question but studying a primary role for TEAD1 CAD causality at this point would be a substantial genome editing endeavor. We are working on CRISPR screens of CAD loci genes and will build guides to target the enhancers in this locus for these future studies. Also, we hope to address the question whether this gene promotes or inhibits the FMC-2 state, but this is difficult to address in the absence of SMC eQTLs, and without more detailed study of allelic variation at the TEAD1 locus.\nThus, we can only point to the interaction of TEAD1 with TCF21 as a mechanism of CAD risk, however there is now considerable evidence that this interaction is substantial. We now show new data that these proteins physically interact (Co-IP, Results – line 452, Materials and Methods - line 972, Figure 5G), in addition to colocalization across the genome in discrete enhancer regions (ChIP and PLA), and overlapping CAD GWAS loci enrichment (GWASAnalytics), modulate the binding of each other across pseudotime and promote the transcriptional activity of each other (luciferase assays).\nFurther, in new analyses we provide additional data that TCF21 risk is accentuated through interactions with TEAD1. Employing ChIPseq binding data and GWASAnalytics results for GWAS SNP localization, we found that all TCF21 binding sites, including those with TEAD1 peaks, showed high level enrichment for CAD (-log p-value 10), but in the absence of TEAD1 peaks showed only low-level enrichment for cancer and metabolism (S-Fig. 9F, I). Also, TEAD1 peaks including TCF21 peaks showed strong enrichment for hypertension (-log p-value 10), while TEAD1 peaks only showed low level enrichment of SNPs for diabetes and hypertension (S-Fig. 9G, H). Taken together, these data suggest that the interaction of TCF21 and TEAD1 play a significant specific role toward CAD risk.\nThese analyses are now included in the Results section of the manuscript.\nResults, line 429:\nInterestingly, performing enrichment analysis for GWAS SNPs at TCF21-TEAD1 co-bound loci using ChIPseq binding data and GWAS SNP localization with the GWASAnalytics package, we found that all TCF21 binding sites, including those with TEAD1 peaks, showed high level enrichment for CAD (-log p-value 10) (Fig. 5E, S-Fig. 9F), but in the absence of TEAD1 peaks showed only low level enrichment for cancer and metabolism (S-Fig. 9I). Also, TEAD1 peaks including TCF21 peaks showed strong enrichment for hypertension (-log p-value 10), while TEAD1 peaks only showed low level enrichment of SNPs for diabetes and hypertension (S-Fig. 9G, H). Taken together, these data suggest that the interaction of TCF21 and TEAD1 plays a significant specific role toward CAD risk. This TCF21-TEAD1 relationship corresponded with murine trajectory analysis nominating Tead1 as an FMC-2 driver (Fig. 2L) and making it an intriguing Tcf21-interacting partner for further study. These results for TCF21-TEAD1 were contrasted with similar studies for CEBPB where inclusion of TCF21 loci detracted from the CEBPB metabolism signal (S-Fig. 9J, K). Analyses for HNF1A peaks served as negative control for these studies (S-Fig. 9L, M).\nThe use of CellOracle provided a rich set of candidate TFs for driving SMC phenotypic switching. With the exception of Tcf21, some (even limited) validation of additional TFs could strengthen conclusions made.\nWe thank the reviewer for this recommendation to expand the validation of our analytical algorithms. We performed siRNA knockdown of three highly ranked transcription factors predicted from our CellOracle analysis to affect the transition process, AR, EPAS1 and ZEB1. We performed siRNA knockdown followed by qPCR for contractile markers including ACTA2 and TAGLN as well as FMC markers LUM and PDGFRB. Knockdown of these genes provided robust evidence for increased contractile marker expression across these conditions. Furthermore, we also observed significant downregulation of LUM and PDGFRB for EPAS1 and ZEB1 knockdown while AR knockdown showed a smaller effect size but a trend towards decreased expression levels of these FMC markers. Together, this provides additional functional validation of our predictive modeling and we include these findings in a new supplemental figure (S-Fig. 7A-C).\nResults, line 304:\nTo validate our analytical algorithms, we selected three highly ranked transcription factors predicted from our CellOracle analysis to affect the transition process, AR, EPAS1 and ZEB1 and performed siRNA knockdown followed by qPCR for contractile markers including ACTA2 and TAGLN as well as FMC markers LUM and PDGFRB (S-Fig. 7A-C). Knockdown of these genes provided robust evidence for increased contractile marker expression across these conditions. Furthermore, we also observed significant downregulation of LUM and PDGFRB for EPAS1 and ZEB1 knockdown while AR knockdown showed a smaller effect size but a trend towards decreased expression levels of these FMC markers. Together, this provides additional functional validation of our predictive modeling.\nWe incorporate new methods for in vitro validation of TF SMC function in the Materials and Methods section of the manuscript:\nMaterials and Methods\nsiRNA Experiments, line 936:\nVascular SMCs were transfected with siRNAs targeting TEAD1, AR, EPAS1, and ZEB1(Dharmacon ON-TARGETplus SMARTpool siRNA L-012603-00-0005, L-003400-00-0005, L-004814-00-0005, L-006564-01-0010; Horizon Discovery). Silencer Select negative control siRNA (ThermoFisher 4390844) was used and has been previously tested in our laboratory to ensure no cellular physiology changes. siRNA pool transfection was subsequently performed using Lipofectamine RNAiMax Transfection Reagent (ThermoFisher 13778150) according to manufacturer’s instructions and incubated for 12 hours post-transfection in serum free media followed by an additional 12 hours of recovery in SmBM supplemented media prior to RNA collection.\nRNA Extraction and Quantitative PCR (qPCR) experiments, line 946:\nRNA extraction was performed using RNAEasy Mini kits (Qiagen 74106). 500ng of RNA per sample is then aliquoted for reverse transcription with High-Capacity RNA-to-cDNA Kit (Life Technologies 4388950). Quantitative PCR reactions were conducted with Taqman Universal Master Mix (44440048) and qPCR FAM probes for genes of interest (ThermoFisher) on a QuantStudio 6 Pro Real-Time PCR System. Relative transcript abundance was determined by the comparative Ct (Δ ΔCt) method, using the housekeeping gene UBC.\nSome functional validation in primary SMCs to identify and distinguish FMC1 and FMC2 would strengthen the story inferred from the multi-omics data.\nWe appreciate this suggestion for in vitro functional modeling which can validate the in vivo scRNAseq findings and further facilitate in vitro discovery of molecular mechanisms. To pursue this objective, we used available data from 52 unrelated primary coronary artery smooth muscle cell (HCASMC) “lines” previously published in Liu, et al.26. Using top marker genes from our single cell mouse timecourse, we performed non-negative least squares cellular deconvolution and identified cell lines 2105 and 1508 as exhibiting the top ‘FMC-1-like’ and ‘FMC-2-like’ cell states, respectively (S-Fig. 7D-G). We then treated these primary cell lines with calcium phosphate media24 for seven days to simulate an osteochondrogenic CMC-like transition phenotype. Quantitative PCR identified a greater increase in CMC genes including RUNX2, SOX9, and HAPLN1 in the 1508 cell line, whereas there was minimal change in expression of these genes in the 2105 line (S-Fig. 7I). This is in keeping with our observations that the FMC-2 population has a greater probability of transition to the CMC state. Furthermore, this important finding suggests that there are functional differences in the primary cell lines depending on their origin/disease state and they may harbor differential ability to respond to various simulated vascular stresses. Identification of these primary HCASMC lines with FMC-1 vs FMC-2 phenotype will promote further studies of molecular differences between these cell states. These studies are now described in the revised manuscript as follows.\nResults, line 312:\nTo provide functional validation of the integrity of the FMC-1 versus FMC-2 cellular phenotypes, we characterized cell expression signatures from primary human coronary artery smooth muscle cell (HCASMC) “lines” previously published in Liu, et al.26. Using top marker genes from our single cell mouse timecourse, we performed non-negative least squares cellular deconvolution and identified cell lines 2105 and 1508 as exhibiting the top ‘FMC-1-like’ and ‘FMC-2-like’ cell states, respectively (S-Fig. 7D-G). We then treated these primary cell lines with calcium phosphate media24 for seven days to simulate an osteochondrogenic CMC-like transition phenotype. Quantitative PCR identified a greater increase in CMC genes including RUNX2, SOX9, and HAPLN1 in the 1508 cell line, whereas there was minimal change in expression of these genes in the 2105 line (S-Fig. 7H, I). This is in keeping with our observations that the FMC-2 population has a greater probability of transition to the CMC state. Furthermore, this important finding suggests that there are functional differences in the primary cell lines depending on their origin/disease state and they may harbor differential ability to respond to various simulated vascular stresses. Identification of these primary HCASMC lines with FMC-1 vs FMC-2 phenotype will aid future studies of molecular differences between these cell states.\nWe incorporate new methods for in vitro SMC primary cell line analysis in the Materials and Methods section of the manuscript:\nBulk RNASeq deconvolution, line 880:\nBulk RNASeq fastq files were downloaded from GSE113348 and processed through a uniform pipeline with cutadapt (5.0), STAR (2.7.10b) and featureCounts (2.0.6) for trimming, alignment and generation of count matrices, respectively. Data was then converted to TPM with Kallisto (0.51.1) for comparison between cell lines. Following the prior publication, we removed nine cell lines from the original dataset due to poor data quality for a total of 52 cell lines used. We then applied non-negative least squares deconvolution using the python implementation of nnls from scipy to estimate the contributions of SMC-derived phenotypic clusters based on our tdTomato SMC subset reference. The top 100 marker genes per cluster (ranked by pct.1/pct.2 with at least 25% expression in the reference cluster) were used as features in the cell type signature matrix.\nCalcification assay, line 928:\nPrimary cell lines (1508 and 2105) were split into 6 well plates at 60k cells per well. Cells were allowed to equilibrate over 12 hours prior to media change with SmBM with supplements as noted above. After 24 hours, once cells are noted to be confluent, we applied a calcification media cocktail consisting of SmBM basal media supplemented with 0.4% FBS, (Lonza CC-3182 without additives). Cells were incubated in calcification media for 3 days followed by media change with calcification media and allowed to grow to 7 days prior to RNA collection as previously described24.\nMinor corrections:\nFigure 6 legend: repeat of “and” in title.\nThank you for pointing this out, we have corrected the labeling in the Fig 6 legend as noted.\nReferences for Reviewer Comments\nVascular smooth muscle cell state trajectories mediate molecular mechanisms of coronary disease risk\nResponse to Reviewer Comments:\nReviewer #1 (Remarks to the Author):\nThe authors have responded to all my questions and I have no further concerns. This will be a nice use of single cell technology to identify cellular trajectories.\nWe thank the Reviewer for their insightful suggestions during this process.\nReviewer #2 (Remarks to the Author):\nLi and Kundu and colleagues do an excellent job in responding to all reviewer comments and concerns, have added substantial new genomic and experimental analyses from their own and public data – including elegant spatial data, have performed in vitro SMC functional studies (including clever use of selected human SMC lines of specific molecular phenotype) to validate key new regulators/drivers and SMC-derived cell states, and have revised their model and interpretation of data to reflect more openly the advanced knowledge yet the uncertainties and conflicts in data and possible limitations.\nOverall, this is a very fine paper that adds importantly to the literature. Kudos for this tour de force!\nWe sincerely appreciate these comments by the Reviewer, we express our heartfelt thanks.\nThe remaining concerns are focused on results and discussion – and suggested considerations for enhancements in interpretation, clarity, and what are the future steps in context of the findings, unknowns and limitations.\n1. Although the presentation and framing of complex genomic data findings are very clear mostly the human genetic/GWAS CAD disease inference sections (last 2 paragraphs of results and parts of discussion) remain dense and a little confusing. While appreciating what is meant by genetic analyses linking FMC2 to CAD risk, the framing for the reader initially is to infer FMC2 are associated with CAD risk (many will assume this as higher risk) while the later interpretation makes it clear that it is associated with lower CAD risk. It might make it easier for reader and be clearer if the framing used a little later is used initially and very specifically in terms of protective direction of risk analyses for FMC2 (e.g., lines 573-575). This is particularly important in order to have the reader focus on the uncertainties and surprises of no clear inference from genetic analyses of increased CAD risk for CMCs – and how to resolve the effect of Tcf21 KO producing a decrease in CMCs as well as FMC2s, which is considered in a balanced and open manner in the revised discussion.\nWe appreciate the reviewer highlighting the potential confusion that can arise through the text where we attempt to assimilate multiple concepts related to genetic risk and risk directionality as well as suggestions to improving the framing of our discussion.\nThe assignment of risk in the FMC-2 cellular phenotype is based on scDRS linking genetic risk in SMC loci, as identified with GWAS, with gene expression across the different clusters of the scRNAseq data. This type of analysis does not determine actual direction of effect, and indeed there may be genes that promote and some genes that inhibit disease risk. At an individual locus, if we have evidence which is the causal gene in that locus, i.e., the gene regulating CAD risk, we can use functional genomics data such as expression QTLs to determine whether the disease promoting versus disease protective allele(s) promotes or inhibits expression of this gene. i.e., whether more gene expression leads to greater or lesser risk. An analysis across all genes in the genome that contribute to the FMC-2 phenotype, directionality is complicated, but tools such as PrediXcan, which is the tool used in our studies offers the best approach. It is a foundational Transcriptome-Wide Association Study (TWAS) method that links genetic variants (SNPs) to complex traits/diseases by predicting Genetically Regulated Expression (GReX) (gene expression influenced by genetics) using tissue-specific models, then correlating this imputed expression with disease phenotypes. This analysis thus predicts at a cluster wide level whether genes expressed by these cluster cells are disease inhibitory or promoting. While this is the one of the best current method to ascertain the question of cluster disease directionality there is certainly much room for uncertainty with this model.\nImportantly, this is the case for CAD genes in this cluster such as ZEB2 which has been shown to likely inhibit CAD risk.\nWe have now clarified the explanation of the methods used to determine directionality of the CAD risk, and made clear that risk can be both promoting and inhibitory, early in the results section of these analyses.\nResults, line 483-493:\nBecause scDRS quantifies risk gene enrichment but not directionality, we also computed GWAS risk-weighted average expression by taking average individual gene expression multiplied by its scDRS gene weights and further inferred directionality using updated heritability adjusted S-PrediXcan modeling1 (S-Table 16). Using this framework, we identified FMC-2 having overall statistical enrichment for excess CAD-associated risk gene expression by scDRS (Fig. 6A). Interestingly, neither the transcriptionally similar FMC-1 nor the calcification-associated CMC exhibited significant scDRS enrichment. Next, from the aggregated GWAS risk-weighted expression analysis, of genes with available predicted risk directionality, we calculated greater averaged expression of CAD risk genes in CMC relative to FMC and SMCs. Consistent with this complex biology, FMC-2 cells express genes that both promote or suppress CAD risk with the net predicted risk direction more protective, prompting further investigation of distinguishing expression patterns of these transition phenotypes.\nDiscussion, line 577-586:\nOur high-resolution dataset extends upon previous observations and identifies FMC-2 as the cell state harboring expression of genes that mediate CAD risk. The directionality of this risk is an important consideration, since we have previously shown that TCF21 promotes the FMC phenotype transition and has a protective effect toward risk causality, and it is imperative to know which of the two FMC clusters that we have described mediates this protective effect. The FMC-1 versus FMC-2 cellular phenotypes are quite different and the mechanism for protection could be achieved through a number of different pathways in each cell type. We employed S-PrediXcan for this purpose, which uses composite eQTL data to make this determination (Fig. 6A). This algorithm is able to take GWAS results and predict the effects of each variant on expression levels of genes that are associated with the trait, and it is able to do this for every loci that is associated with the trait throughout the genome.\n2. This work and the excellent work of this group is broadly genetic/epigenetic /genomic/transcriptomic oriented which provides a specific perspective on cell phenotypes, functions, regulators and disease relations. In the discussion, it would be worthwhile to acknowledge existing data and need for future multi-omic approaches that incorporate proteomic data and metabolic/lipidomic data that can add layers of cellular physiology to the genomic analyses, aid in cell type specification, functional inference and ex vivo (human) study as well as incorporate known biology e.g., the important of metabolic state of cells; the emerging knowledge of true SMC-foam cell distribution across multiple SMC-derived cell states but with a lack of understanding of functional and disease consequences for foaminess vs. non -foaminess within specific SMC-derived cell types; such knowledge as well as deep spatial multi-omic cell type knowledge may well help to clarify the uncertainties and unanswered questions and reveal greater functional precision of the cell states inferred in the current work. At the same time, the discussion might tackle how the new knowledge of novel key drivers and regulators of specific cell states and transitions might start to be translated into therapeutics and to propose such non-genetic therapeutic approaches for experimental work.\nWe agree with the Reviewer on these points, and we have incorporated text in the Discussion to address these suggestions. Below, we address components of this comment separately and include manuscript additions at the end.\n… need for future multi-omic approaches that incorporate proteomic data and metabolic/lipidomic data that can add layers of cellular physiology to the genomic analyses, aid in cell type specification, functional inference and ex vivo (human) study as well as incorporate known biology…\nOur focus on transcriptomics and genomics is related to the need for causal information that can only come from human genetics, with our investigation of SMC state transitions described here serving as homework to identify the basic cellular and molecular processes that are the causal mediators of the genetic disease risk. These studies have identified a number of genes and gene programs that regulate SMC cell state changes that are integral to disease risk, and indeed a number of these reside in CAD loci and deserve further validation and mechanistic studies. The overall direction in this laboratory has been the characterization of causal gene regulatory networks that mediate CAD risk, and this goal served as the basis for the studies reported here. The characterization of TCF21 gene regulatory network has led to the identification of additional critical gene and gene programs that can mediate SMC state changes and contribute to SMC phenotypic transitions. These findings will expand the CAD gene regulatory network. However, there are many additional loci that harbor disease risk and will likely require the implementation of expanded and novel approaches with integration of additional data types.\nThe power of orthogonal lipid data was highlighted by the recent paper from Raj Gupta, where his group used an endothelial genetic personalized risk score that was shown to identify individuals who were more likely to benefit from aggressive control of LDL cholesterol level.2 Further, assessing risk at the level of the proteome, which stands between genetic and environmental risk factors, will add significant support to quantifying overall disease risk and point to attractive avenues for manipulating disease mechanisms. Continued integration of CyTOF methods and CiteSeq will improve fidelity and robustness of genetic findings. Indeed we used CiteSeq in our first single cell study of atherosclerosis.3 Integration of single cell methylation data into the genetic workflow will also add significantly to assessment of risk and support identification of disease pathways that introduce novel therapeutic considerations. We have now discussed these considerations in a new section of the Discussion where we make these point suggested by the Reviewer, as detailed below.\n...emerging knowledge of true SMC-foam cell distribution across multiple SMC-derived cell states but with a lack of understanding of functional and disease consequences for foaminess vs. non -foaminess within specific SMC-derived cell types…\nThe field of macrophage function and contribution to disease risk is a complex area that is currently the target of much interest and active research. While considerable single cell work has established the paradigm that foam cells are not associated with disease while inflammatory macrophages are highly correlated with active disease, there is some evidence for a pathogenic macrophage transition from a foamy cell phenotype to an inflammatory lipid-associated phenotype.4 These considerations highlight the profound need for detailed definitive studies that finally resolve the controversy regarding whether vascular SMC can transition to the foam cell phenotype in the atherosclerotic plaque. Initial speculation that SMC could transition to a macrophage phenotype was understandable given extensive data that modulated SMC express a number of macrophage genes such as Lgals3, but advent of single cell studies clearly indicate that this is a very rare event.3,5 Whether SMC transition to foam cells has been more complicated. Lipid loading studies with SMC has shown that they can take up lipid in vitro, and we and others have shown that they can take up lipid in vivo3 but the probability and to what extent vascular SMC can transition to the foam cell phenotype remains unknown. The literature is replete with primary research and review papers discussing and attempting to characterize the contribution of such a cell type to plaque, and indeed its role in disease risk could be substantial and the avenues for manipulation important for therapeutic considerations substantial, but this cell must be characterized with modern genetic and genomic features. A recent comprehensive single cell study of carotid plaque has identified cells with gene expression features of both SMC and macrophage cellular phenotype, and further study of this cluster could address the current need in the field.6 Without question, identification of the pathway by which SMC lineage foam cells might take up lipid and contribute to destabilizing the plaque could have significant importance for targeted therapy development. We now make these points in the Discussion section of the manuscript, as noted below.\n“the discussion might tackle how the new knowledge of novel key drivers and regulators of specific cell states and transitions might start to be translated into therapeutics”\nThere has been considerable discussion regarding approaches that might allow the implementation of therapeutic strategies targeted to the SMC lineage.7,8 However, we have to argue that at this point the field is just not ready to explore this direction, although we are getting closer to making this a reality. It is increasingly clear that SMC can either promote disease risk or inhibit risk, depending on its transitional phenotype. We have learned much from integrating human genetic information and complex mouse genetic models which suggest that increasing transition to the fibromyocyte number will stabilizing this phenotype, while attenuating transition to the CMC phenotype and related intimal calcification will stabilize the intimal plaque and promote fibrous cap formation. The position of this lab at this time, is that we are still in the discovery phase of SMC pathobiology, and that sound approaches to SMC therapeutics will depend on sustained thoughtful experimentation with open data sharing and multidisciplinary collaboration, the combination of lipidomic, metabolomic, and other disease related single cell omics platforms.\nDiscussion, lines 675-696:\nWe are compelled to note that our extensive characterization of SMC phenotypic transitions does not provide evidence that this cell lineage can adopt the macrophage phenotype. Initial speculation that SMC could transition to a macrophage phenotype was understandable given extensive data that lineage traced SMC express a number of macrophage genes such as Lgals3, but single cell studies have indicated that these are a rare events.3 Whether SMC transition to foam cells is currently understudied. Lipid loading studies with SMC have shown that they can take up lipid in vitro,9 and we and others have shown that they can take up lipid in vivo,10.3 but the probability and to what extent vascular SMC can transition to the foam cell phenotype remains unknown. SMC foam cells are reported to have a different phenotype than macrophage foam cells, due to relative deficiency of lysosomal acid lipase in SMC, and retain lipid droplets in their lysosomes rather than in the cytoplasm as seen in macrophages11. This would predict SMC foam cells have an altered gene expression pattern and phenotype different from macrophage foam cells. The literature is replete with discussions regarding the impact of SMC-derived foam cells to plaque, and indeed the role of such cells toward disease risk could be significant and the avenues for manipulation important for therapeutic considerations substantial, but this cell must be characterized with modern genetic and genomic methods. A recent comprehensive single cell study of carotid plaque has identified cells with gene expression features of both SMC and macrophage cellular phenotype, and further study of this cluster could address the current need in the field. 6 Without question, identification of the molecular pathways by which SMC lineage foam cells take up lipid and contribute to destabilizing the plaque could have significant importance for targeted therapy development. Moreover, additional multi-omic approaches incorporating DNA methylation, proteomic, metabolic and lipidomic data will allow us to better interrogate cellular physiology in the context of genomic analyses.\n3. The rigor of lineage tracing settings and cell annotations and clustering is appreciated and important for the field yet likely limits identification of cell states/transitions/cell-cell interactions that exist within the currently described populations. This has consequences – there may be critically important rare/uncommon cell states that are missed (e.g., SMC/FMC as pioneer or primed cells existing in normal media before and early into effect of diet on genetic background) and the more nuanced gene signatures in local cell populations that reflect metabolic or lipid stress may be missed. These issues should be acknowledged while also emphasizing the focus on rigor and reproducibility with stringent parameter settings.\nAgreed, we use stringent cluster methods based on a rigorous workflow that makes certain each cell cluster phenotype is homogeneous and different from other clusters in our analyses. We carefully identify the group of genes that are most representative of the variance across the dataset, and number of principle components that best describe the identified variation. Further, we check the gene ontology terms that are identified by cluster genes, to make sure there is minimal overlap between clusters. The field of single cell transcriptomics is rife with examples of over-clustering and claims of “new progenitor cell or intermediate phenotypes” associated with complex disease processes. But, as noted by the Reviewer, this rigor is sustained at the risk of missing interesting rare/uncommon cell states. In that regard, we have identified early markers of the FMC-1 and FMC-2 that identify these cells at early timepoints and transgenic Cre lines using these markers show differential expansion and migration of these true progenitor cell populations. We are currently mapping the cell state changes that are represented by these genes and their effect on SMC contribution to the plaque and fibrous cap. These studies are beyond the scope of the current work and will take considerable experiments to complete, but it is clear that this type of experiment will be highly instructive regarding the early developmental events for the FMC clusters. We adapt the text to clearly acknowledge the likelihood for the presence of such progenitor genes and their utility in understanding the molecular and cellular mechanisms of FMC-1 vs FMC-2 development and contribution to disease processes.\nWe further agree with the reviewer that -- more nuanced gene signatures in local cell populations that reflect metabolic or lipid stress may be missed. These issues should be acknowledged while also emphasizing the focus on rigor and reproducibility with stringent parameter settings.\nDiscussion, lines 665-674:\nOur study uses a rigorous workflow that identifies reproducible cell cluster phenotypes. However, this rigor is implemented at the risk of losing rare cell states including those that respond to metabolic or lipid stress and most notably, phenotypes that represent progenitor cells such as those that undergo clonal expansion in the disease setting. Future studies will be directed at identification of these putative precursor cells of interest utilizing novel lineage tracing methods (e.g. in vivo barcoding, somatic mutation detection, or novel lineage tracers using putative driver genes) that can better capture these rare populations while using our dataset as a source for validation. This way, descriptive characterization of novel cell phenotypes coupled with experiments using such lineage tracing will allow robust characterization of how these precursor cells contribute to the plaque and fibrous cap, and how the deletion of specific markers alter the course of disease.\n4. Mouse models cannot be relied upon to prove causality of genes in human CAD and CVD – they can be supportive. The statement in discussion on this should be modified. In this context, in the discussion, it would be very helpful for the reader if there was a brief discussion of the importance of complementary (to GWAS / regulatory etc) human genetic studies of rare coding variation in genes to infer causality and directionality in CAD/CVD of inferred and proved loss of function and gain of function mutations – and to note the state of this work related to the current paper e.g., for TCF21 and other proposed drivers and regulatory genes. Again, such data will be a critical complement to understand clear causality in humans. There is value also to acknowledging the need for studies and emerging data on somatic mutations in SMC-derived cells in lesions to understand specific causal SMC genes under the hypothesis of somatic variation.\nWe certainly agree that the mouse models cannot prove causality of specific genes and gene programs, but they can provide direct evidence for disease mechanisms by which human causal CAD genes mediate risk-related genetic and molecular processes in the vessel wall. We further agree that rare coding variation can provide highly useful information, occasionally providing a human disease phenotype and direction of effect for CAD genes that are identified with GWAS studies. A good recent example is the MFGE8 gene which has inframe mutations identified in the FinnGen biobank where these variants are present at high frequency and the inactivation of MFGE8 function has been correlated with greater lifespan.12 It is important to note of course that rare coding variation in the protective gene TCF21 would be expected to increase the risk for CAD, and indeed likely cardiovascular developmental defects that may be inconsistent with fetal survival. Also, TCF21 has important developmental roles in the kidney and lung, and mice lacking Tcf21 die at birth due to failure of lung function. This is the likely reason that we have not found coding region mutations in the TCF21 gene that are associated with protein function. Such mutations are likely removed from the genome through purifying selection.\nFurther, along the lines of using rare coding variation to inform on causality and direction of effect for candidate genes that reside in associated loci for complex human disease, a recently described approach from the Pritchard lab employs loss of function burden tests along with relevant Perturb-seq data to bridge the gap between genetic association and biological mechanism.13 By combining these two forms of data, their approach builds causal graphs in which the directional associations of genes with a trait can be explained by their regulatory effects on gene programs.\nThe possibility of using human somatic mutation in SMC to implicate genes of relevance to vascular disease is very intriguing. The NIH program Somatic Mosaicism Across Human Tissues (SMaHT) Network is focused on cataloging naturally occurring DNA mutations (somatic mosaicism) in healthy subjects to understand aging and disease. We have been in touch with these investigators to determine possible collaborations. While they are certainly interested, they are not actively studying any vascular tissue. Given the logistics of collecting and evaluating tissue from multiple humans and then pursuing a very deep sequencing effort, we have not embarked on these studies. Nevertheless, we are anxious to see such an effort move forward and will mention this need in the Discussion section of the manuscript.\nTo address the Reviewer’s comments we have made the following edits to the manuscript:\nDiscussion, lines 611-638:\nWe have undertaken these studies in order to better understand cellular and molecular mechanisms of CAD risk, by identifying and characterizing individual genes and gene programs that modulate SMC phenotype transitions. While such studies do not prove causality, they inform on possible mechanisms of disease risk that reside in the vessel wall. Although not explored in the work discussed here, there are a number of approaches related to mapping causality in human data that warrant further investigation, and can serve to strengthen the genetic data derived in our single cell studies of mouse models. For instance, more detailed study of rare coding variation, when present, can provide human disease phenotype and direction of effect for CAD genes that are identified with GWAS studies. A good example is the CAD GWAS MFGE8 (lactadherin) gene which was identified in the FinnGen biobank to have an association between an inframe insertion rs534125149 and protection against coronary atherosclerosis. 12 Along this line of reasoning, a recently described approach from the Pritchard lab employs loss of function burden tests along with relevant Perturb-seq data to bridge the gap between genetic association and biological mechanism.13 By combining these two forms of data, their approach builds causal graphs in which the directional associations of genes with a trait can be explained by their regulatory effects on gene programs. It is important to note of course that rare coding variation in the protective gene TCF21 would be expected to increase the risk for CAD, and indeed likely cardiovascular developmental defects that may be inconsistent with fetal survival. Also, TCF21 has important developmental roles in the kidney and lung, and mice lacking Tcf21 die at birth due to failure of lung function. This is one possible reason that we have not found coding region mutations in the TCF21 gene that are associated with protein function. Such mutations are likely removed from the genome through purifying selection.\nAnother approach to enriching the pool of CAD causal genes in SMC transitions is the possibility of using human somatic mutation to implicate genes of relevance to vascular disease. The NIH program Somatic Mosaicism Across Human Tissues (SMaHT) Network is focused on cataloging naturally occurring DNA mutations (somatic mosaicism) in healthy subjects to understand aging and disease. This network is not currently studying arterial tissues, and the tasks of collecting and subjecting tissue from multiple humans and then pursuing a very deep sequencing effort is daunting, but possibly well worth pursuing.\nReviewer #3 (Remarks to the Author):\nThis revised version is marked improvement upon the first submission of the manuscript as the authors have provided more evidence to support the identification of the SMC cell states, the role of Tcf21 as a top regulator of SMC transition to FMC/CMC, and TEAD1 as novel interacting partner of TCF21 at enhancer regions. The addition of spatial Xenium data from ApoE-/- atherosclerotic mice with SMC lineage tracing integrated with single-cell clusters and RNAScope analysis allows for orthogonal validation at high resolution of the previously identified SMC clusters. The authors have also included in vitro perturbation assays to validate their algorithmic findings on highly ranked transcription factors that affect the transition process. Knockdown of the selected TFs resulted in a more contractile SMC state. They also performed additional functional experiments on primary HCASMCs to validate that the FMC-2 population had greater probability to transition to the CMC state.\nWe thank the Reviewer for these positive comments.\nMajor comments:\n1. The study utilizes RNAScope, single-cell RNA-seq, spatial transcriptomics, but protein-level confirmation is still limited. The authors mention that individual antibody-based studies are time-consuming and difficult to scale. It would be helpful to at least stain for critical proteins of interest (ie. TCF21, TEAD1) in the lineage-traced atherosclerotic mouse sections (showing Myh11).\nThank you for this suggestion. Here, we stained for TEAD1 and TCF21 using aortic root serial sections from our Myh11 lineage-tracing mice at 16weeks high fat diet to show protein localization.\nWe trialed multiple conditions as below:\nIF Dilutions Antibody Primary Incubation Antigen\n(Overnight) retrieval\nTead1 1:50; Cell signaling D9X2L; 4C; room No\n1:100 Abcam 133533 temperature\nTcf21 1:50; Abcam 32981; 4C; room No\n1:100 Abcam182134; LS temperature\nBio B3504; Sigma\nHPA13189\nIHC\nTead1 1:100; Cell signaling -D9X2L; 4C; room Yes\n1:300 Abcam 133533 temperature\nTcf21 1:100; Abcam 32981; 4C; room Yes\n1:300 Abcam182134; LS temperature\nBio B3504; Sigma\nHPA013189\nWe were unable to obtain adequate images by immunofluorescence. With antigen retrieval, we were successful in obtaining signals using immunohistochemistry. Unfortunately, this results in loss of the lineage tracing marker. Images below are shown for Tead1 and Tcf21 in serial sections with minimal background as indicated by Rb-IgG control. We have incorporated this as a new supplemental figure (See below and S-Fig. 10).\nS-Fig 10. Immunohistochemistry demonstrating TEAD1, TCF21 protein localization with corresponding IgG control.\nMinor comments:\n1. A figure defining the key marker gene panel for identifying each cluster or top expressed genes would be important to show early. In figure 2, there are key drivers for FMC1 and FMC2 switching but identification markers would be helpful.\nThank you for pointing this out. We agree that early presentation of top marker genes are important. As such, we had included marker genes as an early supplemental figure (S-Fig 1L) in the original edit. We note that because of the high number of markers and space constraints, we felt it was a better fit in the supplement where readers can review both the dot plot (S-Fig 1L) as well as table (S-Table 2) in detail as multiple markers should be utilized to define these complex cell states with high fidelity. This was described early in lines 113-114.\n2. Please amend Figure 1A so it includes spatial transcriptomics using Xenium.\nWe have amended our Figure 1A as suggested.\n3. Having an understanding of the UMAP generated by the Xenium dataset would be important to show the clusters alongside their spatial distribution which is depicted.\nAs recommended, we generated corresponding spatial UMAPs with our transferred labels based on the transcriptomic data, demonstrating distinct clustering. We have incorporated this data as S-Fig 2B.\nReviewer References:"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have responded to all my questions and I have no further concerns. This will be a nice use of single cell technology to identify cellular trajectories.", "reviewer_2": "Li and Kundu and colleagues do an excellent job in responding to all reviewer comments and concerns, have added substantial new genomic and experimental analyses from their own and public data – including elegant spatial data, have performed in vitro SMC functional studies (including clever use of selected human SMC lines of specific molecular phenotype) to validate key new regulators/drivers and SMC-derived cell states, and have revised their model and interpretation of data to reflect more openly the advanced knowledge yet the uncertainties and conflicts in data and possible limitations.\nOverall, this is a very fine paper that adds importantly to the literature. Kudos for this tour de force!\nThe remaining concerns are focused on results and discussion – and suggested considerations for enhancements in interpretation, clarity, and what are the future steps in context of the findings, unknowns and limitations.\n1. Although the presentation and framing of complex genomic data findings are very clear mostly the human genetic/GWAS CAD disease inference sections (last 2 paragraphs of results and parts of discussion) remain dense and a little confusing. While appreciating what is meant by genetic analyses linking FMC2 to CAD risk, the framing for the reader initially is to infer FMC2 are associated with CAD risk (many will assume this as higher risk) while the later interpretation makes it clear that it is associated with lower CAD risk. It might make it easier for reader and be clearer if the framing used a little later is used initially and very specifically in terms of protective direction of risk analyses for FMC2 (e.g., lines 573-575). This is particularly important in order to have the reader focus on the uncertainties and surprises of no clear inference from genetic analyses of increased CAD risk for CMCs – and how to resolve the effect of Tcf21 KO producing a decrease in CMCs as well as FMC2s, which is considered in a balanced and open manner in the revised discussion.\n2. This work and the excellent work of this group is broadly genetic/epigenetic /genomic/transcriptomic oriented which provides a specific perspective on cell phenotypes, functions, regulators and disease relations. In the discussion, it would be worthwhile to acknowledge existing data and need for future multi-omic approaches that incorporate proteomic data and metabolic/lipidomic data that can add layers of cellular physiology to the genomic analyses, aid in cell type specification, functional inference and ex vivo study as well as incorporate known biology e.g., the important of metabolic state of cells; the emerging knowledge of true SMC-foam cell distribution across multiple SMC-derived cell states but with a lack of understanding of functional and disease consequences for foaminess vs. non -foaminess within specific SMC-derived cell types; such knowledge as well as deep spatial multi-omic cell type knowledge may well help to clarify the uncertainties and unanswered questions and reveal greater functional precision of the cell states inferred in the current work. At the same time, the discussion might tackle how the new knowledge of novel key drivers and regulators of specific cell states and transitions might start to be translated into therapeutics and to propose such non-genetic therapeutic approaches for experimental work.\n3. The rigor of lineage tracing settings and cell annotations and clustering is appreciated and important for the field yet likely limits identification of cell states/transitions/cell-cell interactions that exist within the currently described populations. This has consequences – there may be critically important rare/uncommon cell states that are missed (e.g., SMC/FMC as pioneer or primed cells existing in normal media before and early into effect of diet on genetic background) and the more nuanced gene signatures in local cell populations that reflect metabolic or lipid stress may be missed. These issues should be acknowledged while also emphasizing the focus on rigor and reproducibility with stringent parameter settings.\n4. Mouse models cannot be relied upon to prove causality of genes in human CAD and CVD – they can be supportive. The statement in discussion on this should be modified. In this context, in the discussion, it would be very helpful for the reader if there was a brief discussion of the importance of complementary (to GWAS / regulatory etc) human genetic studies of rare coding variation in genes to infer causality and directionality in CAD/CVD of inferred and proved loss of function and gain of function mutations – and to note the state of this work related to the current paper e.g., for TCF21 and other proposed drivers and regulatory genes. Again, such data will be a critical complement to understand clear causality in humans. There is value also to acknowledging the need for studies and emerging data on somatic mutations in SMC-derived cells in lesions to understand specific causal SMC genes under the hypothesis of somatic variation.", "reviewer_3": "This revised version is marked improvement upon the first submission of the manuscript as the authors have provided more evidence to support the identification of the SMC cell states, the role of Tcf21 as a top regulator of SMC transition to FMC/CMC, and TEAD1 as novel interacting partner of TCF21 at enhancer regions. The addition of spatial Xenium data from ApoE-/- atherosclerotic mice with SMC lineage tracing integrated with single-cell clusters and RNAScope analysis allows for orthogonal validation at high resolution of the previously identified SMC clusters. The authors have also included in vitro perturbation assays to validate their algorithmic findings on highly ranked transcription factors that affect the transition process. Knockdown of the selected TFs resulted in a more contractile SMC state. They also performed additional functional experiments on primary HCASMCs to validate that the FMC-2 population had greater probability to transition to the CMC state.\nMajor comments:\n1. The study utilizes RNAScope, single-cell RNA-seq, spatial transcriptomics, but protein-level confirmation is still limited. The authors mention that individual antibody-based studies are time-consuming and difficult to scale. It would be helpful to at least stain for critical proteins of interest (ie. TCF21, TEAD1) in the lineage-traced atherosclerotic mouse sections (showing Myh11).\nMinor comments:\n1. A figure defining the key marker gene panel for identifying each cluster or top expressed genes would be important to show early. In figure 2, there are key drivers for FMC1 and FMC2 switching but identification markers would be helpful.\n2. Please amend Figure 1A so it includes spatial transcriptomics using Xenium.\n3. Having an understanding of the UMAP generated by the Xenium dataset would be important to show the clusters alongside their spatial distribution which is depicted."}, "author_response": "Response to Reviewer Comments:\nReviewer #1 (Remarks to the Author)\nThe authors have responded to all my questions and I have no further concerns. This will be a nice use of single cell technology to identify cellular trajectories.\nWe thank the Reviewer for their insightful suggestions during this process.\nReviewer #2 (Remarks to the Author):\nLi and Kundu and colleagues do an excellent job in responding to all reviewer comments and concerns, have added substantial new genomic and experimental analyses from their own and public data – including elegant spatial data, have performed in vitro SMC functional studies (including clever use of selected human SMC lines of specific molecular phenotype) to validate key new regulators/drivers and SMC-derived cell states, and have revised their model and interpretation of data to reflect more openly the advanced knowledge yet the uncertainties and conflicts in data and possible limitations.\nOverall, this is a very fine paper that adds importantly to the literature. Kudos for this tour de force!\nWe sincerely appreciate these comments by the Reviewer, we express our heartfelt thanks.\nThe remaining concerns are focused on results and discussion – and suggested considerations for enhancements in interpretation, clarity, and what are the future steps in context of the findings, unknowns and limitations.\n1. Although the presentation and framing of complex genomic data findings are very clear mostly the human genetic/GWAS CAD disease inference sections (last 2 paragraphs of results and parts of discussion) remain dense and a little confusing. While appreciating what is meant by genetic analyses linking FMC2 to CAD risk, the framing for the reader initially is to infer FMC2 are associated with CAD risk (many will assume this as higher risk) while the later interpretation makes it clear that it is associated with lower CAD risk. It might make it easier for reader and be clearer if the framing used a little later is used initially and very specifically in terms of protective direction of risk analyses for FMC2 (e.g., lines 573-575). This is particularly important in order to have the reader focus on the uncertainties and surprises of no clear inference from genetic analyses of increased CAD risk for CMCs – and how to resolve the effect of Tcf21 KO producing a decrease in CMCs as well as FMC2s, which is considered in a balanced and open manner in the revised discussion.\nWe appreciate the reviewer highlighting the potential confusion that can arise through the text where we attempt to assimilate multiple concepts related to genetic risk and risk directionality as well as suggestions to improving the framing of our discussion.\nThe assignment of risk in the FMC-2 cellular phenotype is based on scDRS linking genetic risk in SMC loci, as identified with GWAS, with gene expression across the different clusters of the scRNAseq data. This type of analysis does not determine actual direction of effect, and indeed there may be genes that promote and some genes that inhibit disease risk. At an individual locus, if we have evidence which is the causal gene in that locus, i.e., the gene regulating CAD risk, we can use functional genomics data such as expression QTLs to determine whether the disease promoting versus disease protective allele(s) promotes or inhibits expression of this gene. i.e., whether more gene expression leads to greater or lesser risk. An analysis across all genes in the genome that contribute to the FMC-2 phenotype, directionality is complicated, but tools such as PrediXcan, which is the tool used in our studies offers the best approach. It is a foundational Transcriptome-Wide Association Study (TWAS) method that links genetic variants (SNPs) to complex traits/diseases by predicting Genetically Regulated Expression (GReX) (gene expression influenced by genetics) using tissue-specific models, then correlating this imputed expression with disease phenotypes. This analysis thus predicts at a cluster wide level whether genes expressed by these cluster cells are disease inhibitory or promoting. While this is the one of the best current method to ascertain the question of cluster disease directionality there is certainly much room for uncertainty with this model.\nImportantly, this is the case for CAD genes in this cluster such as ZEB2 which has been shown to likely inhibit CAD risk.\nWe have now clarified the explanation of the methods used to determine directionality of the CAD risk, and made clear that risk can be both promoting and inhibitory, early in the results section of these analyses.\nResults, line 483-493:\nBecause scDRS quantifies risk gene enrichment but not directionality, we also computed GWAS risk-weighted average expression by taking average individual gene expression multiplied by its scDRS gene weights and further inferred directionality using updated heritability adjusted S-PrediXcan modeling1 (S-Table 16). Using this framework, we identified FMC-2 having overall statistical enrichment for excess CAD-associated risk gene expression by scDRS (Fig. 6A). Interestingly, neither the transcriptionally similar FMC-1 nor the calcification-associated CMC exhibited significant scDRS enrichment. Next, from the aggregated GWAS risk-weighted expression analysis, of genes with available predicted risk directionality, we calculated greater averaged expression of CAD risk genes in CMC relative to FMC and SMCs. Consistent with this complex biology, FMC-2 cells express genes that both promote or suppress CAD risk with the net predicted risk direction more protective, prompting further investigation of distinguishing expression patterns of these transition phenotypes.\nDiscussion, line 577-586:\nOur high-resolution dataset extends upon previous observations and identifies FMC-2 as the cell state harboring expression of genes that mediate CAD risk. The directionality of this risk is an important consideration, since we have previously shown that TCF21 promotes the FMC phenotype transition and has a protective effect toward risk causality, and it is imperative to know which of the two FMC clusters that we have described mediates this protective effect. The FMC-1 versus FMC-2 cellular phenotypes are quite different and the mechanism for protection could be achieved through a number of different pathways in each cell type. We employed S-PrediXcan for this purpose, which uses composite eQTL data to make this determination (Fig. 6A). This algorithm is able to take GWAS results and predict the effects of each variant on expression levels of genes that are associated with the trait, and it is able to do this for every loci that is associated with the trait throughout the genome.\n2. This work and the excellent work of this group is broadly genetic/epigenetic /genomic/transcriptomic oriented which provides a specific perspective on cell phenotypes, functions, regulators and disease relations. In the discussion, it would be worthwhile to acknowledge existing data and need for future multi-omic approaches that incorporate proteomic data and metabolic/lipidomic data that can add layers of cellular physiology to the genomic analyses, aid in cell type specification, functional inference and ex vivo (human) study as well as incorporate known biology e.g., the important of metabolic state of cells; the emerging knowledge of true SMC-foam cell distribution across multiple SMC-derived cell states but with a lack of understanding of functional and disease consequences for foaminess vs. non -foaminess within specific SMC-derived cell types; such knowledge as well as deep spatial multi-omic cell type knowledge may well help to clarify the uncertainties and unanswered questions and reveal greater functional precision of the cell states inferred in the current work. At the same time, the discussion might tackle how the new knowledge of novel key drivers and regulators of specific cell states and transitions might start to be translated into therapeutics and to propose such non-genetic therapeutic approaches for experimental work.\nWe agree with the Reviewer on these points, and we have incorporated text in the Discussion to address these suggestions. Below, we address components of this comment separately and include manuscript additions at the end.\n… need for future multi-omic approaches that incorporate proteomic data and metabolic/lipidomic data that can add layers of cellular physiology to the genomic analyses, aid in cell type specification, functional inference and ex vivo (human) study as well as incorporate known biology…\nOur focus on transcriptomics and genomics is related to the need for causal information that can only come from human genetics, with our investigation of SMC state transitions described here serving as homework to identify the basic cellular and molecular processes that are the causal mediators of the genetic disease risk. These studies have identified a number of genes and gene programs that regulate SMC cell state changes that are integral to disease risk, and indeed a number of these reside in CAD loci and deserve further validation and mechanistic studies. The overall direction in this laboratory has been the characterization of causal gene regulatory networks that mediate CAD risk, and this goal served as the basis for the studies reported here. The characterization of TCF21 gene regulatory network has led to the identification of additional critical gene and gene programs that can mediate SMC state changes and contribute to SMC phenotypic transitions. These findings will expand the CAD gene regulatory network. However, there are many additional loci that harbor disease risk and will likely require the implementation of expanded and novel approaches with integration of additional data types.\nThe power of orthogonal lipid data was highlighted by the recent paper from Raj Gupta, where his group used an endothelial genetic personalized risk score that was shown to identify individuals who were more likely to benefit from aggressive control of LDL cholesterol level.2 Further, assessing risk at the level of the proteome, which stands between genetic and environmental risk factors, will add significant support to quantifying overall disease risk and point to attractive avenues for manipulating disease mechanisms. Continued integration of CyTOF methods and CiteSeq will improve fidelity and robustness of genetic findings. Indeed we used CiteSeq in our first single cell study of atherosclerosis.3 Integration of single cell methylation data into the genetic workflow will also add significantly to assessment of risk and support identification of disease pathways that introduce novel therapeutic considerations. We have now discussed these considerations in a new section of the Discussion where we make these point suggested by the Reviewer, as detailed below.\n...emerging knowledge of true SMC-foam cell distribution across multiple SMC-derived cell states but with a lack of understanding of functional and disease consequences for foaminess vs. non -foaminess within specific SMC-derived cell types…\nThe field of macrophage function and contribution to disease risk is a complex area that is currently the target of much interest and active research. While considerable single cell work has established the paradigm that foam cells are not associated with disease while inflammatory macrophages are highly correlated with active disease, there is some evidence for a pathogenic macrophage transition from a foamy cell phenotype to an inflammatory lipid-associated phenotype.4 These considerations highlight the profound need for detailed definitive studies that finally resolve the controversy regarding whether vascular SMC can transition to the foam cell phenotype in the atherosclerotic plaque. Initial speculation that SMC could transition to a macrophage phenotype was understandable given extensive data that modulated SMC express a number of macrophage genes such as Lgals3, but advent of single cell studies clearly indicate that this is a very rare event.3,5 Whether SMC transition to foam cells has been more complicated. Lipid loading studies with SMC has shown that they can take up lipid in vitro, and we and others have shown that they can take up lipid in vivo3 but the probability and to what extent vascular SMC can transition to the foam cell phenotype remains unknown. The literature is replete with primary research and review papers discussing and attempting to characterize the contribution of such a cell type to plaque, and indeed its role in disease risk could be substantial and the avenues for manipulation important for therapeutic considerations substantial, but this cell must be characterized with modern genetic and genomic features. A recent comprehensive single cell study of carotid plaque has identified cells with gene expression features of both SMC and macrophage cellular phenotype, and further study of this cluster could address the current need in the field.6 Without question, identification of the pathway by which SMC lineage foam cells might take up lipid and contribute to destabilizing the plaque could have significant importance for targeted therapy development. We now make these points in the Discussion section of the manuscript, as noted below.\n“the discussion might tackle how the new knowledge of novel key drivers and regulators of specific cell states and transitions might start to be translated into therapeutics”\nThere has been considerable discussion regarding approaches that might allow the implementation of therapeutic strategies targeted to the SMC lineage.7,8 However, we have to argue that at this point the field is just not ready to explore this direction, although we are getting closer to making this a reality. It is increasingly clear that SMC can either promote disease risk or inhibit risk, depending on its transitional phenotype. We have learned much from integrating human genetic information and complex mouse genetic models which suggest that increasing transition to the fibromyocyte number will stabilizing this phenotype, while attenuating transition to the CMC phenotype and related intimal calcification will stabilize the intimal plaque and promote fibrous cap formation. The position of this lab at this time, is that we are still in the discovery phase of SMC pathobiology, and that sound approaches to SMC therapeutics will depend on sustained thoughtful experimentation with open data sharing and multidisciplinary collaboration, the combination of lipidomic, metabolomic, and other disease related single cell omics platforms.\nDiscussion, lines 675-696:\nWe are compelled to note that our extensive characterization of SMC phenotypic transitions does not provide evidence that this cell lineage can adopt the macrophage phenotype. Initial speculation that SMC could transition to a macrophage phenotype was understandable given extensive data that lineage traced SMC express a number of macrophage genes such as Lgals3, but single cell studies have indicated that these are a rare events.3 Whether SMC transition to foam cells is currently understudied. Lipid loading studies with SMC have shown that they can take up lipid in vitro,9 and we and others have shown that they can take up lipid in vivo,10.3 but the probability and to what extent vascular SMC can transition to the foam cell phenotype remains unknown. SMC foam cells are reported to have a different phenotype than macrophage foam cells, due to relative deficiency of lysosomal acid lipase in SMC, and retain lipid droplets in their lysosomes rather than in the cytoplasm as seen in macrophages11. This would predict SMC foam cells have an altered gene expression pattern and phenotype different from macrophage foam cells. The literature is replete with discussions regarding the impact of SMC-derived foam cells to plaque, and indeed the role of such cells toward disease risk could be significant and the avenues for manipulation important for therapeutic considerations substantial, but this cell must be characterized with modern genetic and genomic methods. A recent comprehensive single cell study of carotid plaque has identified cells with gene expression features of both SMC and macrophage cellular phenotype, and further study of this cluster could address the current need in the field. 6 Without question, identification of the molecular pathways by which SMC lineage foam cells take up lipid and contribute to destabilizing the plaque could have significant importance for targeted therapy development. Moreover, additional multi-omic approaches incorporating DNA methylation, proteomic, metabolic and lipidomic data will allow us to better interrogate cellular physiology in the context of genomic analyses.\n3. The rigor of lineage tracing settings and cell annotations and clustering is appreciated and important for the field yet likely limits identification of cell states/transitions/cell-cell interactions that exist within the currently described populations. This has consequences – there may be critically important rare/uncommon cell states that are missed (e.g., SMC/FMC as pioneer or primed cells existing in normal media before and early into effect of diet on genetic background) and the more nuanced gene signatures in local cell populations that reflect metabolic or lipid stress may be missed. These issues should be acknowledged while also emphasizing the focus on rigor and reproducibility with stringent parameter settings.\nAgreed, we use stringent cluster methods based on a rigorous workflow that makes certain each cell cluster phenotype is homogeneous and different from other clusters in our analyses. We carefully identify the group of genes that are most representative of the variance across the dataset, and number of principle components that best describe the identified variation. Further, we check the gene ontology terms that are identified by cluster genes, to make sure there is minimal overlap between clusters. The field of single cell transcriptomics is rife with examples of over-clustering and claims of “new progenitor cell or intermediate phenotypes” associated with complex disease processes. But, as noted by the Reviewer, this rigor is sustained at the risk of missing interesting rare/uncommon cell states. In that regard, we have identified early markers of the FMC-1 and FMC-2 that identify these cells at early timepoints and transgenic Cre lines using these markers show differential expansion and migration of these true progenitor cell populations. We are currently mapping the cell state changes that are represented by these genes and their effect on SMC contribution to the plaque and fibrous cap. These studies are beyond the scope of the current work and will take considerable experiments to complete, but it is clear that this type of experiment will be highly instructive regarding the early developmental events for the FMC clusters. We adapt the text to clearly acknowledge the likelihood for the presence of such progenitor genes and their utility in understanding the molecular and cellular mechanisms of FMC-1 vs FMC-2 development and contribution to disease processes.\nWe further agree with the reviewer that -- more nuanced gene signatures in local cell populations that reflect metabolic or lipid stress may be missed. These issues should be acknowledged while also emphasizing the focus on rigor and reproducibility with stringent parameter settings.\nDiscussion, lines 665-674:\nOur study uses a rigorous workflow that identifies reproducible cell cluster phenotypes. However, this rigor is implemented at the risk of losing rare cell states including those that respond to metabolic or lipid stress and most notably, phenotypes that represent progenitor cells such as those that undergo clonal expansion in the disease setting. Future studies will be directed at identification of these putative precursor cells of interest utilizing novel lineage tracing methods (e.g. in vivo barcoding, somatic mutation detection, or novel lineage tracers using putative driver genes) that can better capture these rare populations while using our dataset as a source for validation. This way, descriptive characterization of novel cell phenotypes coupled with experiments using such lineage tracing will allow robust characterization of how these precursor cells contribute to the plaque and fibrous cap, and how the deletion of specific markers alter the course of disease.\n4. Mouse models cannot be relied upon to prove causality of genes in human CAD and CVD – they can be supportive. The statement in discussion on this should be modified. In this context, in the discussion, it would be very helpful for the reader if there was a brief discussion of the importance of complementary (to GWAS / regulatory etc) human genetic studies of rare coding variation in genes to infer causality and directionality in CAD/CVD of inferred and proved loss of function and gain of function mutations – and to note the state of this work related to the current paper e.g., for TCF21 and other proposed drivers and regulatory genes. Again, such data will be a critical complement to understand clear causality in humans. There is value also to acknowledging the need for studies and emerging data on somatic mutations in SMC-derived cells in lesions to understand specific causal SMC genes under the hypothesis of somatic variation.\nWe certainly agree that the mouse models cannot prove causality of specific genes and gene programs, but they can provide direct evidence for disease mechanisms by which human causal CAD genes mediate risk-related genetic and molecular processes in the vessel wall. We further agree that rare coding variation can provide highly useful information, occasionally providing a human disease phenotype and direction of effect for CAD genes that are identified with GWAS studies. A good recent example is the MFGE8 gene which has inframe mutations identified in the FinnGen biobank where these variants are present at high frequency and the inactivation of MFGE8 function has been correlated with greater lifespan.12 It is important to note of course that rare coding variation in the protective gene TCF21 would be expected to increase the risk for CAD, and indeed likely cardiovascular developmental defects that may be inconsistent with fetal survival. Also, TCF21 has important developmental roles in the kidney and lung, and mice lacking Tcf21 die at birth due to failure of lung function. This is the likely reason that we have not found coding region mutations in the TCF21 gene that are associated with protein function. Such mutations are likely removed from the genome through purifying selection.\nFurther, along the lines of using rare coding variation to inform on causality and direction of effect for candidate genes that reside in associated loci for complex human disease, a recently described approach from the Pritchard lab employs loss of function burden tests along with relevant Perturb-seq data to bridge the gap between genetic association and biological mechanism.13 By combining these two forms of data, their approach builds causal graphs in which the directional associations of genes with a trait can be explained by their regulatory effects on gene programs.\nThe possibility of using human somatic mutation in SMC to implicate genes of relevance to vascular disease is very intriguing. The NIH program Somatic Mosaicism Across Human Tissues (SMaHT) Network is focused on cataloging naturally occurring DNA mutations (somatic mosaicism) in healthy subjects to understand aging and disease. We have been in touch with these investigators to determine possible collaborations. While they are certainly interested, they are not actively studying any vascular tissue. Given the logistics of collecting and evaluating tissue from multiple humans and then pursuing a very deep sequencing effort, we have not embarked on these studies. Nevertheless, we are anxious to see such an effort move forward and will mention this need in the Discussion section of the manuscript.\nTo address the Reviewer’s comments we have made the following edits to the manuscript:\nDiscussion, lines 611-638:\nWe have undertaken these studies in order to better understand cellular and molecular mechanisms of CAD risk, by identifying and characterizing individual genes and gene programs that modulate SMC phenotype transitions. While such studies do not prove causality, they inform on possible mechanisms of disease risk that reside in the vessel wall. Although not explored in the work discussed here, there are a number of approaches related to mapping causality in human data that warrant further investigation, and can serve to strengthen the genetic data derived in our single cell studies of mouse models. For instance, more detailed study of rare coding variation, when present, can provide human disease phenotype and direction of effect for CAD genes that are identified with GWAS studies. A good example is the CAD GWAS MFGE8 (lactadherin) gene which was identified in the FinnGen biobank to have an association between an inframe insertion rs534125149 and protection against coronary atherosclerosis. 12 Along this line of reasoning, a recently described approach from the Pritchard lab employs loss of function burden tests along with relevant Perturb-seq data to bridge the gap between genetic association and biological mechanism.13 By combining these two forms of data, their approach builds causal graphs in which the directional associations of genes with a trait can be explained by their regulatory effects on gene programs. It is important to note of course that rare coding variation in the protective gene TCF21 would be expected to increase the risk for CAD, and indeed likely cardiovascular developmental defects that may be inconsistent with fetal survival. Also, TCF21 has important developmental roles in the kidney and lung, and mice lacking Tcf21 die at birth due to failure of lung function. This is one possible reason that we have not found coding region mutations in the TCF21 gene that are associated with protein function. Such mutations are likely removed from the genome through purifying selection.\nAnother approach to enriching the pool of CAD causal genes in SMC transitions is the possibility of using human somatic mutation to implicate genes of relevance to vascular disease. The NIH program Somatic Mosaicism Across Human Tissues (SMaHT) Network is focused on cataloging naturally occurring DNA mutations (somatic mosaicism) in healthy subjects to understand aging and disease. This network is not currently studying arterial tissues, and the tasks of collecting and subjecting tissue from multiple humans and then pursuing a very deep sequencing effort is daunting, but possibly well worth pursuing.\nReviewer #3 (Remarks to the Author):\nThis revised version is marked improvement upon the first submission of the manuscript as the authors have provided more evidence to support the identification of the SMC cell states, the role of Tcf21 as a top regulator of SMC transition to FMC/CMC, and TEAD1 as novel interacting partner of TCF21 at enhancer regions. The addition of spatial Xenium data from ApoE-/- atherosclerotic mice with SMC lineage tracing integrated with single-cell clusters and RNAScope analysis allows for orthogonal validation at high resolution of the previously identified SMC clusters. The authors have also included in vitro perturbation assays to validate their algorithmic findings on highly ranked transcription factors that affect the transition process. Knockdown of the selected TFs resulted in a more contractile SMC state. They also performed additional functional experiments on primary HCASMCs to validate that the FMC-2 population had greater probability to transition to the CMC state.\nWe thank the Reviewer for these positive comments.\nMajor comments:\n1. The study utilizes RNAScope, single-cell RNA-seq, spatial transcriptomics, but protein-level confirmation is still limited. The authors mention that individual antibody-based studies are time-consuming and difficult to scale. It would be helpful to at least stain for critical proteins of interest (ie. TCF21, TEAD1) in the lineage-traced atherosclerotic mouse sections (showing Myh11).\nThank you for this suggestion. Here, we stained for TEAD1 and TCF21 using aortic root serial sections from our Myh11 lineage-tracing mice at 16weeks high fat diet to show protein localization.\nWe trialed multiple conditions as below:\nIF Dilutions Antibody Primary Incubation Antigen\n(Overnight) retrieval\nTead1 1:50; Cell signaling D9X2L; 4C; room No\n1:100 Abcam 133533 temperature\nTcf21 1:50; Abcam 32981; 4C; room No\n1:100 Abcam182134; LS temperature\nBio B3504; Sigma\nHPA13189\nIHC\nTead1 1:100; Cell signaling -D9X2L; 4C; room Yes\n1:300 Abcam 133533 temperature\nTcf21 1:100; Abcam 32981; 4C; room Yes\n1:300 Abcam182134; LS temperature\nBio B3504; Sigma\nHPA013189\nWe were unable to obtain adequate images by immunofluorescence. With antigen retrieval, we were successful in obtaining signals using immunohistochemistry. Unfortunately, this results in loss of the lineage tracing marker. Images below are shown for Tead1 and Tcf21 in serial sections with minimal background as indicated by Rb-IgG control. We have incorporated this as a new supplemental figure (See below and S-Fig. 10).\nS-Fig 10. Immunohistochemistry demonstrating TEAD1, TCF21 protein localization with corresponding IgG control.\nMinor comments:\n1. A figure defining the key marker gene panel for identifying each cluster or top expressed genes would be important to show early. In figure 2, there are key drivers for FMC1 and FMC2 switching but identification markers would be helpful.\nThank you for pointing this out. We agree that early presentation of top marker genes are important. As such, we had included marker genes as an early supplemental figure (S-Fig 1L) in the original edit. We note that because of the high number of markers and space constraints, we felt it was a better fit in the supplement where readers can review both the dot plot (S-Fig 1L) as well as table (S-Table 2) in detail as multiple markers should be utilized to define these complex cell states with high fidelity. This was described early in lines 113-114.\n2. Please amend Figure 1A so it includes spatial transcriptomics using Xenium.\nWe have amended our Figure 1A as suggested.\n3. Having an understanding of the UMAP generated by the Xenium dataset would be important to show the clusters alongside their spatial distribution which is depicted.\nAs recommended, we generated corresponding spatial UMAPs with our transferred labels based on the transcriptomic data, demonstrating distinct clustering. We have incorporated this data as S-Fig 2B.\nReviewer References:"}, {"round": 3, "reviewer_comments": {"reviewer_2": "Authors have addressed all reviewer concerns and questions in full.\nThanks", "reviewer_3": "The authors have done well to answer previous comments and have performed TEAD1 and TCF21 staining using IHC on aortic root serial sections showing localisation. It is understandable that the lineage tracing marker is lost without the use of IF, which was not possible. I acknowledge that due to high number of markers and space constraints, it would be better to fit the full list in the supplemental table and indicate the key marker genes in the supplemental figure. The addition of the xenium platform and amendment to Xenium UMAPs have been well implemented to the figures of this version.\nThe authors have done a remarkable job at responding to all past comments and I have no further comments to add."}, "author_response": ""}]

Decision letter