Immune exclusion and evasion are central barriers to the success of immunotherapies and cell therapies in solid tumors. Here we applied single cell spatial and perturbational transcriptomics alongside clinical, histological, and genomic profiling to elucidate immune exclusion and evasion in high-grade serous tuboovarian cancer (HGSC). Using high-plex spatial transcriptomics we profiled more than 1.3 million cells from 95 tumors and 60 patients, revealing generalizable principles in HGSC tumor t issue organization. Our data demonstrates that effector T cells resist stroma-mediated trapping and sequestration. However, upon infiltration into the tumor, T cells, as well as Natural Killer (NK) cells, preferentially co-localize only with a subset of malignant cells that manifest a distinct transcriptional cell state. The latter consists of dozens of co-regulated genes and is repressed under various copy number alterations. Performing CRISPR Perturb-seq screens in ovarian cancer cells, we identified functionally diverse genetic perturbations 3 including knockout of the insulin sensing repressor PTPN1 and the epigenetic regulator ACTR8 3 that de-repress the proposed immunogenic malignant cell state identified in patients and indeed sensitize ovarian cancer cells to T cell and NK cell cytotoxicity. Taken together, our study uncovered a profound connection between somatic genetic aberrations, malignant cell transcriptional dysregulation, and immune evasion at the cellular and tissue level, allowing us to identify targets that reprogram malignant cell states as an avenue to unleash anti-tumor immune responses.

Multicellular dysregulation plays a key role in the initiation and progression of a wide range of diseases, including cancer, where tumor development and accompanying immune responses depend on (and shape) the location of diverse cell type populations, tissue properties, and organization (137). Cellular and animal models have been instrumental in recovering central immune suppressors (8310) and led to major breakthroughs in cancer patient care. However, many cancer patients do not respond to current immunotherapies (11313), resulting, at least in part, from two central gaps. First , in contrast to the study of cancer genetics, where genome sequencing of tumors across large and diverse patient populations provided a strong foundation to study the genetic basis of cancer and develop targeted therapies, we still lack equivalent maps of tumor tissue organization to study the inherently spatial processes of multicellular dynamics and immune exclusion in patients. Second, identifying the regulators of cell states and reciprocal intercellular interactions poses additional challenges and requires functional interrogation across a larger search space of combinatorial gene-environment perturbations.

In tubo-ovarian high grade serous carcinoma (HGSC) 3 the most common and aggressive form of ovarian cancer ( 14 ) 3 this gap of knowledge is pronounced. HGSC is often diagnosed at adva nced stages, has poor response to current immunotherapies ( 15,16 ), and is prone to chemoresistance, resulting in 5-year survival rate below 50% ( 14 ). Underscoring the need to elucidate the clinically relevant ba rriers that prevent antitumor immunity in HGSC, it is well known that, despite poor response to immunotherapies, abundant tumor infiltrating lymphocytes (TILs) are a robust prognostic marker of better clini cal outcomes in HGSC patients ( 17,18 ). The genetic properties of HGSC have been thoroughly characterized (19 322) 3 demonstrating nearly ubiquitous TP53 mutations, massive copy number alterations (CNA), along with mutation in homologous recombination genes such as BRCA1 and BRCA2 3 and recent single-cell studies provided important resources and insights by characterizing the cellular properties of HGSC in different anatomical sites and genetic backgrounds ( 6,23 ). Yet, the molecular and cellular modalities that promote or suppress immune recruitment and infiltration in HGSC patients remain elusive.

Here, we applied high-plex image-based spatial transcriptomics (ST) wit h subcellular resolution to more than 1.36 million cells across 95 HGSC tumors. Our data demonstrates that effe ctor T cells and NK (T/NK) cells are rarely retained in the tumor stroma. However, T/NK cel l infiltration into the tumor parenchyma is skewed towards subsets of transcriptionally distinct malignant cells, leaving other malignant areas immune deserted. This, together with high-content CRISPR screens, revealed malignant cell transcriptional immunogenicity that is repressed by copy number alterations and can be de-repressed by an array of functionally diverse gene knockouts to sensitize ovarian cancer cells to T and NK cell cytotoxicity. Taken together, our study provides a molecular map of HGSC tumor tissue spatial organization in patients, delineates generalizable principles that predict lymphocyte location and state, and, through integration of spatial and perturbational maps, identified novel targets that reprogram tuboovarian malignant cell states as an avenue to unleash anti-tumor immune responses in this aggressive disease.

To spatially map HGSC in the setting of metastatic disease, we applied in situ imaging with high-plex RNA detection at the single cell resolution to 95 HGSC tumors from a tota l of 60 patients and 136 tissue sections, yielding a total of 1,365,244 high quality single cells9 spatial transcriptomics profiles (Figure 1a, Table S1). Tumor sections were obtained from the adnexa (ovaries/fallopian tube, n = 73), and /or omentum (n = 63), with 37 patient-matched pairs of adnexal and omental tumors. All tum or tissue sections were obtained from debulking surgeries in either the treatment naïve (n = 66) or neoad juvant chemotherapy treatment (n = 70) setting, with associated patient clinica l data including treatment and survival outcomes (Figure 1a, Table S2, Methods). For 40 patients we also obtained DNA sequencing data spanning a 648-gene panel (Figure 1b, Table S1, Methods), focused on actionable single nucleotide variations (SNV), somatic copy number alterations (CNA), chromosomal rearrangeme nts, and tumor mutational burden (TMB), providing a basis to probe the connection between tiss ue structure and somatic genetic aberrations.

The spatial data was collected using three spatial transcriptomic (ST) platforms, allowing rigorous crossplatform validation of these recently developed technologies. A discovery ST dataset spanning 100 tissue sections was profiled with Spatial Molecular Imaging (SMI) ( 24 ), allowing in situ image-based quantification of 960 genes with subcellular resolution (n = 100, formalin-fixed paraffin-e mbedded (FFPE) tissue sections). For comparison and validation, in situ sequencing (ISS via Xenium( 25 ); n = 32, FFPE tissue sections) and MERFISH ( 26 ) (Multiplexed Error-Robust Fluorescence In Situ Hybridization, n = 4, fresh-frozen tissue sections) were applied to profile 280 genes and 140 genes, re spectively (Methods, Table S1b).

Applying a recursive clustering-based cell type annotation procedure (Methods) on processed gene expression profiles (Figure S1a-d, Methods, Supplementary Information) we identified malignant cells (n = 314,191), T and NK cells (T/NK, n = 28,676), B cells (n = 16,373), monocytes (n = 45,549), mast cells (n = 606), fibroblasts/stromal cells (n = 72,861), and endothelial cells (n =

As expected, patients with higher T/NK cell abundance had improved clinical outcomes (p = 5.0*10 -2, Univariate Cox regression Wald Test, p = 3.34*10 -3 log-rank test, Figure S2h), and consistent enrichment of T/NK cells (p = 9.1*10 -3, fisher test) and B cells ( p = 8.29*10 -3, fisher test) is observed in the omental vs. adnexal tumors. Malignant cells and fibroblasts are found to form spatially distinct compartments (i.e., the tumor parenchyma versus the stroma; Figure 1g), with significantly low intermixing between the two cell types (BH FDR < 0.05, hypergeometric test, Figure 1g), such that T/NK cells are preferentially localized in the stroma rather than within the tumor parenchyma (p < 1*10 -4, paired Wilcoxon sum rank test, co-localization quotient, Figure 1h, Figure S2i-k, Methods).

Taken together, these findings demonstrate the quality of our data and validity of our processing pipelines and provide an initial mapping of the tumor organization into spatially segregated compartments (tumor parenchyma and stroma). This dataset sets the stage to probe into cell ular transcriptional states to delineate multi-scale mechanisms underlying immune infiltration and evasion.

Using the ST cohort, we mapped the immune cell intrinsic and extrinsic factors that mark immune infiltration and exclusion. Starting with immune cell intrinsic properties, we mapped immune cell states as a function of their tumor infiltration status, defined based on proximity to malignant cells (Methods). Unsupervised embedding and clustering using single cell expression profiles alone, without any spatial information provided, shows that immune cells residing in the malignant compartment (tumor parenchyma) are transcriptionally distinct from those that reside outside (i.e., i n the fibroblast compartment (stroma), Figure 2a). For each of the five immune cell subtypes robustly represented in the data (CD8 T cells, CD4 T cells, Tregs, NK cells, and monocytes), we ide ntified tumor infiltration associated genes that are significantly (BH FDR < 0.05, mixed-effect, Methods) over or under-expressed as a function of proximity to malignant cells (Figure 2b,c, Figure S3b, Table S4).

The CD8 T cell infiltration program demonstrates that effector and exhausted CD8 T cells are rarely excluded and frequently co-localize with malignant cells (p = 3.24*10 -53, mixed effects). Tumor infiltrating CD8 T cells are characterized (BH FDR < 0.05, mixed-effects) by effe ctor cytotoxicity genes (e.g., GZMB and PRF1) and exhaustion markers ( CTLA4, PD1, TIM3), as well as the pan-cancer exhaustion marker CXCR6 (33335) ( Figure 2b,c) . CD8 T cells that are distant from malignant cells are characterized (BH FDR < 0.05, mixed effects) by naïve and memory T cell markers (IL7R, SELL), overexpress the chemokine receptor CXCR4 (Figure 2b), and reside in the stroma, next to fibroblasts (Figure 2a) . Expanding the CD8 T cell infiltration program to whole-transcriptome level based on scRNA-seq data ( 23 ) ( Methods; Table S4) identified TCF7 3 a central regulator of naïve and resting T cells ( 34,35 ) that directly represses CXCR6 expression 3 as one of the top gene negatively associated with a CD8 T cell infiltration (p < 1*10 -16, rs = 0.23, Spearman correlation).

To investigate the role of the stroma in retaining naïve and memory T cells whilst permitting effector T cells to infiltrate the tumor parenchyma, we integrated sample-matched H&E stains independently annotated by a gynecologic pathologist (Figure S3c), with ST data. Analyzing these data with unsupervised embedding (Figure 2d) and non-linear classifiers ( Methods) revealed two fibroblast subsets, one marking normal adnexal stroma and the other marking desmoplasia (Figure 2d-f, Extended Data Fig 3d-f i.e., a neoplasia-associated alteration in fibroblasts and extracellular matrix with distinct tissue morphology (36340)), which we find to be more prevalent in the omentum ( Figure 2e, Figure S3ej). As expected ( 41,42 ), desmoplastic fibroblasts overexpress collagen fibril organization and extracellular matrix genes (p < 1*10 -2, permutation test, Figure S3d, Table S5), but also upregulate CXCL12 (the cognate ligand of CXCR4, overexpressed in naïve/memory T cells, Figure 2b) and are associated with T/NK rich niches (Figure 2g. p < 1*10 -4, mixed effects).

To systematically map spatially dependent multicellular circuits we identify all the ligand-receptor pairs that show significant (BH FDR < 0.05, mixed effects) spatial co-expression acros s cell types (Methods), revealing suppressive ligand-receptor interactions in the malignant compartment (e.g., CD80/CD86:CTLA4, CD8 T cell:monocyte; TIM3:LAGLS9, CD8 T cell:malignant cel l) and CD8 T cell mediated chemoattraction of other immune cells via CCL2 and CCL5. Co-localization of CXCR6:CXCL16 (CD8 T cell:malignant cell) and CXCR4:CXCL12 (CD8 T cell: fibroblasts) mark chemoattraction cellcell interactions of infiltrating and excluded CD8 T cells, respectively (Figure 2h, BH FDR < 1*10 -10, mixed effects test).

Collectively, these findings demonstrate a differential infiltration process wherein naïve/memory T cells primarily co-localize with the stroma, whereas effector/exhausted T cells reside primarily in the malignant compartment of the tumor (Figure 2a-c). While the data suggests that the stroma is not playing a major role in suppressing or trapping effector T cells ( 43 ) in HGSC patients, we find t hat Tumor Infiltrating Lymphocytes (TILs) are spatially segregated within the tumor parenchyma itself, as described ne xt.

Mapping the spatial distributions of Tumor Infiltrating Lymphocytes (TILs, defined here as both T cells and NK cells) revealed that TILs preferentially co-localize with a subset of malignant cells (Methods, Figure 3a-c, Table S6). Although malignant cell states are highly patient-specific ( Figure S4a) and vary also within patients (Figure S4b-d), the connection between TIL location and malignant cell gene expression appeared repeatedly across the heterogenous tumors in our cohort (Figure 3, Figure S4e-i). Formulating these findings, we identified a Malignant Transcriptional program that robustly marks the presence of Infiltrating Lymphocytes, denoted as MTIL (Figure 3a, Table S6). The program consists of 100 up- and 100 down-regulated genes whose expression in malignant cells is significantly (BH FDR < 0.05, mixed-effects test) positively (M TIL-up) and negatively (M TIL-down) correlated with and predictive of T/NK cell infiltration (Figure 3d,e). M TIL overall expression in malignant cells (Methods) reflects both inter- and intra-sample variation in TIL levels (Figure 3d,e), irrespective of anatomical site ( p < 1*10 -30, mixed effects test; Figure S4g). M TIL continuously increases as a function of T/NK cell abundance and proximity (Figure 3d), also when stratifying the T/NK population into its respective cell subt ypes (Figure seq dataset( 44 ) demonstrates that the M TIL program expression in malignant cells is highest in tumors annotated as <infiltrated=, moderate in tumors annotated as <excluded=, and lowest in tumors annotated as <immune desert= (Figure S4j).

Gene set enrichment analyses demonstrate the connection between MTIL and immune evasion(45350). MTIL-up includes chemokines (e.g., CCL5, CXCL10, CXCL9, and CXCL16 the cognate ligand to CXCR6), and oxidative stress genes (e.g., GPX3, SOD2, Figure 3a,b), and is enriched with multiple immune response genes, including antigen presentation (e.g., B2M, CIITA, HLA-A/B/C), interferon gamma response genes (e.g., IDO1, IFI27, IFIH1, OAS1/2/3, JAK1, STAT1), and cell adhesion molecules (e.g., ICAM1, ITGAV, ITGB2; BH FDR = 1.91*10 -9, 2.86*10-10, 4.59*10 -2, respectively, hypergeometric test, Figure 3b, Table S6b). M TIL-up also includes immune suppression genes, most notable is LGALS9, encoding for galectin 9 3 the ligand of the immune checkpoint TIM3 (i.e., HAVCR2), which is upregulated in the infiltrating T/NK cells (Figure 2h). M TIL-down reflects diverse processes including Wnt signaling (e.g., CTNNB1, FZD3/4/6, SMO, FGFR2, WNT7A), epigenetic regulation ( DNMT3A, HDAC1/11/4/5), as well as genes involved in insulin signaling (e.g., IGFR1, IGFBP5) and cell differentiation (BH FDR < 0.05, hypergeometric test; Figure 3b, Table S6).

The majority of the MTIL genes have never been described in the context of HGSC immune evasion. However, supporting its role in this context, a collection of CRISPR screens(45349) assembled and analyzed here shows that MTIL-up is enriched with genes that sensitize cancer cells to immune mediated selection pressures (including ICAM1, JAK1, NLRC5, SOD2, STAT1, p = 1.82*10 -4, hypergeometric test), while MTIL-down includes genes with desensitizing effects (BCL2, FGFR1, HDAC1, HDAC5, ITGB5, and

Our cohort and independent genomic data suggest that somatic genetic variation intrinsically regulates the MTIL program and in turn impacts T/NK cell influx and non-uniform spread (Figure 4a-c).

First, MTIL overall expression varies across HGSC patients and is associated with improved overall 1*10-30, ANOVA test).

Second, aligned with the finding that malignant cell transcriptomes are tightly associated with CNAs in cis (Figure S3a), M TIL expression strongly correlates with the copy number of multiple genes in our cohort, the top ones being IFNGR2 and IFNAR1 (positively correlated) and TCF7L2, FGFR2, and AXL (negatively correlated, p < 5*10 -3, mixed effects, Methods, Figure 4b-c).

Third, CNAs of MTIL genes are predictive of TIL abundance scores (Methods) in an independent TCGA cohort of 578 HGSC tumors( 20 ) (AUROC = 0.82, on unseen test samples, supervised vector machi nes (SVM) model, Methods), where tumors with amplification of M TIL-down genes (e.g., DNMT3A, FZD3, MYL9, SRC, and TGFB2) or deletion of M TIL-up genes (e.g., CX3CL1, CXCL10, CXCL9, ICAM1, GPX3, NR3C1) have significantly lower TIL abundance scores compared to tumors without these copy number changes (BH FDR < 5*10 -3, one-sided t-test, Figure 4e).

These findings propose a genetic basis to immune evasion and tolerance in HGSC, where the transcriptional CNA-driven malignant cell states can impact cancer-TIL interactions and shape TIL recruitment. To examine this model, we turned to identify regulators controlling the MTIL program and examine their functional impact on cancer cell response and susceptibility to T/NK cell-mediated cytotoxicity. cytotoxicity

To functionally probe the MTIL program genes and examine their effect on cancer cell response to lymphocyte cytotoxicity, we performed high content CRISPR knockout (KO) screens in ovarian ca ncer cells in monoculture and co-culture with cytotoxic lymphocytes, including T Cell Receptor (TCR)engineered CD8+ T cells and NK cells. Using this approach, we sought to functionally identify and distinguish between co-regulated immune response and immune suppressive genes captured by the MTIL program (e.g., ICAM1 and LAGSL9) and identify perturbations that trigger the former.

Instead of targeting only genes in the MTIL program itself, we devised a meta-analysis pipeline to identify program regulators based on available Perturb-seq datasets (Methods). Using four previously published Perturb-seq datasets(51353), we identified 43 and 104 perturbations that result i n significantly higher or lower expression of the program, respectively (Figure 5, Table S6, Methods). Demonstrating the value of this approach, it revealed a wider and more diverse set of regulators, most of which are not included in the MTIL program itself or not included in the spatial data gene panels (Table S1b). Negative M TIL regulators are enriched for chromatin organization (e.g., DNMT1, INO80, TAF10, WDR5), Wnt pathway, Myc targets, and immune resistance genes ( 45349,54 ) (BH FDR < 1*10 -3, hypergeometric test). The top negative regulator identified here is PTPN1, which is supported by both gene activation and inhibition (Figure 5a,c) experiments.

This approach guided our design of the pooled knockout of 74 MTIL genes and regulators (Table S7) in ovarian cancer cells (TYK-nu cell line, Figure 6a, Figure S5-6). Mapping fitness upon genetic perturbations under both innate and adaptive immune selection pressures (Figure 6a,b, BH FDR < 0.05, MAGeCK, Methods) along with Perturb-seq scRNA-seq readouts in monoculture and co-culture with NK cells (Figure 6a,c), allowed us to identify perturbations that activate or repress the program a nd track subsequent effects of these perturbations on immune escape. In total we profiled 18,585 high quality single cell transcriptomes, each assigned to an ovarian cancer cell with a single sgRNA confidently identified, and a median of 4,251 genes detected per cell (Figure 6c, Figure S7a). Differentially expressed genes were identified for each gene knockout across the three conditions (fisher method; Methods), resulting in 74 gene <perturbation signatures= (Methods) that were then used to identify gene knockouts that significantly repress or activated the MTIL program, denoted as -activators and -repressors, respectively (Figure 6d, Methods).

Validating our hypothesis and approach, the top perturbations activating the program 3 PTPN1 and ACTR8 knockouts 3 sensitize malignant cells to T/NK cell cytotoxicity (Figure 6b,d-e, Figure S7b), while the top perturbations that repress the program, IFNGR1 and STAT1 knockouts, allow ovarian cancer cells to resist T cell mediated killing (Figure 6b,d-e, Figure S7b). Knockout of ACTR8 and PTPN1, as well as other top repressors FGFR1, MAPK1, and MED12 were found to sensitize cancer cells to immune elimination also in based on data from previous in vivo CRISPR screens(45349). Moreover, we find that knockout of MTIL repressors (ACTR8, DNMT1, FGFR1, PTPN1, MED12, and MIF) mimics and amplifies the transcriptional responses to NK cells, while knockout of MTIL activators, as STAT1, IFNGR1, INTS2, IRF1, PARP12 and others, represses and counteracts the transcriptional response to NK cells (Figure 6fh, Extended Data 7c-f). Lastly, knockout of specific genes within the program, including GPX3 and

TAGLN show substantial impact on the ovarian cancer cell susceptibility to NK mediated killing (Figure 6b), demonstrating both global and gene-specific effects.

Taken together, coupling HGSC spatial tumor organization with multimodal functional probing, identified new and clinically relevant targets to sensitize ovarian cancer cells to innate and adaptive cytotoxic lymphocytes and demonstrated the role of cancer cell intrinsic transcriptional dysregulation as an important driver dictating the outcomes of the malignant-T/NK cell interplay.

Our study maps the tumor tissue landscape in HGSC patients and reveals generalizable principles of tissue organization that dictate lymphocyte location and state within these aggressive and genetically unstable tumors. It uncovers a profound connection between somatic genetic aberrations, malignant transcriptional dysregulation, and immune evasion at the cellular and tissue level, providing a new perspective to the barriers preventing the anti-tumor immune response in HGSC patients and new leads to derepress HGSC cancer immunogenicity.

Innate and adaptive cytotoxic lymphocytes (CTLs) have a substantial effect on c ancer cell transcriptome (Figure 6c). As shown here, genetic dysregulation that prevents this transcriptional respons e can have significant effects on cancer cell susceptibility to immune elimination even in the highly controlled cocultures as those used here, where CTLs are already primed and activated, and spatial segregation is unlikely to occur. These effects can be amplified in the context of in vivo cancer-immune co-evolution where immune tolerance is reinforced due to positive feedback loops across cells. Indeed, immune checkpoint blockade and other immunotherapies have shown modest effects in tumors with low TIL levels at baseline ( 17,55 ). The data shown here proposes that this may not be only due to i mmune exclusion per se, but also due to cancer intrinsic differences between TIL-rich and TIL-deprive tumors that protect malignant cells even in the presence of targeting CTLs. Our findings and approach open new directions for further investigation of the genetic basis of tumor immune evasion through the lens of spatial organization and put forward a framework to design targeted strategies to counteract or bypass these resistant mechanisms.

More generally, as more spatial datasets become available, there is a growing need to use this rich information to delineate new drivers of complex multicellular processes and phenotypes. Here we show the value of mapping spatial cell states to genetic information across individuals and to design perturbational screens with single cell readouts. Importantly, we show that using existing Perturb-seq datasets to identify latent regulators of gene expression programs is critical and provides a data-driven framework to uncover regulators that are not necessarily included in the program itself. As more Perturbseq datasets, as the one generated here, become available across a more diverse range of cell types and conditions, it will be possible to use this information more effectively to extrapolate from one context to another with increasing accuracy ( 56,57 ).

The key findings from our study can fuel new lines of investigation towards new clinical interventions in HGSC patients. We anticipate that the detailed mapping of HGSC tumors provided here will help inform the design of new T/NK cell engineering strategies to reach better control of cell delivery and location in a more precise manner that is aligned with the tumor cellular and molecular structure in patients. Our findings demonstrate that the stroma forms a differential <filter= that supports differential occupancy of effector T cells in the malignant compartment 3 this calls for dynamic tracking of tumor reactive T/NK cells across non-tumor sites (i.e., in the circulation and lymph nodes) and withi n the tumor to help elucidate this process and examine if T cells can also egress back to the stroma cells to avoid or reverse exhaustion( 58 ) and how to best leverage, as opposed to eliminate or target, the stroma. Our data provides new leads to target HGSC resistance, including epigenetic regulators (e.g., ACTR8 and MED12), fibroblast growth factor receptors ( FGFR1/2), GPX3, and PTPN1. PTPN1, which we found to be one of the most potent MTIL repressors, provides pre-clinical rationale to test new PTPN1/N2 inhibitors (NCT04777994, NCT04417465, phase I clinical trials)(59364) in HGSC patients, and demonstra tes a connection between immune evasion, insulin resistance, and type 2 diabetes. PTPN1 is a negative regulator of insulin and leptin signaling ( 65 ) that has been an attractive drug targe t for treatment of type 2 diabetes and obesity (66369). PTPN19s protein product PTP1b is inactivated by oxidation( 70 ), which may explain MTIL activation under oxidative stress (as indicated by the up regulation of GPX3 and SOD2). Further supporting the connection to insulin resistance, TCF7L2, which we identified as a top gene amplification associated with the repression of the MTIL program in the HGSC cohort (Figure 5b-c) harbors the most significant SNP associated with type 2 diabetes risk ( 71 ).

Taken together, this integrative study provides a blueprint to functionally map and probe the molecular landscape of multicellular interplay in complex biological tissues and reveals unrecognized spatial, molecular, and genetic regulation of immune escape in HGSC, opening new avenues to activate targeted immune responses in this aggressive disease. collected across three platforms: SMI (discovery dataset), ISS (validation da taset 1), and MERFISH (validation dataset 2); n denotes the number of tissue sections profiled. (b) Clinical annotations of the patients and samples included in the cohort. (c) Uniform Manifold Approximation and Projection (UMAP) embedding of cell transcriptomes from the discovery dataset (top left), validati on dataset 1 (top right), and validation dataset 2 (bottom left). Cells are colored according to their cel l type annotations. n denotes number of cells with each cell type annotation. (d) Representative ST images (right) and corresponding H&E (left, where available) depicting cell segmentations with each ce ll colored based its cell type annotations. (e) Co-embedding spatial cell transcriptomes from this study with publicly available scRNAseq datasets ( 27,29332,72,73 ). Unified UMAP of co-embedded cell transcriptomes is s hown with cells colored by cell types (top) and dataset (bottom). ( f) Cell type composition (y axis) per sample (x axis) from this study and in publicly available scRNA-seq HGSC cohorts( 27,29332,72,73 ). ( g) Pairwise colocalization analysis: the number of samples (x axis) where each pair of ce ll types (y axis) 133 shows significantly (BH FDR < 0.05, hypergeometric test) higher (red), lower (blue), or expec ted (grey) colocalization frequencies compared to those expected by random. (h) Log2 Co-localization Quotient (CLQ, y-axis) of T/NK cells with fibroblasts (blue, x axis) and T/NK cel ls with malignant cells (green, x axis) in each tissue section from the discovery dataset ( **** p < 1*10 -4, paired Wilcoxon rank sum test).

Light grey lines connect paired fibroblasts and malignant cells within each tissue section. Boxplots middle line: median; box edges: 25th and 75th percentiles; whiskers: most extreme points that do not exceed ± IQR x 1.5; further outliers are marked individually. considering all genes (top) or only T cell specific genes (bottom, for further confirm ation). Cells are colored according to the frequency of malignant cells ( 1, 2 ) or fibroblasts ( 3 ) in the T cell microenvironment (Methods), the overall expression (OE) of the CD8 T cell infiltration program ( 4, 5 ) or their k-Nearest Neighbor cluster ( 6 ). ( b) CD8 T cell tumor infiltration program, showing the association (p-value and effect size) of each gene (row) with infiltration status, when conside ring only specific immune cell subsets (columns). ( c) Representative ST images from validation dataset 1 depicting the CD8 T cell tumor infiltration program identified in the discovery dataset. Malignant cells are in grey, CD8 T cells are colored according to the Overall Expression (OE) of the infiltrat ion program identified in the discovery dataset (color bar). The respective p-values denote per tissue section if the OE of the CD8 T cell infiltration program is significantly higher in CD8 T cells with a high (above median) vs. low (below median) abundance of malignant cells within a radius of 30¿m (one-sided t-test). ( d) UMAP embedding of fibroblast cell ST profiles colored by stromal morphology (left) and anatomical site (right) annotations. (e) Average gene expression (z-score, red/blue top middle color bar) of the top 50 desm oplasia associated genes (columns) across fibroblasts in each sample (rows), sorted by overall expre ssion score of the 50 genes (Methods, left color bar), and labeled by their anatomical site (middle color bar) and st romal morphology annotation (right color bar). ( f) Representative tissue section (HGSC24, adnexa, discovery dataset) wherein the desmoplasia associated genes capture intra-tumoral differential stromal morphology (p-value = 7.23*10 -80, Wilcoxon rank sum test). Hematoxylin & Eosin Stain (left), Cell Types in situ (middle), and cells plotted in situ with fibroblasts colored according to the overall expression (OE) of desmoplasia associated genes and all other cells in grey. (g) OE score of top 50 desmoplasia associated genes per fibroblast (x axis) as a function of T/NK cell density within a 30-¿m radius (y axis) in the adnexa (left) and in the omentum (right). ( h) Ligand-receptor interactions (lines) consisting of genes from the CD8 T cell infiltration program (up-regulated in yellow, down-regulated in dark purple) and their respective ligand/receptor in the cancer compartment (light blue, i.e., tumor infiltrating programs of other immune cells or genes specific to cancer cells in T cell rich areas) and stroma compartment (light purple, i.e., genes specific to fibroblasts in T cell rich areas). The arrows conne ct each gene to the cell type where it was found to mark the respective spatial pattern, namely, tumor infiltration in immune cells, and colocalization with T/NK cells in the non-immune (fibroblast or malignant) cells. cells. (a) Heatmap of MTIL genes. Average expression (z score, red/blue color bar) of the M TIL genes (rows) across spatial frames (columns), sorted by M TIL overall expression (OE), and labeled by (color bar from top to bottom): anatomical site, sample ID, the detection of different T /NK cell subsets. (b) MTIL gene ontology enrichment analysis. (c) Spatial distribution of T/NK cells (black) and M TIL OE in malignant cells (color bar, top right) shown in representative tumor tissue se ctions from six different patients and anatomical sites; other non-malignant cell types are colored grey. p-values denote if MTIL OE is significantly (one-sided t-test) higher in frames with high vs. low T/NK a bundance (defined based on the median level) in the respective tissue section. ( d) MTIL OE (y axis) in malignant cells, stratified based on the relative abundance of T/NK cells in their surroundings (top) and the prese nce of T/NK cells at different distances (bottom). Middle line: median; box edges: 25 th and 75th percentiles; whiskers: most extreme points that do not exceed ± IQR x 1.5; further outliers are marked individually. ****p < 1*10 -4, mixed effects (Methods). ( e) ROC curves obtained for cross-validated Support Vector Machine classifier using MTIL expression in malignant cells to predict T/NK cell levels, at the sample (black), spatial frames (red) and single cell levels (blue). Abbreviations: AUROC = area under the ROC c urve. (f) Single cells visualized in situ in one representative whole tissue section from validation dataset 2 (left), juxtaposed with magnified region (right). T/NK cells are in black, malignant cells a re colored via normalized MTIL overall expression in the color bar, and non-malignant cells are in grey (MTIL expression in TIL-high versus TIL-low niches, p = 2.2*10 -16, Wilcoxon rank sum test). (OE) in malignant cells (y axis) residing in spatial frames where T/NK cells were not detected, stratified by patients (x axis). ( b) MTIL OE in malignant cells (y axis), stratified by somatic copy number of the respective gene (x axis) based on patient-matched bulk tumor genomic profiles. ( c) Top CNAs showing a significant (BH FDR < 0.05, mixed-effects; Methods) positive (red) or negative (light blue) association with MTIL OE in malignant cells in the discovery spatial cohort. (d) Kaplan Meier Survival curves depicting differential survival probability (y axis) as a function of average M TIL OE in the malignant cells of each patient (log rank test p = 4.09*10 -2) ( e) Deletion (red) of M TIL-up genes and amplification (light blue) of M TIL-down genes (x axis) are signifcantly (BH FDR < 0.05, one side t-test) associated with low T/NK levels (y axis; inferred based on gene expression of T/NK cell signatures) in an independent TCGA

HGSC cohort of 578 patients ( 20 ). Grey distribution depicts the T/NK levels in tum ors without the respective genomic abberation.

A subset of replicated cells was allowed to recover for an additional day until confluency prior to being snap frozen and stored at -80#. Genomic DNA of the snap frozen cells was extracted using the QuickDNA Midiprep Plus Kit (Zymo Research, D4075). sgRNA amplification was perform ed following a previously published protocol( 49 ). Equimolar amounts of indexed libraries were pooled togethe r and sequenced on a MiSeq Nano V2 in a single-end run at the Stanford Genomics Service Center. A second screen was performed for sgRNA sequencing. TYK-nu Cas9 library cells we re seeded in a 6 well dish (Cole-Parmer, 0192770) and were allowed to adhere overnight. NK-92 cells were added at a 2.5-to-1, 5-to-1, and 7.5-to-1 effector-to-target ratios for 48 hours. TYK-nu cells were allowed to recover for three days prior to being snap frozen and prepared for genomic DNA extraction as described above. Each experimental condition was performed in triplicates with > 1000x cells per sgRNA, resulting in 6 and 12 sequencing samples from the first and second screen.

Raw fastq files were processed using the cellranger pipeline (10x Genomics Cell Ranger 7.1.0). Counts were converted to transcript per million (TPM) values. For each condition (monocul ture, 1:1 co-culture, and 2.5:1 co-culture) data was analyzed to remove non-malignant cells. Seurat R package was used for KNN clustering, resulting in a distinct NK cluster in the co-culture conditions, with expression of CD3E and NCAM1. This cluster was removed and only cancer cells with a detection of a single sgRNA were retained for downstream analyses. For each of the three conditions, DEGs were identified for each perturbation using a two-sided t-test comparing the cells with the perturbation to those with NTCs. Fisher test was used to combine the three p-values. Hypergeometric tests were performed to examine if the up or down regulated genes identified for each perturbation were enriched with MTIL-up or MTIL-down genes, or vice versa, and the combined p-value (fisher test) was reported as the final summary statist ics. MAGeCK algorithm was used to compute differential fitness effects on the cancer cells under the monoculture and co-culture conditions, either with the different types of CD8+ T cells or with the NK cells. Each experimental condition was performed in triplicates. First, the sgRNA counts of the different samples were median normalized to adjust for the effect of library sizes and read count distributions. Second, the variance of read counts was estimated by sharing information across the different sgRNAs, allowing to fit a negative binomial (NB) model to test whether sgRNA abunda nce differs significantly between treatments (i.e., co-culture) and controls (i.e., monoculture or co-culture wit h non-specific T cells). Third, sgRNAs were ranked based on p-values calculated from the NB mode l, and an ³ robust ranking aggregation (³-RRA) algorithm was used to identify positively or negatively sel ected genes. The pairwise tests were performed considering each of the co-cultures (CD8 T cell or NK cell) compared to the monoculture and the two co-cultures (with specific and non-specific T cells) to compute a combined Fisher statistics, one for T cell and another for NK cell sensitizing and desensitizing hits.

All the data collected in this study, including spatial transcriptomics data, single-cell Perturbseq data, targeted genomics, deidentified clinical meta-data, and processed tissue images will be deposited and made publicly accessible through Gene Expression Omnibus (GEO), Zenodo, and CEL LxGENE. Processed data in the form of standardized RObjects will be available via Zenodo. Upon publication raw and processed spatial transcriptomics data will be available on CELLxGENE for download in .h5ad format and interactive exploration. Processed gene expression matrices with cell type annotations from 6 scRNAseq studies with HGSC tumor samples were downloaded from publicly available repositories specified in their respective publications ( 23,27332 ). Specifically, preprocessed gene express ion and metadata matrixes of HGSC scRNA-seq data were downloaded from Synapse (syn33521743 ( 23 )), GEO (GSE118828 ( 27 ), GSE173682 ( 29 ), GSE147082 ( 28 ), GSE154600 ( 31 ), GSE146026 ( 80 )), and https://lambrechtslab.sites.vib.be/en/data-access ( 30,32 ). An additional external validation dataset ( 44 ) hosted on the European Genome-Phenome Archive (EGAD00001006973, EGAD00001006974) was made available for this study through a Data Access Agreement with Genentech, Inc..

All data processing and analysis code will be available via GitHub once the paper is published. The study GitHub repository includes documented code for data processing and code required to reproduce the main figures and supplementary tables of the study.

Weeden CE, Gayevskiy V, Marceaux C, Batey D, Tan T, Yokote K, et al. Early immune pressure initiated by tissue-resident memory T cells sculpts tumor evolution in non-small cell lung cancer. Cancer Cell. 2023 May 8;41( 5 ):837-852.e6.

We thank Reece Villarin Akana, Georgia Schmitt, Grace Allard, and Kaitlyn Spees for help with setting up assays and models related to this work; Lettie McGuire for help with artwork; Kai Wucherpfennig and Nathan Mathewson for generously providing NY-ESO-1 constructs; Theodore Roth for discussions on primary T cell culture; Stanford Genomics Service Center (SFGF), Norma Ne ff and the rest of the Stanford Chan Zuckerberg Biohub Center team for ongoing support. L.J. is a Chan Zuckerberg Biohub Investigator and an Allen Distinguished Investigator. L.J. holds a Career Award at the Scientific Interface from the Burroughs Wellcome Fund and a Liz Tilberis Early Career Award from the Ovarian Cancer Research Alliance (OCRA). This study was supported by the OCRA The Liz Tilberis Earl y Career Award Grant 889076. This work was support by Under One Umbrella, Stanford Women9s Cancer Center, Sta nford Cancer Institute, an NCI-designated Comprehensive Cancer Center, and supported in part by funds from the Departments of Pathology and Genetics at Stanford University. Processed spatial and scRNA-seq Perturb-seq data will be available on the Gene Expression Omnibus (GEO) and Zenodo. Figure 1a, 7a and (Licensing

C.Y.Y., K.A., O.L., B.E.H, and L.J., designed the study. C.Y.Y., K.A., O.L., B.E.H, and L.J. conducted the research and interpreted the results. K.A., O.L., and S.K., preformed the experiments with L.J. supervision. A.W., B.L., X.Z., L.H., and B.E.H. collected, annotated, and provided tissue samples for 43 profiling. C.Y.Y., and L.J. performed the computational and statistical analyses. R.V. supported experimental design. M.C.B. supported experimental design and manuscript writing. M.P.S. supported All authors reviewed and approved the manuscript.

M.P.S. is a co-founder and scientific advisor of Personalis, SensOmics, Qbio, January AI, Fodsel, Filtricine, Protos, RTHM, Iollo, Marble Therapeutics and Mirvie. He is a scientific advisor of Yuvan, Jupiter, Neuvivo, Swaza and Mitrix. M.C.B has outside interest in DEM Biopharma. The remaining authors declare no competing interests. Figure S1. Cell segmentation and cell type annotations of spatial transcriptomics. (a) Representative 45 whole-cell segmentation performed for the discovery dataset. Input data includes

Table S3. Cell Type Signature Genes. Includes both gene signature derived from scRNA-seq and CellTypist Immune Encyclopedia (a) and from the HGSC spatial transcriptomics data coll ected here (b). Table S4. Immune tumor infiltration signatures derived from discovery dataset, shown for different immune cell subsets (a) and for the CD8 T cell infiltration program expanded to whole-transcriptome using scRNA-seq data of CD8 T cells (b).

Table S5. Fibroblast desmoplasia associated genes (a) and their Gene Ontology enrichment analysis

Table S6. MTIL genes (a) and its Gene Ontology enrichment analysis (b).