CXCL10 and CXCL9 differential expression profiles of 1,673 patients [mean months ± SD, 42.1 ± 34.2; 8% early stage (I, II), 92% late stage (III, IV)] in four independent cohorts were analyzed (Table 1). The patients in different cohorts were stratified according to CXCL10 and CXCL9 expression-related OS, and the cutoff value was determined by the survminer package. Overall, in four independent cohorts, we found that a lower CXCL10 expression was significantly associated with poorer prognosis (high risk) [standard cohort HR range: 0.69 (95% CI: 0.55–0.86; P = 0.003)]. Thus, we stratified the standard cohort and three test cohort patients into high-risk (lower CXCL10 expression) and low-risk groups (higher CXCL10 expression), respectively, by visually inspecting the curves of estimated survival (Figure 1A). Meanwhile, the lower CXCL9 expression was significantly associated with poorer prognosis in the standard (P = 0.002) and Test 3 groups (P < 0.001) (Figure 1B). The CXCL10 and CXCL9 expression levels showed a positive correlation in four independent groups of ovarian cancer (Figure 1C). Not surprisingly, combination group survival analysis revealed that patients with a high expression of CXCL10 and CXCL9 showed better prognosis (Figure 1D). Next, we verified the consistency of the three independent cohorts and the standard cohort by the multivariable Cox model (HR = high vs low risk) after adjusting the datasets, stages, and grades. The prognostic index of CXCL10-based high-risk classification in three test cohorts was obviously higher than that of The Cancer Genome Atlas (TCGA) standard cohort (Figure 1E). Although the results of subgroup analysis were heterogeneous, high-risk classification was positively associated with malignant neoplastic grades (G2 to G4) and pathological stages (II to III) in the merged group (standard plus three test cohort patients). We also observed that tumor invasion, residual status, and grades in the standard cohort were positively associated with the high-risk classification (Figure 1F).

In order to further characterize and understand the tumor immunobiology with CXCL10 expression, the construction scheme of the ovarian TME cell infiltration pattern was systematically evaluated. We evaluated the association between CXCL10 expression and immune cell populations from transcriptomic data (Figure 1G). Unsupervised hierarchical clustering of 568 samples with matched TME cell microarray probes in the standard cohort was presented. The heatmap depicted a strong correlation between high CXCL10 expression and T cells, cytotoxic cells, T helper 1 (Th1) cells, and different DCs. Likewise, tumor with low CXCL10 expression displayed a general lack of immune infiltration. Moreover, the same observation made across all test cohorts was also examined (Supplementary Figure 1).

Through the pROC package analysis (Robin et al., 2011), we observed that CXCL10 showed a predictive advantage in the infiltration group when compared to CXCL9 in four independent groups of ovarian cancer (Figure 2A, left). In addition, combining CXCL10 and CXCL9 improved the predictive value compared with that of CXCL10 or CXCL9 alone in merged samples (Figure 2A, right). Of note, the CXCL10/CXCL9 showed a robust prediction in the standard group compared with the other two groups (likelihood ratio test, P < 0.0001). Thus, CXCL10 presented a prior role in ovarian cancer when compared to CXCL9. To determine the optimal cluster immune cell types, the TME cell network depicted a comprehensive landscape of high-risk and low-risk classifications, including TME immune cell interactions and cell lineages of patients with ovarian cancer (Figure 2B). Four TME cell infiltration subtypes were shown in schematics (immune clusters A, B, C, and D). Moreover, cluster A was characterized by correlating with the infiltration of CXCL10-related immune cells in both standard and merged groups. However, the patients showed inconsistent survival with two intrinsic phenotypes of CXCL10/CXCL9-related risk and immune cell infiltration (log-rank P < 0.001), which might be attributed to the incongruous function of immune cells in TME (Figure 2C). Given that the key role of individual immune cells in TME can mediate tumor development and predict clinical outcome, we presented 23 types of TME immune cell infiltration as another survival marker in different cohorts of ovarian cancer (Figure 2D). As CXCL10-related risk was correlated with TME immune cell infiltration, we also explored the prognostic value of the merged group. Notably, patients with high-infiltration TME immune cells [DCs, Th1, T helper 2 (Th2), and T cells] and low-risk groups showed a synergistic function in improving OS (log-rank P < 0.01) (Figure 2E).

In addition to patient OS, TME immune cell infiltration tendency, intergenic relationship, and enriched biological processes were also associated with CXCL10 expression. In the standard cohort, as a whole, cytotoxic cells, INF γ-DC (aDC), DCs, Th1, Th2, and T-cell infiltration degrees were positively associated with CXCL10 expression (Figure 3A), and significantly positive correlations were revealed between CXCL10 expression and aDC, DCs, Th2, and T-cell infiltration values (Figure 3B). Furthermore, through overlapping of the immune cells in cluster A and the positive correlation group (Figure 3C), CXCL10 expression was mostly correlated with six of these cells’ infiltration, namely, aDC, DCs, Th1 cells, Th2 cells, T cells, and cytotoxic cells. When taking into account these selected immune cell probes, we utilized a PCA plot to show that samples were clustered primarily by transcription profile in different groups. A closer inspection revealed that the probe expression pattern possessed an obvious separation between high- and low-risk classifications in the standard cohort (Figure 3D). Despite individual variability, the graphics showed appreciable immune cell genetic dysregulation in different risk classifications of ovarian cancer.

Furthermore, we evaluated the enrichment of these selected immune probes in tumor biological function. Gene set variation analysis with two different groups (high risk and low risk) were analyzed by the GSVA package of R software (Hanzelmann et al., 2013), based on the KEGG and Gene Ontology (GO). A direct comparison of these immune probes in high-risk and low-risk classifications revealed that KRAS signaling upregulation (| log fold change| = 0.297, adjusted P-value < 0.0001) was the top enriched signature in high-risk ovarian cancer (Figure 3E). Meanwhile, low-risk-group immune cell signature genes were significantly associated with allograft rejection and IL2-STAT5 hallmark signaling (| log fold change| = 0.346, adjusted P-value < 0.0001). Furthermore, we used gene set enrichment analysis (GSEA) to perform all ovarian cancer transcription analyses from low risk to high risk and high infiltration to low infiltration (Figure 3F). We found that both the low-risk group with a higher CXCL10 expression and the high-infiltration group were closely connected to immune cell function [chemokine signal (ES > 0.62, P-value = 0.001), cytokine receptor interaction (ES > 0.5, P-value < 0.001), Th1/2 differentiation (ES > 0.69, P-value = 0.001), and cytosolic DNA-sensing pathway (ES > 0.7, P-value = 0.001)].

In terms of underlying mechanisms, we tried to investigate genetic instability on transcription. Based on genomic data in the standard cohort, we comprehensively analyzed the genetic alteration and single-nucleotide polymorphism (SNP) in ovarian cancer. We compared copy-number variations (CNVs) between high-risk and low-risk groups’ whole-genome sequencing data from TCGA tumor samples. In general, we observed that low-risk groups showed higher CNV in coverage (Figure 4A). Next, somatic alteration distributions of six immune cell signature genes from the high-risk group and low-risk group were analyzed (Figures 4B and Supplementary Figure 2B). A median of 162 alterations per sample (ranging from 20 to 1,922) in a total of 40,924 coding somatic alterations from 165 valid patient profiles was employed. We found that higher somatic tumor alteration load was obviously associated with high-risk classification patients. In addition, through the comparison between high- and low-risk classifications’ genetic alteration frequency, cytotoxic cell (RORA: 5.2 vs 0%; DUSP2: 3.4 vs 2%; CTSW: 2.6 vs 0%), DCs (PPFIBP2: 4.3 vs 0%), T cells (CD3G: 2.6 vs 0%; CD94: 2.6 vs 0%; CD6: 3.4 vs 0%), Th1 cells (IL12RB2: 6.0 vs 2.0%; ATP9A: 5.2 vs 3.9%; DGKI: 5.2 vs 0%; DPP4: 4.3 vs 2.0%; HBEGF: 2.6 vs 0%) and Th2 cell (CENPF: 5.2 vs 3.9%; CDC7: 3.4 vs 0%; CD3G: 2.6 vs 0%) infiltrations were most probably accompanied with altered genetic functions in high-risk classifications.

Moreover, the overall mutational pattern of these selected samples was dominated by C > T, and alteration patterns were similar between the high-risk and low-risk classifications Figure 4C). Meanwhile, the C > T and T > G average alteration frequency was higher in high-risk patients than in low-risk patients (Supplementary Figure 2C). We extracted five mutational signatures (i.e., signatures 1, 3, 5, 8, and 25) from the above data (Supplementary Figure 2A). The five signatures were annotated against the COSMIC signature version 2 (Alexandrov et al., 2013). Signature 1, featured by C > T transitions at CpG dinucleotides, is thought to be connected with age-related accumulation of spontaneous deamination of 5-methylcytosine. Notably, signature 1 was negatively associated with patient outcomes in human cancer, including ovarian cancer (Essers et al., 2020). The identification of this mutational signature might provide a new perspective to study the mechanism of chemokine-related TME formation and explanation of rising tumor incidence, as well as exploring individual mutations and their roles in cancer immunity and immunotherapy (Table 2).

To better clarify the functionality of the TME signature of CXCL10-mediated risk in ovarian cancer, we tested known TME signatures (Mariathasan et al., 2018) in four independent cohorts (Supplementary Figure 3A). These analyses confirmed that low-risk ovarian cancer patients were positively associated with immune-relevant signatures (antigen-processing machinery, immune checkpoint, and CD8⁺ T effector) (Figure 5A and Supplementary Figure 3A). Consistent with these results, low-risk ovarian cancer was notably linked to higher TME scores (Figure 5B). The Sankey diagram showed the comprehensive correlation of the different clinical parameters, TME characteristics, and ovarian cancer risk (Figure 5C). When examining the association between TME score and survival in four different cohorts, we found that a higher TME score was significantly associated with prolonged OS (Figure 5D and Supplementary Figure 3B). Moreover, in the standard cohort, TME scores remained statistically significant after taking into account CXCL10-related risk classification and pathological stages. As expected, consistent with the outcomes of CXCL10-related infiltrating patterns, the patients with higher immune cell infiltration and CXCL10 expression were positively associated with TME score (Figure 5E and Supplementary Figure 4A). The predictive value of TME score to tumor development was also confirmed in the standard (n = 568) and merged cohorts (n = 1,486) with valid pathological stage data, and we found that a decreased TME score was associated with tumor deterioration (Figure 5F and Supplementary Figures 4B,C).

We retrospectively collected the ovarian cancer gene expression from the Gene-Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases. In total, we gathered four cohorts with 1,673 samples from 10 public microarray datasets, including available follow-up time, survival status and international federation of gynecology and obstetrics (FIGO) stages. Patients without survival information were removed from further evaluation. For microarray data, we downloaded the raw data and clinical information from GEO repository¹ by the GEOquery package (Davis and Meltzer, 2007), then normalizing (Affy package) expression values of all probes and combined the datasets with the same microarray platform. As for the datasets in TCGA, Affymetrix Human Genome U133A Array microarray data of gene expression was downloaded from the Genomic Data Commons (GDC,²) by using the R package TCGAbiolinks, which was specifically developed for integrative analysis with GDC data (Colaprico et al., 2016). The non-biological technical biases batch effects between different datasets within the same platform were adjusted by ComBat algorithm (Johnson et al., 2007). The somatic mutation data was acquired from TCGA database by using the maftools package (Mayakonda et al., 2018). The TCGA ovarian cancer dataset from Affymetrix SNP 6.0 at Genome Analysis Platform of the Broad Institute was downloaded for individual Copy Number Variation (CNV) analysis. Data was analyzed with the R (version 3.6.1) and R Bioconductor packages.

We used Bindea’s immune cells genes signature (Bindea et al., 2013), which is well known as a highly sensitive and specific discrimination of 23 human immune cells phenotypes, including mast cells, DCs, aDC, Natural killer cells (NKs), Macrophages, and T cells subsets, etc. Moreover, we constructed TME signature gene sets from the study of Mariathasan (Mariathasan et al., 2018), and it includes 14 TME categories according to different molecular and signaling pathway functions, such as immune-relevant signatures, mismatch-relevant signatures and stromal-relevant signatures.

To quantify the relative abundance of immune cell infiltration in ovarian cancer TME, we employed single-sample gene-set enrichment analysis (ssGSEA) algorithm that defined the enrichment score of immune cells genes set in each patient within a given dataset. In order to identify the immune cells were associated with CXCL10-determined risk classification, we grouped the cells into four distinct groups, such as Immune cluster-A, Immune cluster-B, Immune cluster-C, and Immune cluster-D. Immune cells differentially enriched scores were determined by using the limma R package (Ritchie et al., 2015), which implemented an empirical Bayesian approach to estimate score changes by using moderated t-tests.

To estimate the risk classification related ovarian cancer TME signature, we constructed a set of scoring system between two risk classifications of patients with ovarian cancer. The construction of TME signature was performed as follows: Two different risk classifications and all samples gene were extracted. Then, a consensus clustering algorithm was used to define TME-relevant genes signature, and principal component analysis (PCA) was conducted. Principal component 1 was extracted to serve as the genes signature score. After that, we evaluated the TME score by using a method similar to GGI (Sotiriou et al., 2006; Zeng et al., 2019).

TME score = ΣPC1_i − ΣPC1_j

where i is the signature score of the low-risk classification, and j is the signature score of the high-risk classification.

To further understand the genes function between two risk classification patterns, we performed a Gene set variation analysis (GSVA) with six immune cells gene sets. GSVA is a non-parametric and unsupervised method, which is generally used to evaluate the variation in pathway (KEGG) and biological process activity (GO) (Hanzelmann et al., 2013). Furthermore, we employed ClusterProfiler R package (Yu et al., 2012) to perform the risk classification associated pathways and biological processes, basing on adjusted expression data of all transcripts. We identified hallmark pathways among low to high risk and high to low infiltration individuals by using gene set enrichment analysis (GSEA) (Subramanian et al., 2005). The gene sets examined in pathway and biological process were obtained from MSigDB database of Broad Institute. Adjusted P-value < 0.05 was considered as statistically significance.

The investigation of CNV was presented based on DNA profiling of individuals on Affymetrix SNP 6.0 Platform. We divided the segment file into two risk classifications and performed with the GISTIC 2.0 algorithm (Mermel et al., 2011). The somatic point, missense, insertion, and deletion mutations were annotated by using information from TCGA somatic mutation dataset. For the two classification mutation plots, we compared 45 immune-cell signature genes alteration levels of selected patients. Then, we employed the framework proposed by Kim et al. (2016) to extract mutational signature. This framework based on Bayesian variant non-negative matrix factorization and it can automatically determine the optimal number of mutational signatures. The mutation portrait matrix was factorized into non-negative matrices and corresponding mutational activities. The columns of non-negative matrix demonstrated the number of extracted alteration, and the rows indicated the 96 mutational contexts, which derived from a combination of six mutational types (i.e., C > A, C > G, C > T, T > A, T > C, and T > G) and their 5′ and 3′ adjacent bases (Dulak et al., 2013; Li et al., 2018). The corresponding mutational matrix revealed the individual mutation activities and the corresponding mutational signatures 1, 3, 6, 8, and 25. These mutational signatures were curated by the Catalog of Somatic Mutations in Cancer (COSMIC) Version 2 (COSMIC, ³).

In this study, statistical analysis was mainly performed by using R software and SPSS software (version 24.0). For the two groups comparison, statistical significance was evaluated by unpaired t tests. The Kruskal–Wallis tests were used to conduct difference comparisons of three or more groups. In the survival analysis, survival curves were generated via the Kaplan-Meier in each data set, and log-rank tests were utilized to determine significance of differences. The cut-off values for risk classification were evaluated based on the association between patient survival and CXCL10 expression in each cohort by using the survminer R package. The survminer can uninterruptedly screen all possible cut points to identify the maximum rank statistic. The associations between characteristics and OS were evaluated by Cox proportional hazard models. Correlation coefficients between TME infiltration score and expression of CXCL10 were computed by Spearman and distance correlation analyses. The waterfall plot of maftools package (Mayakonda et al., 2018) was applied to present the mutation landscape in patients with high and low risk ovarian cancer subtypes in standard cohort. All statistical P value were two-side, with P < 0.05 as statistically significance.