Spatial Heterogeneity of Infiltrating T cells in High-Grade Serous Ovarian Cancer revealed by Multi-omics analysis

doi:10.21203/rs.3.rs-1053777/v1

Download PDF

Article

Spatial Heterogeneity of Infiltrating T cells in High-Grade Serous Ovarian Cancer revealed by Multi-omics analysis

https://doi.org/10.21203/rs.3.rs-1053777/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Tumor-infiltrating lymphocytes (TILs), especially CD8⁺ TILs, represent a favorable prognostic factor in high-grade serous ovarian cancer (HGSOC) and other tumor lineages. The spatial heterogeneity of different TIL subtypes in HGSOC remains to be elucidated. We integrated RNA sequencing, whole-genome sequencing, bulk T cell receptor (TCR) sequencing, as well as single-cell RNA/TCR sequencing to investigate the characteristics and differential composition of TILs across different HGSOC sites. Two immune patterns in ovarian cancer are identified: 1) ovarian lesions with low infiltration of mainly dysfunctional T cells and immunosuppressive Treg cells, 2) omental lesions infiltrated with non-tumor-specific bystander cells. Exhausted CD8 T cells that are preferentially enriched in ovarian tumors are clonal expanded and exhibit evidence for cytotoxic activity. Inherent tumor immune microenvironment characteristics appear to be the main contributor to the spatial differences in TIL status. The landscape of spatial heterogeneity of TILs informs potential strategies for therapeutic manipulation in HGSOC.

TILs

spatial heterogeneity

multiple lesions

multi-omics

High-grade serous ovarian cancer (HGSOC) affects 239,000 women worldwide each year and represents the most lethal type of gynecological cancer¹. Almost 80% of patients are diagnosed as stage III or IV disease, and many succumb to primary treatment-resistance or relapse within 18 months, leading to a 5-year survival rate of about 30%². Unfortunately, the survival odds for HGSOC patients have not improved markedly despite years of extensive biological research and clinical trials and the addition of bevacizumab and PARP inhibitors to the therapeutic armamentarium. There is thus a strong need for effective new therapies including ones that induce effective immune engagement potentially through immune checkpoint blockade (ICB), such as programmed cell death (PD-1) or its ligand (PD-L1) antibodies³with median response rates lower than 15%⁴. The mechanism(s) underlying the lack of response to ICB despite the presence and prognostic impact of T-cell infiltration remains largely unknown.

HGSOC often presents with widespread abdominal cavity dissemination with omentum as the most frequent site of metastasis ^{5, 6}. Multi-site studies, albeit controversial⁷ indicate that genomic inter-lesional heterogeneity⁸ is associated with poor survival⁹. The effect of spatial immunologic variation, especially in T cell infiltration, recognition and expansion, across various tumor foci in the ovary (primary) and distant metastatic foci in the peritoneal cavity and their contribution to the limited response to immune therapy in HGSOC remains unexplored.

To provide a detailed analysis of the landscape of heterogeneity of infiltrating T cells in primary and metastatic lesions and their differential characteristics in HGSOC, we performed multi-site sampling and simultaneous RNA sequencing (RNA-seq), whole-genome sequencing (WGS) and bulk T cell receptor (TCR) sequencing as well as single cell RNA sequencing (scRNA-seq) and paired TCR sequencing (scTCR-seq) in 9 patients (48 sites) with untreated primary HGSOC. We identify two different immune patterns in ovarian cancer: 1) ovarian lesions with low infiltration of mainly dysfunctional T cells and immunosuppressive Treg cells. These exhausted CD8 T cells with cytotoxic function are clonal expanded; 2) omental lesions infiltrated with non-tumor-specific bystander cells. Decreased major histocompatibility complex (MHC)-I antigen presentation ability and failure of T cell infiltration into omental tumors may contribute to lack of tumor-specific T cells in omental metastasis. Together these observations may partly explain the poor response of ovarian cancer to current immunotherapy approaches.

Differential transcriptomic profiles across multiple sites in HGSOC

We performed WGS and RNA-seq on 48 sites from nine treatment-naïve pathological HGSOC patients (Supplementary Table 1, Fig. 1a). Each site in the same patient had a similar proportion of tumor (Supplementary Table 1, Methods). Copy-number variation (CNV) and somatic mutations of these tumors were consistent with known HGSOC genomic patterns ¹⁰ (Fig S1a-b) with the exception of OV001, that despite high grade characteristics on histopathology, did not have TP53 mutations and instead with the only tumor with a NF1 aberration. High level gene amplifications present in ovarian cancer, such as CCNE1, MYC were present in 2/9 (Fig S1b). Eight of the nine HGSOC tumors had somatic TP53 mutations, while 2/9 patients harbored germline BRCA1/2 mutations (Fig S1a-b). Importantly, the majority of the CNV and somatic mutation events did not demonstrate spatial genomic heterogeneity among tumor sites (Fig S1b). The detection of TP53 mutations in metastases but not in ovarian sites of OV004 is one key exception. Manual inspection of sequence tracks failed to identify TP53 mutations in ovarian sites. To characterize the relationship between multiple sites in HGSOC, we first performed principal component analysis (PCA) on transcriptomic profiles of primary ovarian (Ov while HGSOC originate in the fallopian tube, the ovary represents the most frequent site of initial seeding consistent with the definition of primary), omental (Om) and other metastatic lesions (Ot). PCA demonstrated two drivers of heterogeneity, patient specific processes with tumors across different sites within a given individual tending to cluster together and tumors within different sites tending to cluster together for different patients (Fig. 1b). For example, while ovarian tumors from OV004, OV005, OV006 and OV008 were clearly separated from omental and other sites, ovarian tumors from OV001, OV002, OV003 and OV009 tended to cluster closer to their metastatic sites than to other ovarian tumors. Decomposition of immune cell proportions using ssGSEA analysis of RNA-seq data recapitulated the heterogeneity observed in the PCA analysis with information content specific to patients and also to tumor site (Fig. 1c). Figure 1c also demonstrates the robustness and the consistency of the analysis with tumors from the left and right ovary from the same patient clustering together and multiple different omental lesions from the same patient clustering together. Similar to the PCA analysis, immune cell-based clustering suggested that while most ovarian tumors were in a single cluster, OV001, OV002, OV003 and OV009 ovarian tumors tended to cluster with their metastatic sites. Interestingly, the estimated proportions of various immune components were low in ovarian tumors (Fig. 1c). In contrast the immune components were markedly higher in most omental sites compared to matched ovarian tumors; with other sites have lower immune content and indeed a subset of the other metastatic sites clustered with the ovarian tumors (Fig. 1c). A panel of 159 genes selected based on six different characteristics of T cell quantity or spatial distribution¹¹ were used to further characterize the samples demonstrating highest T cell infiltration in omental lesions, with most of the ovarian tumors having low levels, indicating a “desert” T cell phenotype (Fig. 1d). In parallel, CD4⁺ and CD8⁺ T cell infiltration into the different lesions were evaluated based on immunohistochemistry (IHC) staining (Fig S1c) and flow cytometry analysis (Fig S1d). Consistent with the transcriptional profiling data, the density of CD4⁺ and CD8⁺ T cells was much lower in ovarian lesions than in omental lesions (Fig. 1e and f), indicating that the ovarian lesions are immune ‘cold’ lesions. FAP, a marker of activated stroma, in contrast, did not vary across lesion location (Fig. 1e).

Distinct characteristics and differential composition of TILs across different lesions in HGSOC by scRNA-seq

To further detail the landscape of infiltrated T cells and explore the heterogeneity among different lesions, we sorted CD45⁺CD3⁺ T cells from single-cell suspensions prepared from 13 ovarian (Ov), 7 omental (Om), 4 other distant metastatic (Ot) sites and 5 PBMCs of patients OV004, OV005, OV006, OV008, OV009, and OV010 and performed scRNA-seq and matched scTCR-seq using the 10x 5’ platform (Fig. 1a and S2a, Supplementary table 1). After removing confounding batch effects and patient-specific variability (see Methods), a total of 227,769 CD45⁺ CD3⁺ immune cells from all subjects were available for analysis.

Using unsupervised clustering of uniform manifold approximation and projection (UMAP), we identified 22 stable clusters, including 7 clusters for CD4⁺ and 15 clusters for CD8⁺ T cells, each with unique signature genes (Figs. 2a-b and S2b-e). In addition to typical CD8⁺ and CD4⁺ T cell clusters including naïve (Tn), effector (Teff), memory (Memory), mucosal associated invariant T cells (MAIT) of blood and tissue, conventional regulatory T (Treg), and dysfunctional “exhausted” T cells (Tex), we also identified two proliferative clusters that highly expressed MKI67: CD8_C05-TYMS expressing markers associated with exhaustion (designated as Tex.prol) and CD4_C04-TYMS expressing markers associated with Treg (designated as Treg.prol) (Fig. 2b and S2d-e). CD8_C03 (Tex) population showed the highest expression of CXCL13, HAVCR2, and the co-inhibitory receptor PDCD1, as well as increased expression of GZMB, GZMA and GZMH, indicating that cells in this cluster potentially have cytotoxic activity in addition to exhaustion features. Furthermore, a pre-dysfunctional cluster (CD8_C02, referred to as “transitional”) was defined by high expression of GZMK¹² and a progenitor exhaustion cluster (CD8_C07) was defined by higher GPR183¹³ (a central memory marker) and lower PDCD1 (an exhaustion marker) than CD8_C03 (Tex). We identified additional CD8 positive subsets, including CD8_C04 (NK-like) and CD8_C15 (γδ-like). CD8_C04 (NK-like) highly expressed KLRD1 and NKG7, known markers of NK¹⁴/NKT¹⁵ cells and CD8_C15 (γδ-like) highly expressed TRDV2 and TRGV9, known markers of γδ T cells¹⁶ (Fig. 2b and S2e).

We next investigated the relative proportions of different clusters between ovarian, omental and other sites and blood (Fig. 2c-d and S2f-i). Interestingly, the proportion of dysfunctional cells, including CD8_C03 (Tex), CD8_C05 (Tex.prol), and CD8_C07 (Tex.prog) and immunosuppression cluster, CD4_C02 (Treg) were significantly enriched in ovarian tumors (Fig. 2d). CD4_C03 (Tex) exhibited a trend to increase in ovarian as compared to omental or other sites. In contrast, naïve, memory and transition functional state clusters were enriched in omental sites (Fig. 2d). Opal multiplex IHC showed that most T cells enriched in ovarian lesions were dysfunctional states (CD8⁺ PD-1⁺), whereas memory T cells (CD8⁺ PD-1⁻ or CD8⁺GZMB⁻) enriched in omental lesions (Fig. 2e). The major cellular composition difference between ovarian and omental lesions was again observed by flow cytometry analysis (Fig. 2f-g and S3a). Treg and Tex cells were significantly increased in ovarian sites, while central memory T (Tcm) were enriched in omental lesions, and other subsets, including naïve, effector memory T and effector memory re-expressing CD45RA T cells (T_EMRA) were comparable in these two sites (Fig S3b). Taken together, increased Tex and Treg is consistent with primary ovarian tumors being immunosuppressed.

Patient derived TMB is associated with skewed T cell differentiation

Tumor mutation burden (TMB), neoantigen burden and high genomic instability, including deficient mismatch repair (dMMR) and homologous recombination deficiency (HRD), have been associated with increased T-cell infiltration and better response to checkpoint inhibitors in some cancer types^{17, 18, 19}. To explore whether heterogeneity in T cell infiltration in different tumor sites or different patients is related to genomic aberrations, TMB²⁰, HRD²¹ score and COSMIC mutational signature²² of each sample were assessed according to previous analysis pipelines (Fig S4a-b). Concordant with a previous study in NSCLC²³, the correlation matrix revealed that CD4_C03 (Tex) and CD8_C03 (Tex) clusters correlated with TMB, neoantigen burden and HRD score, suggesting that CD4_C03 (Tex) and CD8_C03 (Tex) may be antigen-engaged T cell subsets (Fig S4c-d). However, the association of TMB, HRD score, and COSMIC mutational signatures with CD4_C03 (Tex) and CD8_C03 (Tex) is observed at the patient level rather than site level within individual patients (Fig S4a-b). In addition, we also constructed multi-region evolutionary trees based on somatic single-nucleotide variants (SNV) and structural variants (SV) across tumor sites (Fig S4e). Compared to PBMC, spatial genomic heterogeneity among tumors within individual patients is low especially between ovarian and omental metastatic tumors. Thus, spatial genomic features including TMB, HRD score and COSMIC mutational signature and evolution trajectory fail to explain the differences in T cell infiltration across different lesion sites within patients.

Tumor-specific but exhausted CD8 ⁺ cells preferentially infiltrate primary ovarian tumors, while non-tumor specific bystander cells are enriched in omentum metastases

To further investigate functional differences of CD8⁺ T cell clusters across locations, we first assessed transcriptional features of terminal exhaustion and effector memory signatures among CD8 T cell clusters by functional scores derived from previous reports^{24, 25} (Fig. 3a-b and S5a). As expected, CD8_C03 (Tex) and CD8_C05 (Tex.prol) had the highest terminal exhaustion characteristic (Fig. 3b), while CD8_C02 (Tex,trans), CD8_C04 (NK-like) and CD8_C06 (Teff) had features associated with effector and memory (Fig S5a). Given that exhausted T cells are frequently generated as a consequence of persistent antigen exposure²⁶, we next tested whether CD8_C03 (Tex) and CD8_C05 (Tex.prol) transcriptionally resemble neoantigen-reactive populations using a tumor-specific signature²⁷. Consistent with the concept that exhausted cells have undergone chronic antigen stimulation, the tumor-specific signature was significantly enriched in these two exhausted T subsets (Fig. 3c). Intra-tumoral T cells can also be CD39⁻ bystanders that recognize virus rather than tumor antigens^{28, 29}. Bystander signatures, including virus-specific and CD39⁻ CD69⁻ signatures, were dramatically increased in CD8_C02 (Tex,trans), CD8_C04 (NK-like) and CD8_C06 (Teff) that are enriched in omental tumors (Fig S5b-c). Collectively, CD8_C03 (Tex) and CD8_C05 (Tex.prol) that are enriched in ovarian tumors exhibited high exhaustion, tumor-specific score and low bystander score, whereas CD8_C02 (Tex,trans) and CD8_C04 (NK-like) that are enriched in omental tumors exhibited the opposite characteristics (Fig. 3d-e). Overall, tumor-specific signatures were strongly positively correlated with a terminal exhausted signature and were negatively associated with a bystander signature (Fig. 3f-g). Spatially, ovarian lesions had profoundly higher tumor-specific and terminal exhaustion scores than omental samples (Fig. 3f and h). Conversely, omental lesions exhibited markedly higher bystander scores (Fig. 3g and h). A heatmap of all signature scores showed the same distribution (Fig. 3i). Consistently, flow cytometry analysis confirmed more CD8⁺ tumor-specific T cells, expressing CD39⁺, were enriched in ovarian lesions, whereas more CD8⁺ bystander T cells enriched in omental lesions (Fig. 3j and S5d).

In addition, we reconstructed CD8 T cell antigen receptor (TCR) sequences from the scTCR-seq data. More than 70% of cells in all the tumor subsets had matched TCR information, with the exception of the NK-like subsets indicating limited drop out (Fig S5e). Given that peptide-MHC complex (pMHC) are recognized by specific TCRs, neoantigen and associated TCRs should be present in the same tissue³⁰. Accordingly, we first selected TCRs that had the same distribution as neoantigens and excluded neoantigen/TCR pairs identified in only one sample (Fig S5f-g). Peptide motifs in CDR3 are important for defining antigen specificity with a single antigen being recognized by multiple related TCRs. Consequently, clustering of CDR3 sequences is characteristic of an antigen-driven T cell response³¹. Thus, we calculated the pairwise similarity of CDR3 sequences between selected TCRs (same distribution as neoantigens) and randomly selected TCRs (Fig S5h). Selected CDR3 had higher similarity in each patient (Fig S5i). Finally, we calculated the proportion of cells corresponding to the selected and unselected TCRs in different clusters (Fig S5j). The proportion of selected cells were highest in CD8_C03 (Tex), followed by CD8_C06 (Teff) and lower in CD8_C02 (Tex,trans) (Fig S5k), which again supports the contention that CD8_C03 (Tex) represent a tumor specific cluster. We also compared bulk TCR data of each sample with three virus-specific TCR libraries (see Methods), with the results showing that omentum samples contained the highest proportion of virus-specific TCR (Fig S5l), further supporting their bystander T cell features.

More importantly, pseudotime analysis showed that omental TIL tends to be in early to mid-differentiation with continued transit, while TIL in ovarian tumors have limited transit consistent with terminally differentiated exhausted T cells (Fig. 3k-l). These results collectively indicated that the T cells infiltrating ovarian lesions were characterized by tumor-specific terminal exhaustion, while the T cells in the omentum were non-exhausted but also non-tumor specific.

Exhausted CD8 T cells enriched in primary ovarian tumors are clonally expanded

As noted above, we identified a proliferative CD8⁺ cluster (CD8_C05 (Tex.prol)) that highly expresses proliferation marker genes, such as TUBB, STMN1 and MKI67 that is enriched in ovarian tumors (Fig S6a-b). To better characterize this cluster, we used label transfer to interrogate the “second best” cluster for each proliferating cell³². Interestingly, the CD8_C05 (Tex.prol) cells were majorly regrouped into CD8_C03 (Tex) or CD8_C02 (Tex, trans) (Fig. 4a). Very few cells were reattributed to naïve or effector memory CD8⁺ T cell populations, suggesting that proliferating cells were transcriptionally closer to late-differentiated T exhaustion cells. Differential expression analyses in the regrouped CD8_C02 (Tex,trans) and CD8_C03 (Tex) cells after label transfer showed that proliferation-related pathways, including G2M checkpoint, mitotic spindle, DNA-repair, oxidative phosphorylation³³, and E2F targets³⁴ pathways were concurrently elevated in this subclass (Fig S6c).

Then we performed differential analysis of functional markers between CD8_C02 (Tex,trans), CD8_C03 (Tex), CD8_C05 (Tex.prol) and CD8_C06 (Teff). As expected, CD8_C02 (Tex,trans) showed increased GZMK, GZMM, GZMA, which are markers of transition status (Fig. 4b and S6d). Compared with CD8_C02 (Tex,trans), CD8_C05 (Tex.prol) had modestly increased levels of co-inhibition and co-stimulation genes (PDCD1, LAG3, TIGIT, CTLA4, TNFRSF4/9/14/18 and ICOS) and transcription factors (TOX, RBPJ, IRF9) (Fig. 4b and S5d), which are necessary and sufficient to induce major features of Tex cells³⁵. Of note, these co-inhibition, co-stimulation, and transcription factors were most highly expressed in CD8_C03 (Tex) consistent with exhaustion status (Fig. 4b and S5d). Cytotoxic markers (GZMB, PRF1, GNLY, GZMH) were low in both CD8_C05 (Tex.prol) and CD8_C02 (Tex,trans) indicating poor cytotoxic effector function. Notably, although weaker than CD8_C06 (Teff) cells, CD8_C03 (Tex) exhibited moderate GZMB, and PRF1 in the context of high FASLG and IFNG effector genes (Fig. 4b and S6d). Consistently, we observed the triple positive (CD8⁺PD-1⁺GZMB⁺) T cells in ovarian sample, but not in omental samples, by using opal multiplex IHC stains from site-matched FFPE sections, indicating the exist of CD8_C03 (Tex) cells exclusively in ovarian lesions and having modest cytotoxic activity (GZMB⁺) despite the expression of exhaustion markers (PD-1⁺) (Fig. 4c). Furthermore, we found that the CD8_C03 (Tex) gene signature score was associated with better overall survival, longer disease specific survival and better predicted response to ICB³⁶ in TCGA ovarian cancer patients (Fig S6e-g), which further suggests that the Tex population in ovarian cancer may have cytolytic activity and may contribute to response to ICB and improved outcomes.

To further explore the relationship between exhausted and cytotoxic functions, we calculated effector and exhaustion scores after label transfer. The positive correlation of exhaustion score and effector score in both proliferating and non-proliferating cells suggests that CD8 T cells in HGSOC concurrently exhibit cytotoxic capacity and exhaustion status (Fig. 4d). Proliferating T cells displayed lower effector and exhaustion scores, with a clone size that was much smaller than that of non-proliferating T exhausted cells (Fig. 4e). The STARTRAC-expansion index¹³ also showed that the CD8_C05 (Tex.prol) subclass had modest clonal expansion while the CD8_C03 (Tex) subclass had the highest degree of clonal expansion (Fig. 4f). When the dysfunction population was divided into decile according to exhaustion score, we found that as exhaustion scores increase, the proportion of proliferating cells first increased slightly, and then decreased sharply (Fig. 4g). The most exhausted cells completely lost proliferative ability (Fig. 4g). These results together are consistent with the exhausted CD8 T cell subclass developing from an early differentiation state with high proliferative capacity. Remarkably, we found there was a higher proportion of proliferating cells in each interval in ovarian tumors than in omental tumors (Fig. 4g).

Exhausted CD8 T cells are a consequence of differentiation

We performed STARTRAC-transition analysis to reveal T cell state transitions among CD8 cells. As expected, the probability of the same TCR being present between CD8_C02 (Tex,trans), CD8_C03 (Tex) and CD8_C05 (Tex.prol) was markedly higher compared to other clusters, indicating their considerable developmental state transitions exist across them (Fig. 5a-d). A UMAP of representative clonal sharing among CD8_C02 (Tex,trans), CD8_C03 (Tex) and CD8_C05 (Tex.prol) is shown in Fig. 5e. To further investigate clonal sharing among CD8_C02 (Tex,trans), CD8_C03 (Tex) and CD8_C05 (Tex.prol), we selected the top 30 clonal TCRs shared between CD8_C02 (Tex,trans) and CD8_C03 (Tex) clusters with or without proliferative status (CD8_C05 (Tex.prol)), and calculated the proportion of clonotype in each subclass. Interestingly, most of the top shared clones across the three subclasses were most frequently expressed as CD8_C03 (Tex), especially the TOP10 clones to the left of the dotted line (Fig. 5f). As the proliferative cells decreased (CD8_C05 (Tex.prol)), T cells in these clone types tend to be more in CD8_C02 (Tex,trans) status (Fig. 5f, left). In particular, the vast majority of TCR clones shared between the non-proliferating CD8_C02 (Tex,trans) and CD8_C03 (Tex) clusters were in the CD8_C02 (Tex,trans) cluster (Fig. 5f, right). On the whole, these results further support that the exhausted CD8 T cells develop following proliferation and clonal expansion. Label transfer of CD8_C02 (Tex,trans) cells showed that many of these cells were regrouped into CD8_C03 (Tex) (Fig. 5g-h) supporting the concept that CD8_C02 (Tex,trans) are transiting to the CD8_C03 (Tex) subclass.

The results presented thus far are consistent with the hypothesis that transition between CD8_C02 (Tex,trans) and CD8_C03 (Tex) clusters occurs while cells are proliferating. To test this possibility, we measured the frequency of proliferative cells among clones shared between the CD8_C02 (Tex,trans) and CD8_C03 (Tex) clusters across different sites (Fig. 5i-j) using CD8_C04 (NK-like) as a comparator. As shown in Fig. 5i, clone sharing predominantly occurred between CD8_C02 (Tex,trans) and CD8_C03 (Tex), with most of these shared clones also being present in proliferating cells. Specifically, compared with omental and other tumor sites, primary ovarian sites had a higher frequency of state transitions driven by proliferation between CD8_C02 (Tex,trans) and CD8_C03 (Tex) (Fig. 5j). Together this suggests that terminal exhaustion T cell differentiation preferentially occurs in primary ovarian sites.

CD4 Treg suppress the immune microenvironment in primary ovarian tumor sites

For CD4⁺ T cells (see Fig. 6a for a UMAP of CD4 T cells), CD4_C02 (Treg) and CD4_C03 (Tex) were enriched in ovarian tumors (Fig. 2e-f), while naïve, memory and transition functional state clusters were mainly present in omental tumors (Fig. 2e-f). We next assessed the expression of tumor-specific and bystander gene signatures in the CD4⁺ clusters (Fig. 6b-c). Notably, tumor-specific signature was significantly enriched in CD4_C03 (Tex), followed by CD4_C02 (Treg) population (Fig. 6b), while bystander signature was enriched in other naïve and effector/memory clusters, including CD4_C01/C05/C06 (Fig. 6c). Similar to CD8⁺ T cells, CD4⁺ cells in primary ovarian tumors displayed the highest tumor-specific and terminal exhaustion scores (Fig. 6d). Again, similar to CD8⁺ T cells, CD4⁺ T cells in omental sites exhibited the highest bystander score (Fig. 6e). Furthermore, CD4⁺, similar to CD8⁺, exhausted clusters expressed co-inhibitory and co-stimulatory receptor genes, including TIGIT, HAVCR2, CTLA4, PDCD1, and TNFRSF14 (Fig. 6f). There were differences with for example the co-stimulatory receptors TNFRSF4/18 and the co-inhibitory receptor LAG3 being highly expressed in CD4 Tex cluster, while TNFRSF9 was enriched in the CD8 Tex cluster (Fig. 6f). Of note, unlike CD8 Tex cells, almost all cytotoxic makers, including GZMA, GZMB, PRF1, GZMK, GNLY and CCL5, were absent in CD4_C03 (Tex), indicating a lack of cytotoxic activity (Fig. 6f).

A CD4 T cell cluster with proliferation characteristics expressed MKI67 and FOXP3 (Fig. 2b and S7a). Unlike proliferative CD8_C05 (Tex.prol), label transfer of CD4_C04 (Treg.prol) showed that these cells are exclusively related to CD4_C02 (Treg) but not exhausted CD4_C03 (Tex) cells (Fig. 6g). So, in addition to the lack of cytotoxicity noted above (Fig. 6f), CD4_C03 (Tex) did not exhibit proliferative capacity which was further supported by the relatively small clone size compared to exhausted CD8 T cells (CD8_C03) (Fig S7b). TCR similarity analysis by STARTRAC-transition showed that CD4_C04 (Treg.prol) shared TCRs with CD4_C02 (Treg) rather than CD4_C03 (Tex) (Fig. 6h, Fig S7c). Moreover, the transition between CD4_C02 (Treg) and CD4_C04 (Treg.prol) mainly occurred in ovarian tumors, represented by the green line (Fig S7d). Monocle 2 reconstructed a trajectory capturing the progression of CD4 reprogramming with a root at the highest naïve state (CD4_C07) and ending with two termini (Treg (CD4_C02 and CD4_C04) and Tex (CD4_C03)) corresponding to two distinct reprogramming outcomes (Fig. 6i). More importantly, while the terminal differentiated T cell clusters were enriched in ovarian tumors, early differentiated T cells were more frequent in omental tumors (Fig. 6j). Meanwhile, we computed the Treg score for each cell in Treg cells and calculated the proportion of proliferating cells in each score interval (Fig. 6k). Within the Treg cell pool in ovarian tumors, as the Treg score increased, the proportion of proliferating cells decreases sharply (Fig. 6k). More importantly, the proportion of proliferating cells in each interval is higher in ovarian tumors than that in omental tumors (Fig. 6k).

CD4 T cells can support effective anti-tumor CD8 function, but their cross-talk within the TME is not well characterized. To investigate molecular links underlying the intercellular communication of CD4⁺ and CD8⁺ T cell in HGSOC, CellphoneDB analysis³⁷ was used to identify molecular interactions between ligand-receptor pairs and major cell types in order to construct cellular communication networks. We found that interactions between Treg clusters, including CD4_C02 (Treg) and CD4_C04 (Treg.prol), and CD8 dysfunctional clusters, such as CD8_C03/05/07 rather than CD8_C01/02/04/06/09 non-dysfunctional subsets were commonly observed (Fig S7e). We subsequently analyzed detailed reciprocal connections between CD4_C02 (Treg) and all CD8 populations and identified markedly different ligand-receptor pairs between ovarian and omental tumors (Fig S7f). Notably, the KLRC1-HLA-E axis, a novel checkpoint in the TME³⁸ was exclusively enriched in ovarian tumors, whereas ICAM1/ICAM2 that has been characterized as a site for the cellular entry of human rhinovirus³⁹ and production of proinflammatory effects⁴⁰ was enriched in omental tumors which is in line with the bystander characteristics of omental tumors.

Inherent TME characteristics contribute to spatial differences of TIL status

To explore mechanisms underlying differences in infiltration of the Tex classes in tumor lesions, we performed pairwise STARTRAC-migration analysis of CD8_02/03/05 clusters between different lesions. We did not find evidence for T cell clusters in omental or other tumors preferentially migrating to ovarian tumors or vice versa (Fig. 7a). Moreover, migration of T cells between blood and different tumor lesions was extremely low with no evidence for preference for different tumor sites (Fig S8a). Therefore, spatial specific migration of individual T cell clusters is limited or absent. We subsequently analyzed the top 10 TCR clones per cluster in blood (Fig. 7b, top) or in tumors (Ovarian and Omental tumors) (Fig. 7b, bottom) for potentially transcriptional reprogramming between blood and different tumor foci. Notably, the top 10 TCR clonotypes from CD8_C03 (Tex) and C05 (Tex.prol) exhausted clusters were not observed in blood, consistent with these clones expanding intratumorally. Together the data argue that the preferential infiltration of CD8_C03 (Tex) and C05 (Tex.prol) in ovarian tumors is not due to migration from blood or other tumor sites.

Intra-tumoral T cell dysfunction has recently been suggested to be associated with reactivity to tumor antigens^{12, 41}. Consistent with this concept, TMB correlated with the proportion of dysfunctional T cells in primary ovarian tumors (Fig S3c-d), the differentiation process being associated with neoantigen recognition. We computed the CDR3 sequence similarity to investigate whether this differentiation process is antigen driven (not all patients shown, Fig S8b). Shared CDR3 sequences between transition states (CD8_C02) and exhausted states (CD8_C03 (Tex) and CD8_C05 (Tex.prol) were significantly elevated compared to unshared CDR3 sequences in different patients but not in different tumor regions (Fig. 7c and S8c). This suggests that while neoantigen may drive differentiation towards exhausted states, this does not explain the differences in exhausted T cells between different tumor sites.

We next performed differential analysis of signaling pathways between primary ovarian tumors and omental tumors on bulk transcriptomic profiles. Compared with omental tumors, proliferation-related pathways (G2M checkpoint, mitotic spindle, DNA-repair, oxidative phosphorylation, and E2F targets) and interferon signaling were concurrently increased in ovarian tumors (Fig. 7d). Consistently, these pathways were enriched in total (Fig. 7e, left), CD8⁺ (Fig. 7e, middle) or CD4⁺ T cells (Fig. 7e right) in primary ovarian tumors compared to omental tumors, indicating inherent TME characteristics contribute to spatial differences of TIL status. Notably, proliferation-related pathways, oxidative phosphorylation, glycolysis were all associated with T cell proliferation and function.

In contrast, consistent with the decreased interferon signaling in omental metastasis, MHC-I in tumor area detected by IHC was lower in omental metastasis (Fig. 7f and S8d). To further validate the differential expression of MHC-I in tumor cells of different lesions, we collected an additional 6 patients (36 samples, including 13 ovarian, 8 omental, 7 ascites and 8 other metastasis lesions) to perform scRNA-seq using the BGI droplet platform (Fig S8e, Supplementary table 2, see Methods). After quality control (see Methods), a total of 158,620 cells were available for analysis (Supplementary table 2). Using UMAP, we identified 7 stable clusters, including DCs, plasma, macrophages, T, endothelial, fibroblast and epithelial cells, each with unique signature genes (Figs S8e-f). In consonance with the differential expression of MHC-I detected by IHC staining, genes encoding for MHC-I processing and presentation were lowest in omental lesions, while HLA-B/C/E/F/G were highly expressed in ovarian lesions and HLA-A, CALR, PSMB6/8/9 and so on were highest in tumor cells of other lesions (Fig. 7g). As interferon increases antigen presentation⁴² and MHC-I reflecting antigen presentation ability and provides a marker of inflamed T cell infiltration ¹¹, the results suggested that omental tumors have lower antigen presentation ability.

Meanwhile, IHC staining of T cells exhibited that both CD4 and CD8 T cells were preferentially located in stroma than in tumor of omental tumors (Fig. 7h-i and S8g), indicating that most of the T cells in omental masses are excluded from contact with tumor cells.

T cells represent a major contributor to anti-tumor activity, a concept that is supported by the observation that intratumoral TIL are associated with an improved outcome in multiple diseases including ovarian cancer. The major components of the intratumoral T cell compartment include naïve, effector, memory, Treg and exhausted or dysfunctional T cells⁴³. To explore potential mechanisms underlying the limited response to immune checkpoint blockade in ovarian cancer, we used scRNA-Seq and TCR sequence analysis to determine immune contexture across different tumor sites and across different ovarian cancer patients. We supplemented these platforms with IHC analysis to provide spatial analysis. Together, this study provided a detailed analysis of the immune landscape across different lesions (Fig. 8d). Importantly, we found that the immune contexture in different tumor sites and in particular the two most common sites of ovarian cancer, the ovary and the omentum were markedly different. Ovarian tumors were characterized by an immunosuppressive environment consisting of Tregs and three different populations of exhausted CD8⁺ T cells as well as an exhausted CD4⁺ T cell population that likely acquired the exhausted phenotypes through interaction with tumor antigens in the local ovarian ecosystem. In contrast, TIL in omental lesions appear to consist primarily of non-tumor specific bystander cells with little evidence for response to tumor specific antigens. Differences in tumor mutation burden or in tissue specific immune cell migration do not appear to underly the diversity of TIL lineages in ovarian and omental lesions. Decreased MHC-I levels and antigen presentation could contribute to the low levels of exhausted of T cells and the decreased differentiation of T cells in omental tumors. While the exact underlying mechanisms remain to elucidated, decreased MHC-I antigen presentation and interferon signaling, oxidative phosphorylation and failure of T cell infiltration into omental tumors may contribute to lack of tumor specific T cells in omental metastasis and thus immune evasion⁴⁴.

The exhausted T cell state in ovarian tumors is likely a consequence of antigen stimulation leading to effector T cells eventually becoming exhausted due to prolonged antigen stimulation. We provide evidence for transition between the three types of exhausted T cells in ovarian tumors. Further, the exhausted T cells retain a number of markers that suggest that they could retain some degree of T cell killing activity. This may contribute to an elevated CD8_C03 (Tex) terminal exhausted signature score being associated with better prognosis in TCGA ovarian cohort. This may not be unique to ovarian cancer as Zhang et al. ⁴⁵ reported that CD8-CXCL13 and CD4-CXCL13 T cells, that are proposed to represent exhausted T cells, predict effective responses to PD-L1 blockade in breast cancer. However, other studies suggest that dysfunctional T cells can no longer be reversed and activated by PD-1 therapy⁴⁶.

Recently, Luca et al. developed a machine learning based algorithm, the EcoTyper, to deconvolve cell states and ecotypes⁴⁷, which identified two T cell associated carcinoma ecotypes (CE) across many tumor types, CE9 and CE10, wherein CE9-T cells express activation and exhaustion markers, similar with our CD8_C03 (Tex) cluster, and CE10-T cells that express GZMK and other naïve and memory markers, similar with our CD8_C02 (Tex,trans) cluster. In agreement with our study, CE9-T cells, characterized by higher immunoreactivity, preferentially infiltrated tumors compared to CE10-T cells and were strongly associated with longer overall survival. Additionally, characteristics of T cells identified in our studies are also recapitulated in a recent report which showed that antigen presentation gene sets, IFN gene sets and oxidative phosphorylation are enriched in infiltrated compared to excluded tumor cells⁴⁸. Similar to the CD8_C03 (Tex) subset in our study, they defined a CD8⁺ GZMB T subpopulation enriched in T cell infiltrated tumors that simultaneously exhibited exhaustion and cytotoxic characteristic, such as PRF1, GZMB, LAG3, CTLA4, PDCD1 and HAVCR2. Similar to the CD8_C02 (Tex,trans) cluster in our study, they identified a CD8⁺ GZMK T cell subpopulation that also lacks CD39 and thus likely represent a bystander population are enriched in stroma and are likely tumor excluded. Interestingly consistent with this concept, a GZMK/CD8⁺ ratio, which may represent a bystander signature, was significantly associated with shorter PFS. Both of these studies were based on single site sampling and thus did not observe the spatial heterogeneity of the immune contexture in ovarian cancer that may contribute to the limited response to immune therapy in HGSOC. Critically, the number of ovarian cancer samples analyzed by single cell sequencing was limited and thus this study combined with the recently published data provides an extended dataset that will greatly enhance our understanding of immune contexture across lesions in ovarian cancer and potentially contribute to development of effective immune therapy approaches in ovarian cancer.

In line with previous reports^{12, 49}, our data show that the dysfunctional T cell populations in ovarian cancer do not form a discrete cell population but rather develop from a precursor state with proliferative capacity. As cells differentiate into an exhausted state, they lose proliferative capacity. The proliferating precursor population has evidence for replication stress, high DNA repair capacity and oxidative phosphorylation, properties that have been observed in other tumor lineages^{12, 50}.

Although we used multiple complementary approaches (including IHC, genomic, bulk and single-cell transcriptional and TCR data), these approaches are mainly based on computational inference from static molecular snapshots. However, the derivation of the dysfunctional CD8 T cell state is likely the consequence of a dynamic process that occurs during tumor development. Indeed, as all of the tumors in this analysis were late stage with extensive spread, the “molecular snapshot” likely represents an immune status that is permissive for tumor growth potentially contributing to the extensive exhaustion states. To fully elucidate, the underlying mechanism both a dynamical study and an analysis of tumors at different stages of development will likely be needed. Furthermore, the status of these clusters was inferred by the expression of marker genes rather than by functional assays. Future studies incorporating lineage tracing, and single cell spatially resolved analysis will be needed to elucidate underlying mechanisms. One key clinical question will be to determine how to convert the immunosuppressive and exhausted environment to one that favors tumor clearance. In particular it will be important to determine whether the exhausted T cell state can be reversed to a functional state or whether the exhausted T cells are in an irreversible terminal state or trajectory. If this is the case, effective ovarian cancer immunotherapy may require use of modified T cells, such as CAR-T and TCR-T, combined with ICB. The lack of tumor reactive cells in omental tumors will likely require different approaches to induce immune engagement that in ovarian tumors. Indeed, the marked different in the immune contexture in ovarian and omental sites may be the major reason for failure of current immunotherapy approaches. Approaches that are effective in ovarian tumors, may not have significant activity in omental tumors and vice versa.

Our results, including trajectory analyses, TCR sharing, and cross-tissue comparisons, are most consistent with the final model: 1. Ovarian lesions have a tumor immunosuppressive environment with a high proportion of exhausted T cells and Treg. 2. The majority of tumor-specific TILs in ovarian lesions are exhausted as a consequence with interaction with tumor antigens. 3. TIL in omental lesions were primarily tumor non-specific and non-exhausted. Moreover, the decreased MHC-I antigen presentation and failure of T cell infiltration into tumors maybe associated with immune evasion of omental metastases. These results deepen our understanding of the poor response to ICB therapy in ovarian cancer while concurrently providing information that could improve our ability to engage the immune system in ovarian cancer.

Clinical sample collection

This study was reviewed and approved by the Institutional Review Board of Tongji Hospital, Tongji Medical College, Huazhong University of Science and Technology (TJ-IRB20190320). All the enrolled patients signed an informed consent form, and all the blood samples were collected using the rest of the standard diagnostic tests, with no burden to the patients.

Nine primary, untreated HGSOC patients who were pathologically diagnosed were enrolled in this study. Their ages ranged from 25 to 70 years old. For patients OV001, OV002 and OV003, their tumor tissues across all sites and PBMC isolated from blood were collected for FFPE, WGS, RNA-seq, bulk TCR-seq. For patients OV004, OV005, OV006, OV008, OV009 and OV010, their PBMC isolated from blood and tumor tissues across all sites were obtained for above sequencing and tissue dissociation to sort and obtain CD45⁺ CD3⁺ single cell suspension for scRNA-seq (Supplementary Table 1). Another 12 primary, untreated HGSOC patients who were pathologically diagnosed were enrolled for validation in this study. of which, the tumor tissues (ovarian and omental samples) of 5 patients were collected for flow cytometry analysis, whereas the tumor tissues of the other 7 patients (32 samples, including 9 ovarian, 6 omental, 6 ascites and 11 other metastasis lesions ) were dissociated for performing scRNA-seq (Supplementary Table 2).

Single cell collection

Peripheral blood mononuclear cells (PBMCs) were isolated from fresh peripheral blood by Ficoll-Paque Plus medium (GE Healthcare) according to the manufacturer’s instructions. Briefly, 4ml of fresh peripheral blood was collected during surgery in EDTA anticoagulant tubes and mixed with Ca/Mg-free PBS 1:1, then gently slowly layered onto 8ml Ficoll. After centrifugation, lymphocyte cells remained at the medium layer between plasma and Ficoll and were carefully transferred to a new tube, red blood cell lysis (Solarbio) was performed as appropriated, and then washed twice with PBS. Cell pellet was resuspended with sorting buffer (PBS with 0.5% Bovine Serum Albumin)

Fresh tumor tissues were cut into approximately 1-mm³ pieces and single-cell suspension was obtained by Tissue Dissociation Kits (Miltenyi, 130110201) together with the gentleMACS™ Dissociators and gentleMACS C tubes (Miltenyi, 130093237) according to the protocols. Briefly, tissue pieces were mixed using 5ml enzyme mix (4.7ml RPMI1640 + 200ul enzyme H + 100ul enzyme R + 25ul enzyme A) per C tube. After running the gentleMACS program h_tumor_01, incubate sample for 30min at 37℃ with continuous rotation using rotator. And again run the program h_tumor_01 and incubate sample for 30min at 37℃ using rotator. Finally run the program h_tumor_01 and collected the cell pellet to resuspend and filter through a 40µm cell-strainer ⁵¹ until uniform cell suspensions were obtained. Then the pelleted cells were suspended in red blood cell lysis buffer and washed twice, resuspended in sorting buffer.

Single cell sorting and scRNA library construction sequencing by 10x genomics

Based on FACS analysis, T cells (CD45⁺CD3⁺, BD Biosciences, 340943, 130-113-700) were sorted into tubes containing 0.5% BSA-PBS and stained with0.4% Trypan blue and examined by microscope. When the viability of cells was higher than 80%, use ChromiumTM Controller and ChromiumTM Single Cell 5’ Reagent Version 2 Kit (10x Genomics, Pleasanton, CA) for library construction experiments. In short, GemCode Technology was used to encapsulated sorted cells, reagents and Gel Beads containing barcoded oligonucleotides into nanoliter-sized GEMs. Lysis and barcoded reverse transcription of polyadenylated mRNA from single cells were performed within each GEM. Post RT-GEMs were cleaned up and cDNA were amplified. cDNA was fragmented and repaired at the end of the fragments, and an A-tail was added to the 5’ end. The adaptors were ligated to fragments which were double sided SPRI selected. After sample index PCR, another double sided SPRI selecting was performed. The final library was quality and quantitated using real-time quantitative PCR (TaqMan Probe). The final products were sequenced using the Xten-PE151 platform (BGIShenzhen, China).

TCR V(D)J sequencing

According to the manufacturer’s protocol (10x Genomics), the Chromium Single-Cell V(D)J Enrichment kit was used to enrich the full-length TCR V(D)J segments from amplified cDNA from 5’ libraries via PCR amplification.

Bulk DNA and RNA isolation and sequencing

Genomic DNA of peripheral blood and tissue samples were extracted using the QIAamp DNA Mini Kit (QIAGEN) according to the manufacturer’s instruction. Use the Qubit HsDNA Kit (invitrogen) to quantify the DNA concentrations and use agarose gel electrophoresis to evaluate the DNA quality. The exon library was constructed using the SureSelectXT target enrichment system for the illumine Double-End Multiplexed Sequencing Library Kit (Agilent). The samples were sequenced on the illumine Hiseq 4000 sequencer, and the paired-end read was 150 bp.

RNA of tumor samples was extracted by RNeasy Mini Kit (QIAGEN). The concentration of RNA was quantified by the NanoDrop instrument (Thermo) and the fragment analyzer (AATI) was used to evaluate the quality of RNA. Libraries were constructed using NEBNext Poly(A) mRNA Magnetic Isolation Module Kit (NEB) and NEBNext Ultra RNA Library Prep Kit for illumine Paired-end Multiplexed Sequencing Library (NEB). Samples were sequenced on the illumine Hiseq 4000 sequencer with 150bp paired-end reads.

Bulk TCR sequencing

RNA was extracted as described above and quantity were determined using Nanodrop. HTBI primers and Arm-PCR from iRepertoire were used to construct the libraries including PCR1 and PCR2, inclusively and semi-quantitatively. 5 cycles were used to amplify CDR3 fragments during the first round of PCR1, using the specific primers against each V and J genes. And in the second round, PCR was performed using universal primers.

PCR1: RNA reverse transcription and amplification of the T-cell receptor β CDR3 using the HTBI primers (Huntsville, Alabama, America) was carried out using Qiagen OneStep RT-PCR. The first round of PCR was performed using 200 ng of total RNA mixed with 4 µl random iRepertoire primers, 5 µl 5 × buffer, 1 µl dNTP mix, 0.25 µl RNasin (40 U/µl), and 1 µl enzyme mix, with nuclease-free water added to reach a total volume of 25 µl. After mixing and centrifugation, the reactions were transferred to a thermal cycler that carried out the following program: one cycle of 50°C for 40 min; one cycle of 95°C for 15 min; 15 cycles of denaturation at 94°C for 30 s, annealing at 60°C for 40 min, and extension for 30 s at 72°C; 10 cycles of denaturation at 94°C for 30 s, annealing and extension at 72°C for 2 min; and a final extension at 72°C for 10 min. The samples were then held at 4°C.

PCR2: A 2 µl sample of the PCR1 product was used as template for a second step of amplification following the addition of 5 µl communal primers, 25 µl Multiplex MM prepared using the Multiplex PCR Kit (Hilden, Nordrhein-Westfalen, Germany), and 18 µl nuclease-free water to reach a total volume of 50 µl. The reactions were then transferred to a thermal cycler that carried out the following program: one cycle of 95°C for 15 min; 40 cycles of denaturation at 94°C for 30 s, annealing at 55°C for 30 s and extension at 72°C for 30 s; and final extension at 72°C for 5 min. The samples were then held at 4°C. Size selection was used to purify 250-bp PCR products on magnetic beads (Agencourt No. A63882, Beckman, Beverly, MA, USA). After gel purification, the PCR product was subjected to HTS using the Hiseq PE151 platform.

scRNA library construction sequencing by BGI DNBelab C4 platform

Single-cell RNA-seq libraries were prepared using DNBelab C4 system as previously described⁵². Barcoded mRNA capture beads, droplet generation oil, and the single-cell suspension were loaded into the corresponding reservoirs on the chip for droplet generation. The droplets were gently removed from the collection vial and placed at room temperature for 20 minutes. Droplets were then broken and collected by the bead filter. The supernatant was removed, and the bead pellet was resuspended with 100µl RT mix. The mixture was then thermal cycled as follows: 42°C for 90 minutes, 10 cycles of 50°C for 2 minutes, 42°C for 2 minutes. The bead pellet was then resuspended in 200µl of exonuclease mix and incubated at 37°C for 45 minutes. Afterward, the PCR master mix was added to the beads pellet and thermal cycled as follows: 95°C for 3 minutes, 13 cycles of 98°C for 20 s, 58°C for 20 s, 72°C for 3 minutes, and finally 72°C for 5 minutes. Amplified cDNA was purified using 60µl of AMPure XP beads. The cDNA was subsequently fragmented to 400–600 bp with NEBNext dsDNA Fragmentase (New England Biolabs) according to the manufacturer’s protocol. Indexed sequencing libraries were constructed using the reagents in the C4 scRNA-seq kit following the steps: (1) post fragmentation size selection with AMPure XP beads; (2) end repair and A-tailing; (3) adapter ligation; (4) post ligation purification with AMPure XP beads; (5) sample index PCR and size selection with AMPure XP beads. The barcode sequencing libraries were quantified by Qubit (Invitrogen). The sequencing libraries were sequenced using the DIPSEQ T1 sequencer at the China National GeneBank. The read structure barcode 2, and 10 bp unique molecular identifier, and Read 2 containing 100 bases of transcript sequence, and a 10 bp sample index.

Cell clustering of BGI DNBelab C4 scRNA-seq dataset

Clustering analysis of the scRNA dataset was performed using Seurat (v4.0.4) and the R program, and the parameters were manually curated to portray an optimal classification of cell types with empirical knowledge. Specifically, low quality cells were filtered with fewer than 500 detected genes or above 6,000, as well as with higher than 20% mitochondrial counts in data preprocessing, and all query genes were guaranteed to be expressed in at least three cells for further use. The top 2,000 highly variable genes were then selected according to their mean-variance ratio on expression levels after log1p normalization. For downstream clustering and visualization, principal component analysis (PCA)-based dimension reduction was initially generated, and the first 30 PCs were extracted for subsequent Louvain clustering to define the cell types (the resolution was set to 0.5). The clustering result was finally characterized in a two-dimensional space using the UMAP technique, and the cell types were annotated by known biomarkers that were more highly expressed in a particular cluster (via FindAllMarkers function with default parameters).

Immunohistochemistry

The specimens were collected within 30 min after the tumor resection and fixed in formalin for 48hr. Paraffin-embedded tissues were subsequently cut into 4µm slides and mounted on glass slides. Tissues were subjected to deparaffinization and then rehydrated in 100%, 90%, 70% alcohol successively. Antigen was retrieved prior to antibody staining, and then endogenous peroxidase was inactivated by incubation in 3% H₂O₂ for 30min. After 10% normal goat serum blocking nonspecific sites for 1hr, 37℃, slides were stained with primary antibody overnight at 4℃ (anti-CD8 antibody, Maxim biotechnologies; anti-CD4 antibody, Maxim biotechnologies; anti-FAP, 1:250, abcam, ab207178; anti-MHC-I, 1:100, abcam, ab52922). Negative controls were treated identically, but with normal serum. After the sections were washed with PBS twice for 5min, the antigenic binding sites were visualized using secondary antibody (Beyotime, A0208 and A0216). After mounting, slides were scanned, and pictures were taken. The number of positive cells was counted in random five areas that included both tumor and stroma at a high-power field (HPF, 200x). TILs densities for each image were calculated by normalizing validated TIL counts by total area covered by tissue in the image (in units of cells/HPF). H-score analysis was performed on FAP and MHC-I IHC described by Melanie¹¹. The H-score was calculated by adding up the percentage of cells in each scoring category multiplied by the corresponding score using a semiquantitative five category grading system: 0, no staining; 1, 1–10% staining; 2, 11–50% staining; 3, 50–75% staining; and 4, >75% staining. resulting in scores on a scale of 0-400. Staining score was determined separately by two experts under the same conditions, while discordant scores were reevaluated by another expert.

Opal multiplex IHC

Deparaffinization of formalin fixed paraffin embedded (FFPE) sections was done through xylenes. Rehydration was done through decreasing graded alcohol. AR9 buffer (Perkin Elmer, AR900250ML) was used for retrieving antigens in a microwave oven and a hydrophobic pen was used to circle tissue sections. Before primary antibody incubation, tissue sections were blocked with blocking/antibody diluent (Akoya antibody diluent/block, ARD1001) for 30 min at RT. The tissue sections were incubated with primary antibodies for 60min at RT. Then washed in 1×TBST and incubated with secondary antibody (Perkin Elmer opal polymer HRP Ms + Rb ARH1001EA, 30min at RT). The HRP-conjugated secondary antibody polymer was detected using fluorescent tyramide signal amplification using Opal dyes 520, 540, 620 and 690 (Perkin Elmer FP1487001KT, FP1494001KT, FP1495001KT, FP1497001KT) for 10min at RT. The covalent tyramide reaction was followed by heat induced stripping of the primary/secondary antibody complex using Perkin Elmer AR9 buffer (AR900250ML) at 100°C for 15min preceding the next cycle (each cycle for each marker). After 4 sequential rounds of staining, sections were stained with DAPI (Perkin Elmer, FP1490A) to visualize nuclei. Five color multiplex-stained slides were imaged using the Vectra Multispectral Imaging System version 2 (Perkin Elmer). Scanning was performed at 20×(200× final magnification). Filter cubes used for multispectral imaging were DAPI, FITC, Cy3, Texas Red and Cy5. A spectral library containing the emitted spectral peaks of the fluorophores in this study was created using the Inform analysis software (Perkin Elmer). Using multispectral images from single-stained slides for each marker, the spectral library was used to separate each multispectral cube into individual components allowing for identification of the five marker channels of interest using Inform 2.4 image analysis software. anti-CD8 (MAB-0021, 6ml volume, Maxim Biotechnology, 1:5, 60min, opal540), anti-PD-1 (ab137132, Abcam, 1:1500, 60min, opal620), anti-GZMB (ab4059, Abcam, 1:1500, 60min, opal690) and anti-pan-cytokeratin (pan-CK) (RAB-0050, 6ml volume, Maxim Biotechnology, 1:10, 60min, opal520) respectively at RT sequentially.

Gene expression quantification

Raw paired-end reads are filtered to remove adapter sequence using pipeline in-house. And then align to reference genome hg19 by STAR⁵³. RSEM⁵⁴ was used to quantify gene expression based on uniquely mapped reads. GENECODE V19 is used for annotation.

Somatic mutation calling

Raw reads are pre-processed to remove adapter sequences and low-quality reads. The processed clean reads are mapped to hg19 using BWA⁵⁵ with the default parameter. Picard are used to mark duplicates; GATK4 are employed for base quality correction and realignments. Mutect2⁵⁶ are used for somatic SNV/Indel calling. Mutations were filtered with supported reads ≥ 4 (≤ 2) and coverage ≥ 10 in tumor and (normal tissue), whereas indels were filtered with supported reads ≥ 5 (≤ 1) and coverage ≥ 10 in tumor and (normal tissue. Moreover, somatic mutations and indels were annotated by Oncotator⁵⁷.

Copy number calling and Tumor purity estimation

We estimated copy number profiling over 200bp bins using Patchwork⁵⁸, and then calculated the normalized ratio of standardized, average depth between normal tissue and tumor tissue. Fifty bins are further merged into 10kb windows. Segmentation performed all the 10kb windows. After that, tumor ploidy and purity were quantified using Patchwork based on the VAF of each somatic SNV and the copy number status of each segment.

BRCA germline variants

Germline indels of BRCA1/2 were called by SvABA⁵⁹ using default parameters. Germline deleterious SNVs of BRCA1/2 were selected with annotating be pathogenic in the ClinVar database.

Mutational process

We applied the R package deconstructSigs⁶⁰ to estimate the contributions of 30 mutational signatures documented by the COSMIC⁶¹ for each sample. The 30 signatures are annotated as mutagenesis forms based on COMSIC.

Somatic SV detection

We applied SvABA⁵⁹ to predict somatic SVs and their breakpoints using the suggested parameters. SV with Q value less than ten is filtered. SVs that are marked by TSI-L are omitted.

Chromothripsis

We used four criteria to infer chromothripsis proposed by Campbell⁶²: A). the number four types of SV type (tail-to-tail, head-to-head, head-to-tail, tail-to-head) are comparable. B). the number of segments involved in chromothripsis is more than 5. C). the copy number oscillated between 2 or 3 copy number states. D). there are interspersed LOH within affected regions.

BFB detection

We inferred BFB events by detecting fold-back inversion and telomere loss which is introduced by Campbell⁶³. Fold-back inversions were detected based on three criteria: 1) the single inversions were without reciprocal support-read clusters, 2) the inversion associated with a copy number change (q < 0.001), and 3) the two ends of the breakpoints had to be separated by 30 kb.

Neoantigen identification

The HLA type for each sample was inferred using HLA-VBSeq v2⁶⁴ by optimizing read alignments to HLA allele sequences and abundance of reads on HLA alleles by variational Bayesian inference. MHC binding affinities are inferred as IC50 values for each peptide sequence and patient HLA type. All mutant peptide sequences considered to be neoantigens meet a standard cut-off: the IC50 of mutant peptide < 500nM and the IC50 of the wild peptide > 500nM. NetMHCpan4.0⁶⁵ predicted Peptide-MHC class I binding affinity, while NetMHCIIpan-3.0⁶⁶ is applied to identify peptides that bind to MHC-II molecules. Neoantigens with at least three RNA-seq reads contained the mutated base was considered to be expressed.

Differentially expressed genes between Omen and Ovary

We used a rank-sum test (Python package stats) to compare gene expression between samples of omen and ovary groups. P-values are adjusted using FDR, and genes with FDR < 0.1 are regarded as significant. GSEA is used for pathway enrichment.

Bulk TCR-seq data analysis

Raw sequencing data were processed by the tool IMonitor (v1.4.1)⁶⁷. Briefly, the raw paired-end (PE) reads were merged to one sequence by the overlapped region. Low-quality sequences were filtered out. The clean sequences were aligned to reference that including V, D and J germline sequences (www.imgt.org). Originated V, D and J genes were determined for sequences and CDR3 regions were identified. Sequencing errors in CDR3 sequences were corrected according to the base sequencing quality and CDR3 frequency. Nucleotide CDR3 sequences were translated to amino acid sequences. Finally, multiple diversity indexes were calculated, and Figs were generated to display the TCR repertoire.

TCR repertoire annotated by disease associated TCR database

Disease-associated TCR sequences consisted by three published databases, including VDJdb^{68, 69}, McPAS-TCR⁷⁰ and TBAdb⁷¹. The three databases includes 116,875 records of TCR β sequences covering 43 sorts of infectious diseases. Only high-quality records (VDJdb: Score > = 2; McPAS-TCR: Antigen.identification.method < = 2.5; TBAdb: Grade > = 4) were selected for further analysis. For each TCR sample, CDR3 amino acid sequences were compared with the CDR3s in three databases by Levinshtein distance. The CDR3 sequence in the sample was supposed to be related to the disease if the Levinshtein distance < = 1 between CDR3 in sample and CDR3 associated with a disease in database. At last, the proportion of disease-related CDR3s, the number of disease-related CDR3s divided by total number of CDR3 in a sample, was calculated. For bulk TCR-seq sample, top 9000 CDR3 sequences ranked by frequency。

Inferring neoantigen associated TCRs

Both MHC I and MHC II restricted neoantigens peptides were predicted from the data of WGS. As the peptide-MHC complex (pMHC) can be specifically recognized by TCRs, neoantigen and associated TCRs are supposed to present at the same regional tissue. Thus, according to the location consistency between them, associated TCRs could be inferred by the distribution of neoantigens. Specifically, for a patient, if the neoantigen was identified in multiple samples, the TCRs that were observed in all the same samples and were not observed in other samples were regarded as the neoantigen associated TCRs. The neoantigens identified in only one sample were excluded. MHC I restricted neoantigens were used to find associated CD8 TCRs while MHC II restricted neoantigens were used for CD4 TCRs. To reduce potential error of TCR, the TCR with at least two cells was regarded as presence in the sample. Additionally, for each patient, the equal number of TCRs as control were selected randomly from all remaining TCRs that didn’t include inferred neoantigen associated TCRs and were from at least two cells.

Pairwise similarity calculation between TCRs

Both α and β CDR3 amino acid sequences were used to calculate the CDR3 sequence was deconstructed into series of contiguous triplet amino acids, started from the first amino acid and with stride 1. For example, the length of CDR3 with 15 amino acids could be deconstructed into 13 triplets. The similarity is equal to the number of shared triplet amino acids divided by total number of triplets from pair TCRs:

$$similarity= \frac{{2*n}_{\alpha }+{2*n}_{\beta }}{{N}_{\alpha }+{N}_{\beta }}$$

where, ${n}_{\alpha }$ is the number of shared triplet amino acids from α CDR3, ${n}_{\beta }$is the number of shared triplet amino acids from β CDR3, ${N}_{\alpha }$is the total number of triplet amino acids from pairwise α CDR3s, ${N}_{\beta }$ is the total number of triplet amino acids from pairwise β CDR3s.

scRNA/TCR-seq data processing

Single-cell expression was analyzed using the Cell Ranger Single Cell Software Suite (v3.0.2, 10x Genomics) to perform quality control, sample demultiplexing, barcode processing, and single-cell 3′ and 5′ gene counting. Sequencing reads were aligned to the GRCh38 human reference genome. Seurat v4 (version 4.0.4) R package was used to analyze the scRNA-seq data. To integrate and embed single cells from different samples into a shared low-dimension space (Fig S2a), we utilized integrated analysis (CCA) by the Seurat function IntegrateData⁷².

Signature gene sets

Terminally exhausted CD8⁺ signature, T effector memory signature, CD39⁻CD69⁻ signature, and tumor specific CD8⁺ signature was got from previous studies^{24, 25, 27}. These gene sets were used as modules for the AddModuleScore function in Seurat.

Proliferation state definition

The average expression of known proliferation-related genes (ZWINT, E2F1, FEN1, FOXM1, H2AFZ, HMGB2, MCM2, MCM3, MCM4, MCM5, MCM6, MKI67, MYBL2, PCNA, PLK1, CCND1, AURKA, BUB1, TOP2A, TYMS, DEK, CCNB1 and CCNE1) was defined as the proliferation score⁷³.

Identification of signature genes and TCGA data analysis

We identified differentially expressed genes (DEGs) based on the FindAllMarkers function of seurat by using wilcox test⁷². The DEGs with FDR < 0.01 and log₂(FC) > 1 were selected as the signature genes of CD8_C03 terminal exhausted cluster. The TCGA bulk RNA-seq and clinical data were obtained from UCSC Xena (https://xenabrowser.net/datapages/). The mean value of the expression (log2(tpm + 0.001)) of the CD8_C03 signature genes (CD8A, CXCL13, DUSP4, LAG3, GZMB, CCL5, CCL3, CCL4, NKG7) was calculated as the signature score. Kaplan-Meier survival curves were used to show the survival differences between different groups (high group, greater than or equal to the median signature score, vs. low group, less than the median signature score). The R packages survival v3.2-13 and survminer v0.4.9 were used to perform all survival analyses. TIDE (Tumor Immune Dysfunction and Exclusion) was used to predict the immunotherapy responses as described in a previous study³⁶.

TCR analysis

TCR-seq data for each sample were processed using Cell Ranger software (v3.0.2), with the command “cellranger vdj” using the human reference genome GRCh38. To integrate TCR results with the gene expression data, the TCR-based analysis was performed only for cells that were identified as T cells. T cells with TCR information were used to perform the STARTRAC analysis as we previously described¹³.

Gene Set Enrichment Analysis

Different gene expression between T cells from ovarian (Ov) and T cells from omental (Om) were calculated based on the FindAllMarkers function of seurat by using wilcox test. Sorted (by log fold change) different expression gene list was used to perform the gene set enrichment analysis (GSEA) by using clusterProfiler (version3.18.0) package.

Ligand-receptor expression analysis

To analyze cell-cell interactions between clusters of interest, CellPhoneDB³⁷ (v2.1.1) was used to identify significant ligand-receptor pair in our data. We randomly selected 10% cells per cluster to perform the analysis. Potential ligand- receptor interactions were identified based on specific expression of a ligand by one T cell cluster and the corresponding receptor by another. The ligand-receptor expression analysis of cells from different lesion sites were performed separately.

Trajectory analysis

To compute pseudotime alignment of our transcriptomes, Monocle2 (v2 2.4.0) was used by using the first 30 PCs of the integrated matrix to preform preprocessing and UMAP reduction. DDRTree algorithm was then used to reconstruct the tree embedding.

Software versions

Data were collected using Cell Ranger software (10x Genomics) v3.0.2 and analyzed using R v.4.0.3, and the following packages and versions in R for analysis: Seurat v4.0.4, clustree v0.4.3, and cluster v2.1.2 two-dimensional gene expression maps, were generated using coordinates from the UMAP algorithm using the R package uwot v0.1.10 implementation. Figs were produced using the following packages and versions in R: ggplot v3.3.5, ComplexHeatmap v2.8.0, ggchicket v0.5.2, patchwork v1.1.1, circlize v0.4.13, ggtern v3.3.5, ggpubr v0.4.0, igraph v1.2.7, and RColorBrewer v1.1-2.

Statistical analyses

Python (v3.6) package sklearn is used to fit Gaussian Mixture Models (GMMs). p-values based on two groups are computed using python package stats. Plots are mainly based on matplotlib and seaborn. Paired t test was used to compare differences between two matched groups (Fig. 7d). Two-sided student t test or wilcox.test was used to compare differences between two groups of disease stage. If the multiple groups data followed a normal distribution, we used ANOVA test for multiple comparisons. IHC staining data was plotted and multiple compared by Tukey’s test using GraphPad Prism 8.0.2 software. Data are presented as means ± SEM and p<0.05 was considered significant. Correlation between groups was determined by Pearson correlation test. ANOVA was used to compare differences among multiple groups.

Acknowledgements

We thank the anonymous referees for their useful suggestions and all the enrolled patients for their dedication to science. This study is supported by a grant Nature and Science Foundation of China (81572569, 81874106, 81402163, 81974408), National Key R&D Program of China (2016YFC1303100) and Wuhan Municipal health Commission (WX18Q16).

Author contributions: C.S. and G.C. conceived of and supervised the project. B.Y. collected the samples and performed IHC staining. X.L., W.Z. and J.F. performed the bioinformatics analyses. B.Y. designed, carried out and interpreted the wet laboratory experiments. B.Y. wrote the manuscript. B.Y., W.L and J.Y. scanned the IHC Figs and calculated the score. E.G., X.L., Y.F., S.L., R.X., D.H., X.Q., F.L., Z.W. and T.Q. collected the patient’s clinical information and contributed to data processing and analyses. Q.Z., D.M., S.L. and G.B.M. provided expertise and feedback. G.B.M., C.S. and G.C. contributed to project management and provided valuable critical discussion.

Conflict of interests

GBM has licensed an HRD assay to Myriad Genetics and on Digital Spatial Profiling to Nanostring; is a SAB member/consultant with Amphista, AstraZeneca, Chrysallis Biotechnology, GSK, ImmunoMET, Ionis, Lilly, PDX Pharmaceuticals, Signalchem Lifesciences, Symphogen, Tarveda, Turbine, and Zentalis Pharmaceuticals.

Ferlay J, et al. Cancer incidence and mortality worldwide: sources, methods and major patterns in GLOBOCAN 2012. International journal of cancer 136, E359-E386 (2015).
Torre LA, et al. Ovarian cancer statistics, 2018. CA: a cancer journal for clinicians 68, 284–296 (2018).
Sharma P, Allison JP. The future of immune checkpoint therapy. Science 348, 56–61 (2015).
González-Martín A, Sánchez‐Lorenzo L. Immunotherapy with checkpoint inhibitors in patients with ovarian cancer: still promising? Cancer 125, 4616–4622 (2019).
Jayson GC, Kohn EC, Kitchener HC, Ledermann JA. Ovarian cancer. The Lancet 384, 1376–1388 (2014).
Ma X. The omentum, a niche for premetastatic ovarian cancer. Journal of Experimental Medicine 217, (2020).
McPherson A, et al. Divergent modes of clonal spread and intraperitoneal mixing in high-grade serous ovarian cancer. Nature genetics 48, 758–767 (2016).
Schwarz RF, et al. Spatial and temporal heterogeneity in high-grade serous ovarian cancer: a phylogenetic analysis. PLoS medicine 12, e1001789 (2015).
Ewing A, et al. Structural Variants at the BRCA1/2 Loci are a Common Source of Homologous Repair Deficiency in High-grade Serous Ovarian Carcinoma. (2021).
Network CGAR. Integrated genomic analyses of ovarian carcinoma. Nature 474, 609 (2011).
Desbois M, et al. Integrated digital pathology and transcriptome analysis identifies molecular mediators of T-cell exclusion in ovarian cancer. Nature communications 11, 1–14 (2020).
Li H, et al. Dysfunctional CD8 T cells form a proliferative, dynamically regulated compartment within human melanoma. Cell 176, 775–789. e718 (2019).
Zhang L, et al. Lineage tracking reveals dynamic relationships of T cells in colorectal cancer. Nature 564, 268–272 (2018).
Bongen E, Vallania F, Utz PJ, Khatri P. KLRD1-expressing natural killer cells predict influenza susceptibility. Genome medicine 10, 1–12 (2018).
Zhou L, Adrianto I, Wang J, Wu X, Datta I, Mi Q-S. Single-cell RNA-Seq analysis uncovers distinct functional human NKT cell sub-populations in peripheral blood. Frontiers in cell and developmental biology 8, 384 (2020).
Pizzolato G, et al. Single-cell RNA sequencing unveils the shared and the distinct cytotoxic hallmarks of human TCRVδ1 and TCRVδ2 γδ T lymphocytes. Proceedings of the National Academy of Sciences 116, 11906–11915 (2019).
Germano G, Amirouchene-Angelozzi N, Rospo G, Bardelli A. The clinical impact of the genomic landscape of mismatch repair–deficient cancers. Cancer discovery 8, 1518–1528 (2018).
McGranahan N, et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science 351, 1463–1469 (2016).
Pellegrino B, et al. Homologous recombination repair deficiency and the immune response in breast cancer: a literature review. Translational oncology 13, 410–422 (2020).
Choucair K, Morand S, Stanbery L, Edelman G, Dworkin L, Nemunaitis J. TMB: a promising immune-response biomarker, and potential spearhead in advancing targeted therapy trials. Cancer gene therapy 27, 841–853 (2020).
Lotan TL, et al. Homologous recombination deficiency (HRD) score in germline BRCA2-versus ATM-altered prostate cancer. Modern Pathology 34, 1185–1193 (2021).
Tate JG, et al. COSMIC: the catalogue of somatic mutations in cancer. Nucleic acids research 47, D941-D947 (2019).
Ghorani E, et al. The T cell differentiation landscape is shaped by tumour mutations in lung cancer. Nature cancer 1, 546–561 (2020).
Sade-Feldman M, et al. Defining T cell states associated with response to checkpoint immunotherapy in melanoma. Cell 175, 998–1013. e1020 (2018).
Van der Leun AM, Thommen DS, Schumacher TN. CD8 + T cell states in human cancer: insights from single-cell analysis. Nature Reviews Cancer 20, 218–232 (2020).
Wherry EJ, Kurachi M. Molecular and cellular insights into T cell exhaustion. Nature Reviews Immunology 15, 486–499 (2015).
Oliveira G, et al. Phenotype, specificity and avidity of antitumour CD8 + T cells in melanoma. Nature 596, 119–125 (2021).
Simoni Y, et al. Bystander CD8 + T cells are abundant and phenotypically distinct in human tumour infiltrates. Nature 557, 575–579 (2018).
Kim T-S, Shin E-C. The activation of bystander CD8 + T cells and their roles in viral infection. Experimental & molecular medicine 51, 1–9 (2019).
Zhang AW, et al. Interfaces of malignant and immunologic clonal dynamics in ovarian cancer. Cell 173, 1755–1769. e1722 (2018).
Joshi K, et al. Spatial heterogeneity of the T cell receptor repertoire reflects the mutational landscape in lung cancer. Nature medicine 25, 1549–1559 (2019).
Stuart T, et al. Comprehensive integration of single-cell data. Cell 177, 1888–1902. e1821 (2019).
Jones N, et al. Metabolic adaptation of human CD4 + and CD8 + T-cells to T-cell receptor-mediated stimulation. Frontiers in immunology 8, 1516 (2017).
Zhu JW, et al. E2F1 and E2F2 determine thresholds for antigen-induced T-cell proliferation and suppress tumorigenesis. Molecular and cellular biology 21, 8547–8564 (2001).
Khan O, et al. TOX transcriptionally and epigenetically programs CD8 + T cell exhaustion. Nature 571, 211–218 (2019).
Jiang P, et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nature medicine 24, 1550–1558 (2018).
Efremova M, Vento-Tormo M, Teichmann SA, Vento-Tormo R. CellPhoneDB: inferring cell–cell communication from combined expression of multi-subunit ligand–receptor complexes. Nature protocols 15, 1484–1506 (2020).
Borst L, van der Burg SH, van Hall T. The NKG2A–HLA-E axis as a novel checkpoint in the tumor microenvironment. Clinical Cancer Research 26, 5549–5556 (2020).
Abraham G, Colonno RJ. Many rhinovirus serotypes share the same cellular receptor. Journal of Virology 51, 340–345 (1984).
Etienne-Manneville S, Chaverot N, Strosberg AD, Couraud P-O. ICAM-1-coupled signaling pathways in astrocytes converge to cyclic AMP response element-binding protein phosphorylation and TNF-α secretion. The Journal of Immunology 163, 668–674 (1999).
Thommen DS, Schumacher TN. T cell dysfunction in cancer. Cancer cell 33, 547–562 (2018).
Hervas-Stubbs S, Perez-Gracia JL, Rouzaut A, Sanmamed MF, Le Bon A, Melero I. Direct effects of type I interferons on cells of the immune system. Clinical Cancer Research 17, 2619–2627 (2011).
Zheng C, et al. Landscape of infiltrating T cells in liver cancer revealed by single-cell sequencing. Cell 169, 1342–1356. e1316 (2017).
Dhatchinamoorthy K, Colbert JD, Rock KL. Cancer Immune Evasion Through Loss of MHC Class I Antigen Presentation. Frontiers in immunology 12, 469 (2021).
Qiu S, et al. A single-cell immune atlas of triple negative breast cancer reveals novel immune cell subsets. BioRxiv, 566968 (2019).
Blank CU, et al. Defining ‘T cell exhaustion’. Nature Reviews Immunology 19, 665–674 (2019).
Luca BA, et al. Atlas of clinically distinct cell states and ecosystems across human solid tumors. Cell, (2021).
Hornburg M, et al. Single-cell dissection of cellular components and interactions shaping the tumor immune phenotypes in ovarian cancer. Cancer Cell, (2021).
Azizi E, et al. Single-cell map of diverse immune phenotypes in the breast tumor microenvironment. Cell 174, 1293–1308. e1236 (2018).
Guo X, et al. Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing. Nature medicine 24, 978–985 (2018).
Parvathareddy SK, Siraj AK, Al-Badawi IA, Tulbah A, Al-Dayel F, Al-Kuraya KS. Differential expression of PD-L1 between primary and metastatic epithelial ovarian cancer and its clinico-pathological correlation. Scientific Reports 11, 1–9 (2021).
Liu C, et al. A portable and cost-effective microfluidic system for massively parallel single-cell transcriptome profiling. bioRxiv, 818450 (2019).
Dobin A, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 (2013).
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC bioinformatics 12, 323 (2011).
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754–1760 (2009).
Cibulskis K, et al. Sensitive detection of somatic point mutations in impure and heterogeneous cancer samples. Nature biotechnology 31, 213–219 (2013).
Ramos AH, et al. Oncotator: cancer variant annotation tool. Human mutation 36, E2423-2429 (2015).
Mayrhofer M, DiLorenzo S, Isaksson A. Patchwork: allele-specific copy number analysis of whole-genome sequenced tumor tissue. Genome biology 14, R24 (2013).
Wala JA, et al. SvABA: genome-wide detection of structural variants and indels by local assembly. Genome research 28, 581–591 (2018).
Rosenthal R, McGranahan N, Herrero J, Taylor BS, Swanton C. DeconstructSigs: delineating mutational processes in single tumors distinguishes DNA repair deficiencies and patterns of carcinoma evolution. Genome biology 17, 31 (2016).
Tate JG, et al. COSMIC: the Catalogue Of Somatic Mutations In Cancer. Nucleic acids research 47, D941-D947 (2019).
Korbel JO, Campbell PJ. Criteria for inference of chromothripsis in cancer genomes. Cell 152, 1226–1236 (2013).
Campbell PJ, et al. The patterns and dynamics of genomic instability in metastatic pancreatic cancer. Nature 467, 1109–1113 (2010).
Nariai N, et al. HLA-VBSeq: accurate HLA typing at full resolution from whole-genome sequencing data. In: BMC genomics). Springer (2015).
Jurtz V, Paul S, Andreatta M, Marcatili P, Peters B, Nielsen M. NetMHCpan-4.0: improved peptide–MHC class I interaction predictions integrating eluted ligand and peptide binding affinity data. The Journal of Immunology 199, 3360–3368 (2017).
Nielsen M, Andreatta M. NetMHCpan-3.0; improved prediction of binding to MHC class I molecules integrating information from multiple receptor and peptide length datasets. Genome medicine 8, 1–9 (2016).
Zhang W, et al. IMonitor: A Robust Pipeline for TCR and BCR Repertoire Analysis. Genetics 201, 459–472 (2015).
Shugay M, et al. VDJdb: a curated database of T-cell receptor sequences with known antigen specificity. Nucleic Acids Res 46, D419-D427 (2018).
Bagaev DV, et al. VDJdb in 2019: database extension, new analysis infrastructure and a T-cell receptor motif compendium. Nucleic Acids Res 48, D1057-D1062 (2020).
Tickotsky N, Sagiv T, Prilusky J, Shifrut E, Friedman N. McPAS-TCR: a manually curated catalogue of pathology-associated T cell receptor sequences. Bioinformatics 33, 2924–2929 (2017).
Sun C, et al. Correction to: A reactive oxygen species scoring system predicts cisplatin sensitivity and prognosis in ovarian cancer patients. BMC cancer 20, 210 (2020).
Hao Y, et al. Integrated analysis of multimodal single-cell data. Cell, (2021).
Whitfield ML, George LK, Grant GD, Perou CM. Common markers of proliferation. Nature Reviews Cancer 6, 99–106 (2006).

Yes there is potential Competing Interest. GBM has licensed an HRD assay to Myriad Genetics and on Digital Spatial Profiling to Nanostring; is a SAB member/consultant with Amphista, AstraZeneca, Chrysallis Biotechnology, GSK, ImmunoMET, Ionis, Lilly, PDX Pharmaceuticals, Signalchem Lifesciences, Symphogen, Tarveda, Turbine, and Zentalis Pharmaceuticals.

Supplementarydata.docx
Supplementary Figures and Tables

Download PDF

Version 1

posted

You are reading this latest preprint version

Spatial Heterogeneity of Infiltrating T cells in High-Grade Serous Ovarian Cancer revealed by Multi-omics analysis

Status:

Version 1

Abstract

Figures

Introduction

Results

Differential transcriptomic profiles across multiple sites in HGSOC

Distinct characteristics and differential composition of TILs across different lesions in HGSOC by scRNA-seq

Patient derived TMB is associated with skewed T cell differentiation

Tumor-specific but exhausted CD8 + cells preferentially infiltrate primary ovarian tumors, while non-tumor specific bystander cells are enriched in omentum metastases

Exhausted CD8 T cells enriched in primary ovarian tumors are clonally expanded

Exhausted CD8 T cells are a consequence of differentiation

CD4 Treg suppress the immune microenvironment in primary ovarian tumor sites

Inherent TME characteristics contribute to spatial differences of TIL status

Discussion

Materials And Methods

Clinical sample collection

Single cell collection

Single cell sorting and scRNA library construction sequencing by 10x genomics

TCR V(D)J sequencing

Bulk DNA and RNA isolation and sequencing

Bulk TCR sequencing

scRNA library construction sequencing by BGI DNBelab C4 platform

Cell clustering of BGI DNBelab C4 scRNA-seq dataset

Immunohistochemistry

Opal multiplex IHC

Gene expression quantification

Somatic mutation calling

Copy number calling and Tumor purity estimation

BRCA germline variants

Mutational process

Somatic SV detection

Chromothripsis

BFB detection

Neoantigen identification

Differentially expressed genes between Omen and Ovary

Bulk TCR-seq data analysis

TCR repertoire annotated by disease associated TCR database

Inferring neoantigen associated TCRs

Pairwise similarity calculation between TCRs

scRNA/TCR-seq data processing

Signature gene sets

Proliferation state definition

Identification of signature genes and TCGA data analysis

TCR analysis

Gene Set Enrichment Analysis

Ligand-receptor expression analysis

Trajectory analysis

Software versions

Statistical analyses

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1

Tumor-specific but exhausted CD8 ⁺ cells preferentially infiltrate primary ovarian tumors, while non-tumor specific bystander cells are enriched in omentum metastases