Single-Cell Analysis Identifies NOTCH3-Mediated Interactions between Stromal Cells That Promote Microenvironment Remodeling and Invasion in Lung Adenocarcinoma

Lung adenocarcinoma (LUAD) is the most prevalent form of primary lung cancer in the United States, accounting for 40% of all lung cancer cases (1). The treatment of advanced LUAD has undergone a paradigm shift with the advent of immune checkpoint blockade (ICB), which inhibits immunosuppressive signaling pathways such as PD-1 (2). However, the majority of patients either do not respond or relapse (2). Transcriptomic studies using patient-derived samples have shown that the enrichment of a stromal or mesenchymal gene signature is strongly associated with resistance to ICB (3, 4). This suggests that mesenchymal cells, including cancer-associated fibroblasts (CAF), play a crucial role in mediating immunosuppression and treatment resistance in the tumor stroma.

Single-cell RNA sequencing (scRNA-seq) has revolutionized our understanding of tumor mesenchymal cell and CAF biology, providing valuable insights into CAF interactions with the immune system. In pancreatic cancer, which is characterized by a high degree of desmoplasia, two major CAF subsets have been identified through scRNA-seq studies: “myofibroblastic CAFs” (myCAF) with high expression of alpha smooth muscle actin (αSMA or ACTA2) and extracellular matrix (ECM) molecules, and “inflammatory CAFs” (iCAF) with high levels of cytokines and chemokines (5, 6). The selective depletion of LRRC15⁺ myCAFs has been shown to enhance the efficacy of ICB in a mouse model and is being explored clinically, indicating mesenchymal cells with a “myCAF” phenotype suppress tumor immunity (7). In LUAD, scRNA-seq has been leveraged to study tumor microenvironment (TME) heterogeneity, cell–cell interactions, and the immunophenotypes of tumor-infiltrating immune cells (8–10). However, previous scRNA-seq studies investigating mesenchymal cell or CAF biology in LUAD have been limited in terms of the number of profiled stromal cells, their relationship to previously defined “myCAF” and “iCAF” phenotypes, and their ability to capture cellular spatial architecture and associated signaling pathways.

To overcome these limitations, we harnessed scRNA-seq and imaging mass cytometry (IMC) to better decipher stromal heterogeneity and spatial localization. An in-depth investigation of the stromal interplay by ligand–receptor analysis revealed the activation of the NOTCH pathway in tumor-associated pericytes and CAFs driven by endothelium-derived NOTCH ligands. Leveraging publicly available patient sequencing data, we illustrated the ability of NOTCH3 in combination with several ICB response biomarkers to predict survival. Our study provides insights beyond previous scRNA-seq studies in understanding tumor stroma biology, NOTCH-driven intercellular cross-talk, and the clinical relevance of NOTCH3 as a potential stromal target.

All tissue samples from 9 patients for the scRNA-seq (cohort A) were obtained from Brigham and Women's Hospital, with written informed consent from the patients. The study was conducted in accordance with recognized ethical guidelines of the U.S. Common Rule and Belmont Report and was approved by the Institutional Review Board (IRB) protocol DFCI #98–063. The patient samples were sequenced using the OncoPanel platform (11) and the mutation status of 4 common LUAD tumor suppressor and oncogenes (TP53, EGFR, K-RAS, and MET) was provided in Supplementary Table S1. LUAD samples from a secondary cohort (B) of 7 patients for IMC were purchased from BioIVT, ISpecimen, and Discovery Life Sciences. These samples were acquired with written informed consent from the patients approved by the commercial manufacturers’ IRB protocols or under their partners’ IRB protocols, in accordance with recognized ethical guidelines of the U.S. Common Rule and Belmont Report. Please refer to the Supplementary Table S1 for more details.

Tumors (T) or adjacent nontumor (ANT) tissues from cohort A were minced with scissors and digested with a human tumor dissociation kit (Miltenyi Biotec, 130–095–929), as per the manufacturer's protocol with some modifications. For digesting lung samples for MRC001–004, the digestion cocktail from the kit was used, whereas the digestion cocktail with the addition of 0.8 mg/mL dispase II (Sigma, 4942078001) was used to digest MRC006–10. In brief, each lung sample in 5 mL of digestion cocktail buffer was put into a gentleMACS Octo dissociator (Miltenyi Biotec) for 30 minutes to 1 hour depending on the tumor size in a “37C_h_TDK_3 digestion” program. 20 mL of 10% FBS DMEM (Thermo Fisher Scientific) media was added into the digestion to stop the reaction and the material was further filtered through a 70-μm filter (Miltenyi Biotec). The remaining clumps on filters were grounded using the bottom of a 1-mL syringe (BD Biosciences). Cell suspensions were spun down, lysed using ACK buffer (Thermo Fisher Scientific), and spun down again to obtain single-cell suspensions ready for downstream experiments.

To prepare viable cells for scRNA-seq, single-cell suspensions from MRC001–003 were stained with a live/dead dye (Thermo Fisher Scientific, L34975) over ice for 15 minutes in the dark. Single cells from donor MRC004, 006–010 were stained with a cocktail of live/dead dye, an Fc blocker (BioLegend, 422302), and an anti-CD45 antibody (BD Biosciences, 560976) over ice for 30 minutes before proceeding to cell sorting. Viable cells or viable CD45⁺/CD45⁻ cells were sorted out using a FACSAria III cell sorter (BD Biosciences).

Sorted viable cells were processed by Chromium Single Cell 3′ Reagent v3 kits (10X Genomics, 120234) to generate single-cell cDNAs and prepare barcoded libraries, as per the manufacturer's protocol. Sorted live cells were suspended into a cell density of around 1 × 10⁶ cells/mL in 0.1% BSA of PBS. To minimize doublet percentages, 4,000 cells for donor MRC001–003 or 6,000 cells for the remaining donors were loaded into each lane of a 10× chip. Cells were then partitioned into single-cell gel beads in emulsions (GEM) inside the Chromium instrument, where full-length cDNA synthesis occurred. Cleaned-up cDNAs were then amplified, fragmented, and attached with 5′-adaptor and sample index. Libraries were sequenced using a 150-bp paired-end configuration.

Tumor samples were fixed in 10% formalin, stored in 70% ethanol before paraffin embedding, and cut into 5-μm sections. IHC staining was performed as previously described (12). Slides were deparaffinized, rehydrated, boiled with DC NxGen (Biocare Medical), incubated in Peroxidazed 1 (Biocare Medical, PX968), and blocked with Background Punisher (Biocare Medical, IP974G20) before staining with primary antibodies: anti-TTF1 (Abcam, ab133638), and anti-P63 (Biocare Medical, CM 163A). Sections for IF staining were deparaffinized, rehydrated, and boiled as described above. Sections were then stained with primary antibodies: NOTCH3 (Abcam, ab23426) and MCAM (Abcam, ab210072).

Six-micron–thick tissue sections from cohort B were stained with 18 antibodies (Supplementary Table S2), which include markers covering mesenchymal cells (FAP, MCAM, ACTA2, CD90, and collagen I), ECs (CD31), tumor/epithelial cells (E-cadherin), nonepithelial cells (Vimentin), and immune cells (CD3, CD4, CD8, CD68, CD14, CD33, CD15, CD16, CD11C, and CD20). CD14, CD11C, and CD33 were removed for further analysis due to poor staining quality. Hyperion Imaging System (Fluidigm) was utilized to acquire a randomly selected tumor area on slides for each LUAD section. The process and analysis of IMC images were followed as previously reported (13). In brief, using Fiji software, IMC mcd. Files were transformed into stack tiff files, which were then imported into Ilastik (Version 1.3.2) to produce cell segmentation masks followed by generation of probability masks using CellProfiler (Version 3.1.5). HistoCAT (version 1.75) used all files created above to produce t-SNE plots for projecting high-dimensional single-cell data into two dimensions to enable visualization of marker-specific cell types as well as a Phenograph to define complex phenotypes shared across tumors based on the staining intensity of tested cell markers. The phenograph was used to perform neighborhood analysis with histoCAT default settings, in which a permutation test to compare the number of interactions between all cell types in a given image to that of a matched control containing randomized cell phenotypes was used to determine the P value.

Mesenchymal cells were expanded from freshly resected tumor samples obtained from MRC002, -003, and -004 in cohort A, following a previously described protocol (12). The tumors were minced and digested into approximately 1-mm³ fragments. These fragments were then placed in 6-well tissue culture plates containing DMEM supplemented with 10% heat-inactivated (HI) FBS. Mesenchymal cells were allowed to extravasate from tissue fragments and were expanded until reaching confluency before passaging. Cells were stained with antibodies: NOTCH3 (BD, 745463) and MCAM (BioLegend, 361004) for flow cytometry analysis. These cells were used in subsequent assays within 10 passages. D4A1 cancer-associated fibroblasts were purchased from Bio IVT. The cells were derived from a stage II lung adenocarcinoma patient. The donor number associated with D4A1 fibroblasts is 426674A1, and the lot number is DT01086P1. Human pulmonary artery ECs were purchased from Lonza and cultured in the conditioned media per Lonza's instruction. NCI-H1299 cells were obtained from ATCC and cultured in F12/K and RPMI1640 (Gibco, Thermo Fisher Scientific) supplemented with 10% HI FBS. The cell line was authenticated by short-tandem repeat profiling and tested negative for Mycoplasma. NCI-H1299 cells were infected with Incucyte Nuclight Lentivirus (Sartorius) and selected with puromycin to establish a mKate2⁺ stable cell line for the invasion assay. All cells were maintained at 37°C and 5% CO₂ in an incubator. For compound treatment, MRK-003 (Merck & Co., Inc.; ref. 14), a γ-secretase inhibitor, was used.

For MRK-003 treatment, 1 × 10⁵ ECs and 1 × 10⁵ mesenchymal cells were cocultured and seeded into one well of a 6-well plate. The endothelial conditioned medium was used for the coculture. On the next day, cells were changed to fresh media containing DMSO or 1 μmol/L MRK-003. For siRNA transfection, 2 × 10⁵ ECs and 2 × 10⁵ mesenchymal cells were cocultured in one well of a 6-well plate. On the next day, 200 pmol of siNOTCH3_1 (Thermo Fisher Scientific, HSS107256), siNOTCH3_2 (Thermo Fisher Scientific, HSS107254), or nontargeting control (siNT) was transfected into each well using Lipofectamine RNAi Max (Thermo Fisher Scientific). RNA and supernatants were collected 3 days after treatment for subsequent analysis. RNAs were extracted using the RNeasy Plus Mini Kit (Qiagen), and 1 μg of RNA was used for cDNA synthesis using iScript cDNA Synthesis Kit (Bio-Rad). TaqMan Human Extracellular Matrix & Adhesion Molecules Arrays (Thermo Fisher Scientific) and following TaqMan primers were used: HEYL (Hs01113778_m1), HEY1 (Hs05047713_s1), HES1 (Hs00172878_m1), HES4 (Hs00368353_g1), COL1A1 (Hs00164004_m1), RPL30 (Hs00265497_m1), RPLP0 (Hs00420895_gH). QuantStudio 6 Flex Real-Time PCR System (Thermo Fisher Scientific) was used to perform RT-qPCR. ΔΔC_t was calculated by normalizing treated groups to corresponding controls, and 2^(–ΔΔC_t) was the relative fold change.

CRISPR gene editing was conducted using the Alt-R CRISPR-Cas9 System (IDT). RNP complexes were prepared by combining 2 μL of Alt-R S.p. HiFi Cas9 Nuclease V3 (IDT, 1081060) with 3 μL of 200 μmol/L customized sgRNAs, and subsequently mixed with 1 × 10⁵ D4A1 mesenchymal cells suspended in 20 μL of P3 Primary Cell 4D-Nucleofector X Kit solution (Lonza, V4XP-3032). The nucleofection procedure was performed utilizing the CM-138 program with the 4D-nucleofector Core Unit (Lonza, AAF-1003B). The NOTCH3 target sequences were: sgRNA1 (GCCACTATGTGAGAACCCCG) and sgRNA2 (AGGGTGCACAGGGCACCGCG). The sequence of nontargeting sgRNA (NT_sgRNA) was CGTTAATCGCGTATAATACG.

To determine the cell surface NOTCH3 expression, cells were detached using a cell scraper, stained with a live/dead dye and anti-NOTCH3 antibody (clone MHN3–21) for more than 30 minutes, and analyzed by flow cytometry. Data analysis was performed using FlowJo (Version 10). Cells were first gated on the basis of forward (FSC-A) and side (SSC-A) scatters to exclude debris. Single cells were then selected on the basis of SSC-A versus SSC-W parameters. Dead cells were excluded on the basis of the positive staining of the live/dead dye. The positive cell-surface staining in gated live cells was determined by comparing them to fluorescence minus one as a negative control.

The ability of cells to invade the surrounding matrix was assessed based on the IncuCyte S3 3D Spheroid Invasion Assay (Sartorius). Briefly, 5,000 cells derived from MRC002, 003, or 004 were added into one well of an ultra-low attachment plate (S-Bio) and treated with DMSO or 10 μmol/L MRK-003 for 3 days to form spheroids at 37°C and 5% CO₂ in an incubator. 50% GFR Matrigel (Corning) containing DMSO or MRK-003 was added on top of the spheroids. For D4A1 alone or D4A1 and mKate2⁺ H1299 (1:1) cocultures, a total of 2,000 cells were added into one well of an ultra-low attachment plate (Nexcelom Bioscience) and treated with DMSO or 10 μmol/L MRK-003 as describe above to allow spheroid formation. Supernatants were collected for assessing COL1A1 production. For CRISPR-edited D4A1 cell, or the coculturing mKate2⁺ H1299 with nontargeting control (NT_sgRNA) or NOTCH3 knockout (NOTCH3_sgRNA1, and NOTCH3_sgRNA2) D4A1 cells at a 1:1 ratio, a total of 1,000 cells were plated for 3 days to form spheroids. For the treatment of the additional γ-secretase inhibitors DAPT (Tocris Bioscience; ref. 15) and MRK560 (Tocris Bioscience; ref. 16), 1,000 cells without drug treatment were first added into each well of an ultra-low attachment plate to form spheroids for 3 days. GFR Basement Membrane Matrix (Gibco, Thermo Fisher Scientific) was added onto spheroids in each well. Plates were subsequently centrifuged at 300 × g for 5 minutes at 4°C and incubated for 30 minutes at 37°C to promote polymerization. Culture medium containing DMSO or γ-secretase inhibitors at a concentration of 10 μmol/L was added into each well post polymerization and the plates were incubated and monitored by IncuCyte S3 placed in an incubator at 37°C with 5% CO₂ for indicated hours. Cells that invaded the surrounding matrix were observed using a phase-contrast inverted microscope (magnification, ×4) and images were captured. Data were analyzed using IncuCyte S3 Spheroid Software Module (version 2021A).

The concentration of type 1 collagen in supernatants collected above was measured using the human COL1A1 AlphaLISA Detection Kit (PerkinElmer). Standards or samples were added into white Opaque 96- or 384-well microplates (PerkinElmer) and then biotinylated antibodies and beads were added according to the kit instruction. The results were measured using EnVision.

In the cell invasion assay, the cell viability of spheroids derived from MRC002, MRC003, and MRC004 was assessed using the CellTiter-Glo 3D Cell Viability kit (Promega). In brief, the plate and reagents were equilibrated to room temperature for 30 minutes. The CellTiter-Glo 3D reagent was added to each well of a 96-well plate containing the spheroids. The plate was then shaken for 5 minutes to induce lysis of the cells. The plate was incubated at room temperature for additional 30 minutes before measuring the luminescence signal with Envision multilabel reader (PerkinElmer).

Incucyte Caspase-3/7 Dyes at a final concentration of 5 μmol/L were used to assess apoptotic effect on H1299-mKate2 expressing cells or D4A1 mesenchymal cells following the compound treatments using IncuCyte S3. Cells were imaged using both phase contrast and green fluorescence channels with a magnification of 10x. Data were analyzed using IncuCyte S3 Cell-by-Cell Analysis Software Module (version 2021).

scRNA-seq reads were mapped to the human genome (GRCh38) using CellRanger v1.1.0 and then processed through Seurat v3.0 (17). Cells with mitochondria percentage over 10% or extreme unique gene totals (less than 500 or over 5,000 per cell) were removed from the analysis. The data after QC inspection was then normalized using sctransform v0.3.2 (18), and all libraries were merged after batch effect removal through Harmony v1.0 (19). Following the standard protocol provided by Seurat authors, principal components were computed and used for the Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction. Cell clusters were identified at resolution 0.3 and annotated on the basis of prior knowledge. Of note, one cluster, accounting for 0.9% of the total cells, was identified as the proliferating cell cluster. This cluster showed high expression of the proliferation consensus signature (3) and cell-cycle regulation genes such as STMN1. Because of the low percentage of proliferating cells within the total cell population, we did not perform cell-cycle correction. Two clusters with low cell count were expressing markers of more than one primary cell type and were removed from downstream analysis.

Mesenchymal, lymphatic, and vascular endothelial clusters were annotated on the basis of canonical markers such as COL1A2, PROX1, and RAMP2, respectively, for further dimensionality reduction using the FindClusters function in the Seurat package. The differentially expressed genes of mesenchymal cells compared with all other cells were listed in Supplementary Table S3. Mesenchymal cells were further subclustered with resolution of 0.3 and annotated as four subgroups with FAP, PDPN, MCAM, and ACTA2 expression levels based on previous reports (6, 20, 21). Vascular endothelial cells were subclustered with resolution of 0.3. Selected markers for EC subclusters’ annotation were curated based on a previous report (22). The marker gene lists for mesenchymal cells and endothelial cells were generated using the FindAllMarkers function in Seurat, and the results are provided in Supplementary Table S3 and Table S4, respectively.

UMAPs, violin plots and bubble plots visualizing expression of genes and markers were generated using Seurat 4.0.6 and ggplot2 v3.3.5, and the color palettes were loaded from ggsci v2.9. Unless otherwise specified, statistical testing in differential expression analysis was the Wilcoxon rank-sum test, and the P values for multiple comparisons were adjusted using the Benjamini–Hochberg method. Heat maps were generated using pheatmap R package v1.0.12.

Subpopulation abundance comparison between T and ANT: The relative abundance of each subpopulation over total mesenchymal or endothelial cell counts was calculated in T and ANT samples, respectively. The Wilcoxon rank-sum test was applied to determine the significance of relative abundance differences between T and ANT in each subpopulation. The figures were styled using ggprism v1.0.3.

Cell–cell communications between mesenchymal and endothelial subtypes were inferred based on the analysis of differential expressions of known ligand–receptor pairs between T and ANT samples, which was accomplished by following the official workflow of CellChat v1.1.3 (23). The gene expression of endothelial cells and mesenchymal cells after sctransform processing were grouped by T and ANT and loaded into CellChat separately. These two objects went through preprocessing with the following functions using standard parameters: “identifyOverExpressedGenes,” “identifyOverExpressedInteractions,” and “projectData.” Then the communication probabilities of T and ANT were analyzed separately using the core functions “computeCommunProb,” “computeCommunProbPathway,” and “aggregateNet” with the standard parameters and merged into one object for T vs. ANT comparison. Figures were generated using CellChat functions: “netVisual_diffInteraction” for overall pathway enrichment in T versus ANT, “netVisual_aggregate” for circle plots, “netAnalysis_signalingRole_heatmap” for determining signal senders and receivers of top pathways in mesenchymal–endothelial interactions, and “netVisual_chord_gene” for chord plot of NOTCH pathway information flow from endothelial cells to mesenchymal cells.

To eliminate the potential bias from one donor, we applied a leave-one-out strategy to the differential expression analysis between T and ANT in endothelial and mesenchymal subtypes. One sample was removed in each iteration, and Wilcoxon rank-sum tests were performed between T and ANT to determine the significance of expression enrichment in T. The median and IQR of the P values from the bootstrap analysis were collected for visualization. For expression of each gene by a specific cell type, if the Benjamini–Hochberg P value was over 0.05 in the all donor–included iteration, the data would be filtered out.

Expression correlation analysis was performed in LUAD tumor samples only, and Spearman correlation between the average expression of signature genes and scRNA-seq–derived fibroblast subset fingerprints was calculated. Differential expression analysis with Benjamini–Hochberg P_adj value was performed in TCGA LUAD samples between T and ANT, and the ROC-AUC was calculated by the performance of classifying T using the expression of each tested gene. The Gene set pathway enrichment analysis was performed using the combination of Kyoto Encyclopedia of Genes and Genomes (KEGG; ref. 24) and GeneGo MetaCore Pathways. The above statistical analysis, and survival analysis, including the Cox proportional hazard model and Kaplan–Meier survival analysis, were performed in MATLAB 2021a.

The consensus signatures, including stroma/EMT/TGFβ, angiogenesis, and the 18-gene T-cell inflamed GEP, were published previously (3). The smooth muscle signature was obtained from the Human Protein Atlas (HPA), including 34 genes in the core smooth muscle cell transcriptome that were predicted as enriched in smooth muscle cells in 50%–75% profiled tissues (25). DLL4 and JAG1 stimulation marker signatures were obtained from a publication (26). Fingerprint genes of iCAFs, myCAFs, and antigen-presenting CAFs (apCAF) in Elyada and colleagues were obtained from the publication (5). Fingerprint genes of apCAFs in Kerdidani and colleagues were antigen presenting MHCII⁺ fibroblast-enriched genes collected from the differential expression analysis in the study (27). The gene lists are provided in Supplementary Table S5.

Statistical analysis for scRNA-seq and bulk RNA-seq was performed using the Wilcoxon rank-sum test. Statistical tests for functional assays were conducted using Student t test. Additional statistical tests are indicated in the legends or corresponding method sections for specific details. Differences were considered statistically significant when P ≤ 0.05.

The scRNA-seq data generated and analyzed in this study are publicly available in the NCBI BioProject database with the accession number PRJNA1055415, in Gene Expression Omnibus (GEO) at GSE253013, and within the article and its supplementary data files. Bulk RNA-seq databases from the collaboration between Merck & Co., Inc. and H. Lee Moffitt Cancer Center (the “Collaboration”) and The Cancer Genome Atlas (TCGA) were used in this study. The Collaboration dataset contains 20,000 tumor samples and represents more than 25 different cancers, including 1,434 LUAD samples (3). Patient samples were obtained by the H. Lee Moffitt Cancer Center, and the gene expression data was generated using HuRSTA-2a520709 GeneChips (Affymetrix). Additional details regarding the Collaboration dataset and the assay platform (GPL10379) have been published (3) and can be accessed in GEO at GPL10379. TCGA data were obtained from the Genomic Data Commons Data Portal, and information is available on the TCGA website. Raw data was processed as previously described (3). All other raw data are available upon request from the corresponding author.

We performed scRNA-seq on a cohort of 9 treatment-naïve patients, of which, 6 had paired adjacent nontumor (ANT) tissues (Fig. 1A). To obtain a representative number of nonimmune and immune cells, we sorted live CD45⁻ (nonimmune) and CD45⁺ immune cells from 6 out of 9 freshly resected LUAD tumors prior to scRNA-seq. Aiming to match the current standard of fibroblast/stromal-centric scRNA-seq studies, we obtained 256,379 total cells and 13,857 mesenchymal cells, respectively, after quality control and doublet removal. The cells were visualized as a dimensional reduction plot using UMAP (Fig. 1B). Eight distinct cell types were identified by using RNA expression of canonical markers (Fig. 1B and C). In-depth investigation of nonstromal cell types identified diverse cell subclusters in the TME that are similar to previous scRNA-seq studies in non–small cell lung cancer (NSCLC; Supplementary Fig. S1A–S1D; 8, 10, 28). We next sought to decipher mesenchymal cell heterogeneity and its intercellular network to interrogate dysregulated pathways in the tumor stroma.

To characterize the heterogeneity of mesenchymal cells in LUAD, we first subclustered them and observed four distinct clusters: FAP⁺PDPN⁺ CAFs (5, 6, 20, 29), MCAM⁺ACTA2⁺ pericytes (21), ACTA2⁺ smooth muscle cells (SMC), and lung resident fibroblasts (LRF; Fig. 2A and B). LRFs were predominantly derived from ANT tissues (Fig. 2C; P < 0.001), and SMCs comprised of cells from both tumors and ANTs. Remarkably, MCAM⁺ACTA2⁺ pericytes and FAP⁺PDPN⁺ CAFs were almost exclusively originated from tumors (P = 0.002, for both subtypes). We confirmed the presence of FAP⁺PDPN⁺ CAFs and MCAM⁺ pericytes by flow cytometry (Fig. 2D). Our annotation of the four clusters was further confirmed by investigating additional marker gene sets expressed by these cells (Fig. 2E; Supplementary Fig. S2A–S2C and Supplementary Table S3). For instance, pericytes expressed the well-established marker gene PDGFRB. SMCs were positive for previously described genes that have functional impacts on vascular SMCs, such as MFAP4. In contrast to other clusters, FAP⁺PDPN⁺ CAFs showed a gene expression spectrum characterized by a high expression of collagens and matrix-degrading enzymes, indicating these cells may have a myofibroblast phenotype of promoting ECM deposition and remodeling (Supplementary Table S3).

We next explored the phenotypes of mesenchymal cell clusters by correlating the marker gene sets of the four clusters with the consensus bulk RNA-seq signatures derived from much larger lung adenocarcinoma cohorts, including the Collaboration dataset, TCGA, and the Human Protein Atlas (3, 25). The stroma/EMT/TGFβ and angiogenesis signatures represent canonical biological processes that are common across different tumor types and are associated with resistance to ICB monotherapy (3). The marker gene sets of both CAFs and pericytes were strongly correlated with the Stroma/EMT/TGFβ signature compared to the other two clusters, whereas pericyte marker genes were closely related to the angiogenesis signature (Fig. 2F). Consistent with our annotation, the marker gene set of SMCs showed the highest correlation with the muscle signature. These observations were further confirmed in the TCGA LUAD patient cohort (Supplementary Fig. S2D). Given that the stroma/EMT/TGFβ signature is associated with a “myCAF” phenotype (6, 7), we further explored the previously defined “myCAF” and “iCAF” phenotypes by mapping the related gene signatures to the four clusters. The myCAF marker genes were most enriched in CAFs and then in pericytes, whereas iCAF marker genes were enriched in LRFs (Fig. 2G). Similarly, a gene signature of apCAFs, derived from a prior NSCLC study (27), was also enriched in LRFs and showed the highest expression in myeloid cells (Supplementary Fig. S2E and S2F). Given that CAFs are characterized by their diverse phenotypes and heterogeneity, we compared the expression levels of “myCAF” and “iCAF” signatures within the CAF subclusters (Supplementary Fig. S2G). CAF-S3 exhibited the highest expression of the “myCAF” signature (Supplementary Fig. S2H). On the other hand, although CAF-S7 showed the highest expression of the “iCAF” signature among all CAF subclusters, its expression was lower compared to the levels observed in LRFs (Supplementary Fig. S2H and S2I). These findings indicate the existence of heterogeneity within the CAF subclusters and suggest that the CAF cluster as a whole adopts a myCAF-like phenotype in this cohort of LUAD. Notably, this cell cluster also has the highest correlation with stromal/EMT/TGFβ bulk RNA-seq signature, which has demonstrated clinical association with ICB resistance (3).

Next, we sought to understand how the spatial organization of cell types dictates intercellular interactions that create and maintain a reactive tumor stroma. To dissect cell neighbors interacting with the three major mesenchymal subpopulations: FAP⁺PDPN⁺ CAFs, MCAM⁺ACTA2⁺ pericytes, ACTA2⁺ SMCs, we exploited a high-dimensional imaging technique known as IMC and multiplexed 18 markers in an independent cohort of 7 histologically confirmed LUAD samples (Fig. 3A and B). We found that MCAM⁺ and ACTA2⁺ cells comprised a thin layer surrounding CD31⁺ blood vessels and were indicative of MCAM⁺ACTA2⁺ pericytes or ACTA2⁺ vascular SMCs identified in the scRNA-seq data. In contrast, FAP⁺ collagen I⁺ CAFs resided in stromal regions between vascular zones and E-cadherin⁺ tumor parenchyma. Cell segmentation masks and t-SNE plots for each marker were generated from 7 samples for further analysis and visualization (Fig. 3C; Supplementary Fig. S3A).

To enable quantitative spatial interrogation of cell-to-cell interactions, we used PhenoGraph (30) and identified 25 clusters characterized by specific epitopes across 7 donors (Fig. 3D). The annotation of each cluster can be found in Supplementary Table S6. Clusters 8 and 12 were FAP⁺ CAFs with uneven expression of collagen I and CD90. We identified a CD31⁺MCAM⁺ population (cluster 7) that likely contained both CD31⁺ ECs and surrounding MCAM⁺ pericytes. The mixture of the lineage markers was likely due to the low resolution of IMC to separate ECs and surrounding MCAM⁺ pericytes. This result indicates the geographic proximity of CD31⁺ and MCAM⁺ cells. Cluster 16 was a mixture of SMCs and pericytes with ACTA2 and MCAM expression, whereas cluster 18 expressed MCAM as well as additional mesenchymal cell markers CD90 and Vimentin. We then performed neighborhood analysis to investigate significant interactions or avoidances of cell–cell neighbors in tumors (Fig. 3E; Supplementary Fig. S3B). In more than 50% of tested images, the FAP⁺ CAF cluster 8 significantly interacted with CAFs (cluster 8 and 12), CD4⁺ T cells (cluster 2), CD8⁺ T cells (cluster 10), pericytes, and vSMCs (cluster 16, 18, and 7), and ECs (cluster 6 and 7; Fig. 3E, P < 0.01). Similarly, pericytes and vSMCs (cluster 16 and 7) had robust intercellular interactions with FAP⁺ CAFs and T cells in 50% of tested images. Overall, the neighborhood analysis reveals a close topographical association between tumor-associated mesenchymal cells and endothelial/T cells.

Previous studies have focused on delineating the interactions between mesenchymal and T cells to sculpt an immunosuppressive TME (4, 31, 32). Although ECs are emerging as a key player in pathogenic stromal remodeling in noncancerous indications (26, 33, 34), little is known about how ECs cooperate with mesenchymal cells to remodel and activate tumor stroma in LUAD. Thus, we investigated the endothelial-mesenchymal cross-talk and dissected the underlying molecular determinants. We identified both vascular and lymphatic ECs. The lymphatic ECs formed a unique cluster in the UMAP containing all primary cell types (Fig. 4A), and the vascular ECs were reclustered into 7 subclusters. Subcluster N1, N2, and N3 were significantly enriched in ANT samples, while subcluster T1, T2, and T3 were enriched in tumors (Fig. 4B). Similar to previous findings (22), Endo T1 was an activated postcapillary vein EC cluster expressing ACKR1 and high-endothelial venule EC marker IGFBP7; Endo T2 and T3 were PLVAP^high immature neovascular ECs; endo G1 expressed several arterial EC markers, including DKK2, GJA5, FBLN5, and SERPINE2, and thus was an arterial EC cluster (Fig. 4C; Supplementary Table S5).

To understand intercellular communications based on the differential expression of ligand–receptor pairs between mesenchymal and ECs, we performed statistical analysis using CellChat, with a focus on predicted cell–cell contact pathways that were specifically enriched in tumors (23). We prioritized 11 signaling pathways with the highest statistical measurements of information flow in tumors versus ANTs (Fig. 4D). The top three pathways were ESAM, MHC-II and NOTCH signaling, which were also among the top 5 when considering all potential ligand–receptor interactions in the CellChatDB (Supplementary Fig. S4A). Some of the 11 pathways, such as ESAM, SEMA6, and ADGRE5, were predominantly driven by predicted intracellular interactions within ECs or mesenchymal cells (Fig. 4E; Supplementary Fig. S4B). To understand how the predicted interactions between cell types may modulate the phenotype of mesenchymal cells, we filtered for predictions of ligands from ECs and receptors from mesenchymal cells. The top two pathways that met the requirements were MHC-II and NOTCH (Fig. 4E; Supplementary Fig. S4B and S4C).

Given the established involvement of NOTCH signaling in fibrosis (26, 33), we sought to determine its prevalence in the tumor stroma of LUAD. Multiple endothelial subtypes contributed to the ligand expression, whereas NOTCH3 was the major receptor expressed mainly by pericytes, and by CAFs at a lower level (Fig. 4F; Supplementary Fig. S4D). To validate that the significant expression of NOTCH pathway ligand–receptor pairs in tumors is not determined by a single donor, we performed a leave-one-donor-out analysis. The distributions of P values of ligand–receptor pairs in each subtype were visualized (Fig. 4G and H). Given that the number of pericytes in ANT was too low to perform statistical analysis (n < 10; Fig. 2C), we compared the enrichment of NOTCH receptor genes in all mesenchymal cells and other three subtypes and found the significant enrichment of NOTCH3 across donors (Fig. 4G). Endo G1 had significant enrichment of NOTCH ligands DLL4, JAG1, and JAG2 across donors, whereas the increased ligand expression in other subtypes had a larger donor-driven variation (Fig. 4H). We further confirmed the specific NOTCH3 expression on MCAM⁺ pericytes within the TME by immunofluorescence staining (Fig. 4I). Given that there are limited donors in our scRNA-seq dataset, we evaluated the differential NOTCH receptor expression in a larger cohort utilizing the TCGA dataset, which includes 530 LUAD T samples and 59 LUAD ANT samples. NOTCH3 is the only NOTCH receptor that significantly upregulated in tumors (Fig. 4J).

To further explore the potential downstream effects of NOTCH dysregulation in LUAD tumor samples, we evaluated the activation of this pathway by cross-referencing the NOTCH activation scores developed by Wei and colleagues, which was generated by in vitro JAG1 and DLL4 stimulation of synovial fibroblasts (26), and the canonical NOTCH response markers in our scRNA-seq dataset (35, 36). NOTCH activation scores and the canonical NOTCH-responsive markers were significantly enriched in all mesenchymal cells from tumor samples (Supplementary Fig. S4E). Similarly, by investigating mesenchymal subtypes, CAFs and pericytes had significant enrichment of both activation scores compared with LRFs (Supplementary Fig. S4F). Our observation indicates that the increased expression of NOTCH ligand–receptor pairs translates into NOTCH signaling activation in mesenchymal cells in tumors.

To interrogate the role of NOTCH signaling in regulating the function of tumor-associated mesenchymal cells, we expanded stromal cells from freshly-resected tumor samples and were able to propagate MCAM⁺NOTCH3⁺ cells from three independent donors as indicated by flow staining (Fig. 5A). We then cocultured these cells with primary human pulmonary artery ECs that served as the source of NOTCH ligands. Subsequently, the coculture was treated with MRK003, a well-characterized γ-secretase inhibitor known to block the NOTCH pathway and function as a pan-NOTCH inhibitor (14, 37). MRK003 treatment robustly reduced the expression of several downstream NOTCH target genes in the mesenchymal and EC cocultures (Fig. 5B). After confirming the inhibition of NOTCH pathway, we further investigated whether MRK003 treatment affects the expression of genes associated with ECM remodeling and cell adhesion in the coculture. We found significant downregulation of several fibrillar collagen genes, including COL1A1, as well as the fibrotic gene CTGF (Fig. 5C). We also confirmed that MRK003 treatment significantly reduced COL1A1 at protein levels (Fig. 5D). To further investigate the role of NOTCH3 in regulating collagen expression, we utilized two siRNAs to downregulate NOTCH3 mRNA levels by 70% in the EC and mesenchymal cell cocultures. Consequently, we observed reduced cell surface NOTCH3 expression, accompanied by a significant reduction in COL1A1 expression (Supplementary Fig. S5A-S5E).

Given that enhanced collagen deposition promotes cell invasion in the TME (38, 39), we further explored how NOTCH inhibition affects the invasiveness of mesenchymal cells. Cells were embedded in three-dimensional matrix in the presence of DMSO or MRK003, and the cell invading area was quantified. MRK003-treated cells showed a reduction in the invading area compared with control cells in all donors tested, with significant inhibition seen in two of the three donors (Fig. 5E and F). We further tested a high NOTCH3-expressing mesenchymal cell line derived from an independent donor (D4A1) to confirm the inhibitory effect of MRK003 on cell invasion (Fig. 5G and H). We tested two additional γ-secretase inhibitors, DAPT (15), and MRK560 (16), that block NOTCH pathway signaling. MRK560 significantly reduced mesenchymal cell invasion, whereas DAPT treatment showed a trend of suppressing cell invasion that did not reach statistical significance (Supplementary Fig. S6A). Because previous reports show that tumor-associated mesenchymal cells promote tumor cell invasion and metastasis (40, 41), we investigated the potential of inhibiting the NOTCH pathway to suppress their proinvasive effect. We cocultured a fluorescence-labeled metastatic NSCLC line H1299 with or without the D4A1 mesenchymal cells in the matrix and tracked the tumor cell invasion area reflected by fluorescence intensity. The presence of D4A1 mesenchymal cells increased the invasive capacity of H1299, and MRK003 treatment robustly suppressed the invasion of H1299 cocultured with D4A1 cells to the same level as H1299 alone (Fig. 5I). We also observed that MRK003 treatment directly suppressed H1299 invasion (Fig. 5I).

To determine whether the impaired invasion is associated with changes in cell growth or ECM production, we quantified viable cells or monitored cell death in mesenchymal cells over time. We did not observe any significant changes in the number of viable cells or the occurrence of a cell apoptosis phenotype, as indicated by caspase 3 and 7 signals (Supplementary Fig. S6B–S6D). However, MRK003 treatment did significantly change COL1A1 production in mesenchymal cells (Supplementary Fig. S6E). In contrast, MRK003 treatment in H1299 cells significantly increased caspase 3 and 7 signals, indicating the induction of cytotoxic effects (Supplementary Fig. S6F). The cytotoxic effects led to a decrease in cell numbers, as evident by the reduced red signal intensity of the H1299 spheroids (Supplementary Fig. S6G). This data suggests that MRK003 suppresses the invasive capacity of H1299 by inducing cell death. To further understand whether mesenchymal cells can promote tumor cell invasion in a NOTCH3-dependent manner, we utilized CRISPR gene editing to generate NOTCH3 knockout D4A1 mesenchymal cells and found NOTCH3 knockout cells exhibited reduced invasive capacity (Supplementary Fig. S6H and S6I). Subsequently, we cocultured H1299 cells with mesenchymal cells that were either NOTCH3 wild-type or knockout and observed reduced invasion in H1299 cells cocultured with NOTCH3 knockout mesenchymal cells (Supplementary Fig. S6J). The data suggest that mesenchymal cells can promote tumor invasion through a NOTCH3-mediated effect.

To gain better insight into the function of NOTCH3 in tumor cells, we assessed its expression in H1299 cells. We found limited cell surface expression of NOTCH3 in H1299 cells (Supplementary Fig. S6K). In addition, by analyzing our scRNA-seq dataset, we observed minimal NOTCH3 expression in epithelial cells (Supplementary Fig. S6L and S6M). On the basis of the limited expression of NOTCH3 in tumor cells, we did not further investigate its role in regulating cell invasion of H1299 cells. The observed effect of MRK003 treatment could potentially be attributed to the presence of NOTCH1 expression in H1299 cells (42). Consequently, MRK003 may also inhibit NOTCH1 signaling in H1299 cells, leading to reduced cell death and invasion. In summary, our in vitro data demonstrates the functional role of NOTCH signaling in promoting mesenchymal cell invasion, potentially through the regulation of ECM molecule production, such as COL1A1.

Our in vitro functional data suggest that the activation of NOTCH pathway can regulate collagen production and ECM remodeling. To explore the clinical relevance of this observation, we evaluated the association of NOTCH3 expression with global gene expression in two independent LUAD patient cohorts (Fig. 6A; Supplementary Table S7). COL1A1 was among the top genes with the highest correlation with NOTCH3. We further selected the genes with Spearman correlation coefficients over 0.4 and subjected them to pathway enrichment analysis (Fig. 6B). These genes were enriched in pathways related to cell adhesion, ECM remodeling, TGFβ-induced fibroblast migration and ECM production.

To examine the impact of NOTCH signaling on patient prognosis, we analyzed the association between the expression of each NOTCH pathway gene and overall patient survival in 1,434 LUAD samples in the Collaboration dataset. While NOTCH ligands JAG1, JAG2, and DLL3 were associated with worse patient survival outcome, NOTCH3 expression had no correlation (Fig. 6C). Given that NOTCH3 is a mesenchymal cell–specific gene, we hypothesized that its predictive value for survival may be limited to stromal-rich patients. After dichotomizing patient samples into stroma-low and stroma-high subgroups at the median expression level of the consensus stroma/EMT/TGFβ signature, we observed a lack of NOTCH3 association in the stroma-low group (Fig. 6D). However, in the stroma-high group, NOTCH3 and its ligand, JAG1, emerged as worse prognostic markers (Fig. 6E). Given the pivotal role of mesenchymal cells in driving ICB resistance (3, 4, 7, 43), we further examined the impact of NOTCH3 expression on the prognostic value of an 18-gene T cell–inflamed gene expression profile (GEP) that predicts better response to a PD-1 antagonist pembrolizumab across multiple tumor types (43). Again, we binned LUAD samples at the median NOTCH3 expression and found that GEP served as a good prognosis marker only in the NOTCH3-low subgroup (Fig. 6F and G). This suggests that high NOTCH3 expression potentially impedes the antitumor effects mediated by T cells, which are crucial for effective responses to ICB like PD-1 antagonists. To further explore the clinical relevance of NOTCH3 expression and its potential association with ICB resistance, we examined the correlation between NOTCH3 expression and two consensus gene signatures linked to resistance to anti–PD-1 monotherapy in solid pan-tumor types: stroma/EMT/TGFβ and angiogenesis (3). Our analysis revealed a robust positive correlation between NOTCH3 expression and both gene signatures (Fig. 6H; Supplementary Table S7).

In addition, we examined the relationship between oncogenic mutations and NOTCH3 expression using both the TCGA and the Collaboration datasets. Our finding revealed a significant increase in NOTCH3 expression among patients with KRAS or KEAP1 mutation compared with wild-type in both datasets (Supplementary Fig. S7A and S7B). Together, these results underscore the potential clinical relevance of NOTCH3 as a stromal-specific target. Further investigations utilizing preclinical models are warranted to fully comprehend the crucial role of NOTCH3 in tumor stroma remodeling and to evaluate the therapeutic potential of combining NOTCH3 antagonism with ICB monotherapy.

In this study, we have resolved the cellular composition of epithelial, immune, and stromal compartments in freshly resected tumor tissue from a cohort of nine treatment-naïve lung adenocarcinoma patients with a focus on deconvoluting stromal heterogeneity. Using this approach, we have identified the enrichment of FAP⁺PDPN⁺ CAFs and tumor-associated MCAM⁺ACTA2⁺ pericytes in tumors compared with adjacent nontumor tissues. We further map out the spatial architecture of these three major compartments by imaging mass cytometry and demonstrate a close spatial relationship between CAFs, pericytes, and ECs. Ligand–receptor interaction analysis among stromal lineages reveals the increased expression of NOTCH3 in both CAFs and pericytes, whereas several NOTCH ligands are elevated in tumor-associated EC subsets, highlighting a dysregulated NOTCH3 pathway in the tumor stroma (Fig. 6I). The multimodal characterization of the TME in LUAD underscores the importance of the “myCAF”-like FAP⁺PDPN⁺ CAFs, and necessitates a more comprehensive investigation into the critical involvement of NOTCH3 in tumor stroma remodeling, such as genetically engineered mouse models.

NOTCH signaling is a conserved pathway that plays a critical role in developmental cell-fate decisions and has been linked to multiple diseases including cancer (44, 45); however, NOTCH-targeted therapy is not clinically successful (46). One potential reason lies in the inadequate characterization of the expression and functions of each NOTCH receptor in pathologic contexts, as different NOTCH receptors have contradictory biological effects (46). For example, activation of NOTCH1/2 in CD8 T cells triggers a robust and sustained antitumor response, resulting in increased IFNγ production, and reduced tumor burden (47). Thus, nonselective inhibition of the NOTCH pathway will potentially curb antitumor immune responses. Recent scRNA-seq studies underscore the importance of NOTCH3, among other NOTCH receptors, during pathologic tissue remodeling in human liver cirrhosis and rheumatoid arthritis (26, 33). Tuning down NOTCH signaling in human hepatic stellate cells decreases fibrillar collagen production (33), whereas blocking NOTCH3 signaling either by genetic or pharmacologic inhibition attenuates disease severity in an inflammatory arthritis model (26).

In this study, we build on the body of knowledge that has implicated the function of NOTCH3 during pathogenic stromal remodeling and further demonstrate the dysregulated NOTCH3 signaling in the perivascular niche of the tumor stroma in lung adenocarcinoma. Our data indicates the reciprocal interaction of NOTCH3 expressed by pericytes and CAFs and its ligands, such as JAG1 and DLL4, on tumor-associated ECs is likely to shape the tumor-permissive TGFβ-driven “myCAF” phenotype. We further utilized in vitro functional assays to demonstrate that suppressing NOTCH signaling reduces collagen production and matrix invasion of mesenchymal cells. Mesenchymal cells can promote tumor cell invasion in a NOTCH3-dependent manner. Our finding regarding the potential connection between NOTCH3 and COL1A1 expression has also been observed in several in vivo fibrosis models (48, 49). Further investigation is needed to understand the mechanistic role of NOTCH3 in regulating the expression of COL1A1 or other ECM components and whether increased mesenchymal cell invasiveness is directly or indirectly impacted by COL1A1 within the context of lung cancer in vivo.

By investigating patient-derived transcriptomic data, we find that NOTCH3 expression is associated with poor survival in patients with the high expression of the TGFβ-EMT signature. While the T cell–inflamed GEP biomarker does not exhibit prognostic significance in patients with high NOTCH3 expression, it becomes notably significant in patients with low NOTCH3 expression. Given that the T cell-inflamed GEP is a predictive biomarker for response to immunotherapy (43), investigating the potential of NOTCH3 antagonism as a combination strategy to modulate tumor stroma in conjunction with ICB presents an important avenue that warrants further exploration.

H. Xiang reports other support from Merck & Co., Inc., during the conduct of the study and is an employee of Merck Sharp & Dohme LLC, a subsidiary of Merck & Co., Inc. M.A. Sze reports being an employee of Merck Sharp & Dohme LLC, a subsidiary of Merck & Co., Inc. K.A. Van De Mark is an employee of Merck. C.J. Moure reports personal fees from Merck Research Laboratories during the conduct of the study. J. Clarke reports personal fees from Merck Sharp & Dohme LLC during the conduct of the study and is an employee of Merck Sharp & Dohme LLC, a subsidiary of Merck & Co., Inc. T. Tan is an employee of Merck Sharp & Dohme LLC, a subsidiary of Merck & Co., Inc. K.W. Teng is an employee of Merck Sharp & Dohme LLC, a subsidiary of Merck & Co., Inc. M.N. Ramli is an employee of Merck Sharp & Dohme LLC, a subsidiary of Merck & Co., Inc. C. Zhang is a salaried employee of Merck Sharp & Dohme LLC, a subsidiary of Merck & Co., Inc. C. Zhang holds Merck stock. C.E. Gustafson reports grants from Merck during the conduct of the study. S.L. Korle reports other support from Merck & Co., Inc., during the conduct of the study. R. Bueno reports personal fees and other support from Navigation Sciences (equity/patents), Regeneron, Covidien/Medtronic, DiNAQOR; grants from Verastem, Genentech, Roche, Myriad Genetics, Novartis, Siemens, Gritstone, Epizyme, MedGenome, Merck, Bicycle Therapeutics, Bayer, Intuitive Surgical, Northpond, NCI, NIH, DoD, NIBIB, and grants from NHLBI outside the submitted work. L.Y. Moy reports other support from Merck & Co. outside the submitted work. E.H. Vollmann reports other support from Merck & Co. during the conduct of the study. D.Y. Chiang reports other support from Merck & Co during the conduct of the study and other support from Merck & Co. and Bayer outside the submitted work. P.E. Brandish was an employee of Merck & Co., Inc., at the time the project was initiated and for part of the duration on the work. P.E. Brandish owned stock or stock options in Merck & Co., Inc. A. Loboda reports other support from Merck & Co. during the conduct of the study and from Merck & Co. outside the submitted work. No disclosures were reported by the other authors.

H. Xiang: Conceptualization, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft. Y. Pan: Conceptualization, data curation, formal analysis, validation, investigation, visualization, methodology, writing–original draft. M.A. Sze: Data curation, formal analysis, investigation, visualization, methodology, writing–review and editing. M. Wlodarska: Data curation, formal analysis, supervision, writing–review and editing. L. Li: Data curation, formal analysis, methodology, writing–review and editing. K.A. Van De Mark: Data curation. H. Qamar: Data curation, formal analysis, methodology. C.J. Moure: Data curation, writing–review and editing. D.E. Linn: Data curation, project administration. J. Hai: Data curation, formal analysis, writing–review and editing. Y. Huo: Resources, data curation. J. Clarke: Formal analysis, investigation, writing–review and editing. T.G. Tan: Formal analysis, investigation. S. Ho: Data curation. K.W. Teng: Data curation. M.N. Ramli: Data curation. M. Nebozhyn: Formal analysis, methodology. C. Zhang: Formal analysis, methodology. J. Barlow: Project administration. C.E. Gustafson: Project administration. S. Gornisiewicz: Project administration. T.P. Albertson: Project administration. S.L. Korle: Data curation. R. Bueno: Resources, project administration, writing–review and editing. L.Y. Moy: Resources, supervision, project administration. E.H. Vollmann: Resources, supervision, project administration. D.Y. Chiang: Conceptualization, supervision, project administration, writing–review and editing. P.E. Brandish: Conceptualization, resources, supervision, project administration, writing–review and editing. A. Loboda: Resources, formal analysis, supervision, project administration, writing–review and editing.

We would like to thank Dr. Milena Rosa Hornburg for providing gene signatures and input on the manuscript. We also acknowledge Amber Montano and Dr. Alejandro for sharing protocols. This work was supported by the Merck & Co., Inc., Research Laboratories Postdoctoral Research Program.