The HRD Score was significantly correlated with the prognosis and molecular characteristics of TCGA-OSC cohort.
According to the HRD-algorithm, LOH, TAI and LST were used as the basis for calculating the HRD score. Optimal cutoff scores were determined by assessing the score which had the minimum p value of the log-rank test (Fig. 2A). HRD status was defined as HRD-positive if HRD score was༞57; HRD status was defined as HRD-negative if HRD score was ≤ 57. The Kaplan-Meier survival curve (Fig. 2A) showed that overall survival of patients in the HRD-positive group is much longer than the cases in the HRD-negative group (hazard ratio (HR) = 0.49 (0.37,0.65), log-rank test, P < 0.00001). Subsequently, we investigated the correlation between the HRD score and other hallmarks of genomic instability, including somatic mutation counts, fraction genome altered and MSI. The somatic mutation pattern was significantly different between the HRD-positive group and the HRD-negative group. The median value of somatic cumulative mutations in the HRD-positive group was significantly higher than that in the HRD-negative group (Wilcoxon signed-rank test, P < 0.0001; Fig. 2B). We next compared the fraction genome altered scores between the HRD-positive group and the HRD-negative group. As shown in Fig. 2C, the fraction genome altered scores in the HRD-positive group were higher than those in the HRD-negative group (Wilcoxon signed-rank test, P < 0.05; Fig. 2C). The plot of “somatic mutation counts vs fraction genome altered” clearly showed that the distribution of points in the HRD-positive group was concentrated in the upper right of the coordinate system, while that of the HRD-negative group was scattered (Fig. 2D). The results of transcriptome level analysis were consistent with those at the genome level: through GSEA analysis, it was found that the three signal pathways with the most significant differences between the two groups were DNA replication, homologous recombination and mismatch repair (Supplementary Fig. 1A). There was no difference in the microsatellite instability (MSI) status of the two groups (Supplementary Fig. 1B).
The genomic landscape of the HRD-negative group showed a high proportion of actionable mutations in NF1 and CDK12
Genomic characteristics, such as the oncogene activation (e.g., ERBB2 amplification, EGFR tyrosine kinase mutation) and inactivation of tumor suppressor genes (e.g., MMR, BRCA1/2) have shown a strong correlation with clinical response to target therapy. Therefore, we compared the genomic mutational landscape between the HRD-positive group and the HRD-negative group. The results showed that the genomic landscape of the HRD-negative group was significantly different from that of the HRD-positive group. Only nine of the top twenty genes with the highest mutation rate in the two groups overlapped (Supplementary Fig. 1C). The mutational landscapes of these two subgroups displayed a distinct mutation ratio in TP53 [94.0% (HRD-negative) vs 62% (HRD-positive)], and the mutation classification of the HRD-negative group was more abundant, including a higher proportion of frame shift del, nonsense mutation and so on (Fig. 3A, B). Through the screening of actionable genes in the OncoKB database (https://www.oncokb.org/actionableGenes), among the 20 genes with the highest mutation frequency in the HRD-negative group, two genes were biomarkers for targeted drugs (NF1 and CDK12). Moreover, most of the variant classifications of these two genes were those affecting gene structure. Patients with mutations in these two genes accounted for 13% of the HRD-negative group. The mutational landscapes of HR genes in the two subgroups also exhibited a distinct difference. HR gene mutations in the HRD-positive group were mainly concentrated in BRCA1/2 (7%), while the HRD-negative group was scattered (Supplementary Fig. 2A, B). Besides, we also compared the copy number variation of HR genes in the two subgroups, and observed BRCA2 homozygous deletion was only present in the HRD-positive group (Supplementary Fig. 2C).
CXCL11 expression associated with the HRD score and its prognostic value in OSC
To identify RNAs associated with the HRD score, the TCGA-OSC cohort was sorted in ascending order of HRD scores, and the last 20% (n = 70) and the top 20%(n = 70) of the patients were selected to identify DEGs. Utilizing the egdeR method, a total of 124 DEGs were screened out. Among them, 38 RNAs were found to be upregulated and 86 to be downregulated in the HRD-positive group (Supplementary Table 2). Then, 124 differentially expressed RNAs were used to perform unsupervised cluster analysis on 348 TCGA-OSC samples. As shown in Fig. 4A, we found that not all DEGs clustered the HRD-positive and the HRD-negative groups well in the entire TCGA-OSC cohort. Only the DEGs in the red block region were able to cluster the HRD-positive group and the HRD-negative group well. To further screen out DEGs related to the HRD score and prognosis of the patients, the univariate analysis was conducted in the 124 DEGs for the whole TCGA-OSC cohort. A total of 17 genes with prognostic potentiality were identified by the univariate analysis and log-rank test (P < 0.05) (Supplementary Table 3). The 17 HRD-related genes were then subjected to Lasso-Cox proportional hazards regression and tenfold cross-validation to identify the best gene model. The Lasso coefficient profile plot was produced against the log(lambda) sequence, and the optimal lambda method resulted in only one optimal coefficient (Fig. 4B, C). The only one HRD-related gene was CXCL11. A heatmap of CXCL11 expression and the HRD score and the scatterplot of overall survival (OS) with corresponding risk scores were illustrated in Fig. 4D. Kaplan–Meier analysis displayed that the survival outcomes of TCGA-OSC patients with high CXCL11 expression (CXCL11-positive) were significantly better than patients with low CXCL11 expression (CXCL11-negative) (hazard ratio (HR) = 0.39 (0.29,0.51), log-rank test P < 0.00001) (Fig. 4E). To verify whether the CXCL11 expression signature has similar prognostic value in different OSC cohorts, we further confirmed this phenomenon in two independent OSC cohorts in the GEO database (including GSE140082 and GSE30161): results from Kaplan–Meier analysis also showed that patients in the CXCL11-positive group demonstrated a better prognosis than those in the CXCL11-negative group (Supplementary Fig. 3A, B).
Comparison of immune cells infiltration within the CXCL11-positive group and the CXCL11-negative group
Gene expression of CXCL11 is essential for activating immune cells[25], suggesting that tumor infiltrating immune cells might be different in the CXCL11-positive group and the CXCL11-negative group. To validate this assumption, the TIMER algorithm[26] was applied to estimate enrichment of various immune cell types within different subgroups. We developed a heatmap with TIMER results to visualize the relative abundance of 6 immune infiltrating cell subpopulations from the TCGA-OSC cohort (Fig. 5A). As depicted in the heatmap, there were significant differences in immune cell infiltration between the two subgroups. Anti-tumor lymphocyte cell subpopulations, such as CD4+/CD8+ T cells and dendritic cells were enriched in the CXCL11-positive group (Wilcoxon signed-rank test, P < 0.01). The neutrophils were also enriched in the CXCL11-positive group (Wilcoxon signed-rank test, P < 0.001) (Fig. 5B). We then investigate the correlation of immune cell infiltration with the expression of CXCL11 by spearman correlation coefficients. The results revealed that the expression of CXCL11 was significantly associated with immune cell infiltration in the TCGA-OSC cohort (Fig. 5C). We also further analyzed the correlation between the immune cell infiltration signal and the expression of CXCL11 in the TCGA pan-cancer cohorts, and found similar results (Supplementary Fig. 3C).
Furthermore, GSEA on the gene expression profile of the CXCL11-positive group against the CXCL11-negative group revealed the CXCL11 expression signature related biological signaling pathway. Genes involved in antigen processing and presentation, autoimmune thyroid and cytokine receptor interaction signaling pathways were the most significantly enriched in the CXCL11-positive group (Fig. 5D). However, taste transduction, basal cell carcinoma and hedgehog signaling pathways were enriched in the CXCL11-negative group (Fig. 5E).
CXCL11 expression associated with molecules in antigen processing and presentation pathway
The results from the TIMER and GSEA analysis showed that there were significant differences between the CXCL11-positive group and the CXCL11-negative group in antigen processing and presentation pathway, hinting that the expression of antigen-related genes might be associated with CXCL11 expression. To prove this assumption, we explored the correlation of antigen-related genes with the CXCL11 expression by using the Pearson correlation coefficient. We found that the expression of MHC class I/II (I: HLA-A, HLA-B, and HLA-C; II: HLA-DP, HLA-DM, HLA-DOA, HLA-DOB, HLA-DQ, and HLA-DR) and antigen binding (B2M, TAP1/2 and so on) molecules were highly correlated with the CXCL11 expression signature (Fig. 6A). There were significant differences in the expression of HLA-A, HLA-B and other key antigen presenting molecules between the two subgroups (Wilcoxon signed-rank test, P < 0.001, Fig. 6B). The results were confirmed in the GEO validation cohort (Supplementary Fig. 4). Since antigen processing and presentation pathway plays a crucial role in immune recognition of predict (neo-)antigen produced by cancer cells, we further investigated the relationship between neoantigen load and the CXCL11 expression signature by Pearson correlation coefficient. Correspondingly, predicted neoantigen load was highly correlated with CXCL11 expression (Fig. 6C). These results revealed the correlation between CXCL11 expression and molecules in antigen processing and presentation pathway, suggesting that neoantigen vaccine might be potential therapy for CXCL11-positive ovarian cancer treatment.
CXCL11 expression associated with ICB-related genes
In recent years, ICB therapy, represented by anti-PD-1/L1, has played an increasingly important role in anti-tumor treatment. The characteristics of TIME and immune checkpoint genes in tumor cells have a profound impact on ICB therapy. Therefore, we collected more than 40 common ICB-related gene signatures and analyzed the relationship between CXCL11 expression and ICB-related genes. As displayed by heatmap, CXCL11 expression was significantly correlated with the expression of multiple ICB-related genes (Fig. 7A). Ten of the most relevant ICB-related genes were: LAG3, ICOS, CTLA4, CD48, HAVCR2, PDCD1(PD-1), PDCDILG2(PD-L2), TIGIT, CD274(PD-L1) and CD86, and their expression levels were enriched in the CXCL11-positive group (Wilcoxon signed-rank test, P < 0.001, Fig. 7B). Generally, the key regulatory factors involved in immunity perform similar functions in different tissues. We thus explored CXCL11 expression and ICB-related gene signatures across cancer types. We found that the co-expression of CXCL11 and ICB-related genes was not only present in ovarian cancer, but also in 32 other cancer types (Fig. 7C). Taken together, these results show that CXCL11 expression may be a potential predictor of ICB therapy.
Olaparib-treated ovarian cancer cells up-regulate CXCL11 expression in vivo and in vitro
It has been reported that PARPi triggers robust local and systemic antitumor immunity involving both adaptive and innate immune responses through a STING-dependent antitumor immune response in mice bearing ovarian tumors[27, 28]. To explore the involvement of CXCL11 expression responses to PARPi, we re-analyzed the RNA seq data of high-grade serous ovarian cancer tumor tissues harvested from tumor-bearing mice after 18 days of treatment with olaparib or vehicle (GSE120500). Strikingly, boxplots showed markedly upregulated expression of CXCL11 in tumors treated with olaparib compared with vehicle control (Fig. 8A). We next compared the expression levels of ICB-related genes with the highest correlation with CXCL11 between the olaparib treatment group and the control group. As shown in Fig. 8B, the expression of these genes in the olaparib treatment group was significantly higher than that in the vehicle control group (Wilcoxon signed-rank test, P < 0.01). To further validate that olaparib-treated ovarian cancer cells activate the CXCL11 expression signature, we conducted in vitro experiments. As measured by RT-qPCR (Fig. 8C and D), olaparib treatment caused significant upregulation of CXCL11 mRNA expression in multiple ovarian cancer cell lines. Together, these data indicate that olaparib could upregulate the expression of CXCL11 and its related immune checkpoint genes in ovarian cancer cells in vivo and in vitro.