Fragments per kilobase million (FPKM) data of CC (N = 307) including gene expression profiles, somatic alteration data, and clinical data were downloaded from the Cancer Genome Atlas (TCGA) Genomic Data Commons (GDC) data portal (https://portal.gdc.cancer.gov/). The TCGA database is based on tissue specimens for high-throughput sequencing. The FPKM data were translated into transcripts per kilobase million (TPM). Duplicate recorded samples and overall survival (OS) time or survival status unavailable were excluded. Ultimately, 291 CC patients from the TCGA cohort were enrolled in the analyses, including 241 squamous cell carcinomas, 46 adenocarcinomas, and 4 adenosquamous carcinomas. Specific clinical information on CC patients was supplemented in Data Files S1 .

These 184 ICSs were based on the aggregation of databases such as ImmPort, CIBERSORT, and ImSig (10). Detailed information was listed in Data Files S2 . The normalized enrichment score (NES) generated by single sample gene set enrichment analysis (ssGSEA) with the R package “GSVA” (version 1.42.0) was considered as the infiltrate level of each ICS (11). In this study, only 183 ICSs were evaluated for follow-up analyses due to lack of some marker genes in the transcriptome atlas.

Meanwhile, the IMvigor210 cohort (http://research-pub.gene.com/IMvigor210CoreBiologies/) (N = 348) were also enrolled to validate our findings using the R package “IMvigor210CoreBiologies”, which are patients with metastatic urothelial cancer (mUC) receiving PD-L1 inhibitor (12). Raw count was also translated to TPM to represent gene expression in the IMvigor210 cohort. In addition, to expand our findings in gynecologic tumors, ovarian cancer (OC) and endometrial cancer (EC) data were downloaded and processed from the TCGA database. Finally, 374 OC patients and 539 EC patients were included. The flow of this study was shown in Figure 1 .

The univariate Cox regression was applied to obtain survival-related ICSs (P< 0.050). Subsequently, according to the infiltrate levels of survival-related ICSs, hierarchical agglomerative cluster of CC patients was performed using R package “ConsensusClusterPlus”. This algorithm was repeated 1,000 times to obtain stable classification. The differentially expressed genes (DEGs) between ICSs clusters were analyzed with false discovery rate (FDR)< 0.050 and absolute fold-change > 2, which was implemented by employing the R package “limma”. To identify the genomes and pathways enriched by DEGs, the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway functional enrichment analyses were performed using the R package “clusterProfiler” (version 4.2.2).

In the TCGA cohort, the univariate Cox regression was applied to obtain survival-related genes from DEGs (P< 0.050). Then, the TCGA-CC samples were randomly divided into training (n = 204) and validation (n = 87) sets at 7:3. Least absolute shrinkage and selection operator (LASSO) and multivariate Cox regression analyses were used to construct an ICSs prognostic model in the training set. The calculation formula of Risk score was as follows:

In the formula, β_i was defined as the coefficient of genes correlated with survival and ExpGene _i was the expression value of the corresponding gene in each sample. The cut-off value was determined by the “surv_cutpoint” function of the R package “survminer”, which calculates statistics based on maximally selected rank statistics. The principle of this function to determine the optimal cutoff value is to obtain the two groups with the most statistically significant difference in survival rates through multiple simulations.

Tumor tissues from a total of 23 CC patients were obtained from the Harbin Medical University Cancer Hospital. These samples were obtained from patients who underwent radical surgery, and all patients were untreated prior to radical resection. Specific clinical information is supplemented in Data Files S3 . All specimens were collected in accordance with the ethical standards of the Committee for Human Experimentation. The expressions of S1PR4 and FKBP10 on the CC tissue were performed by immunohistochemistry (IHC) staining. Tissue sections were incubated with a primary antibody against S1PR4 or FKBP10 at 4°C overnight and then incubated with horseradish peroxidase combined with goat anti-rabbit antibody (PV-6001, ZSGB) at room temperature for 30 mins. Tissue sections were stained using DAB and counterstained with hematoxylin. The results of the experiment were analyzed by two doctors and two pathologists. The rules are as follows: Immunoreactive score (IRS) = SI (staining intensity) × PP (percentage of positive cells). SI was assigned as: 0 = negative; 1 = weak; 2 = moderate; 3 = strong. PP was defined as 0 = 0%; 1 = 0–24.9%; 2 = 25–49.9%; 3 = 50–74.9%; 4 = 75–100%. All of the included patients were dichotomized into two groups based on the median score.

To further explore the immune infiltration landscape in two prognostic subgroups, multiple immune-related algorithms were utilized, such as the ESTIMATE, CIBERSORT, and MCPcounter algorithms (13–15). In addition, we also analyzed the distribution of immune subtypes in high and low-risk groups, which were categorized in previous studies (16, 17). To further reveal the underlying mechanisms of CC, gene set variation analysis (GSVA) was performed using the R package “GSVA” (adjusted P< 0.050). The gene set “c2.cp.kegg.v7.4.symbols” was downloaded from Molecular Signatures Database (MSigDB, http://www.broad.mit.edu/gsea/msigdb/).

The main biomarkers for the prediction of immune efficacy included immune mark genes, immune function characteristics, tumor mutational burden (TMB), histocompatibility complex (MHC) molecules, chemokines, cytolytic activity (CYT), and stimulator of interferon genes (STING) (18–20). Therefore, we further explored the differences in these immunotherapy-related biomarkers between high and low-risk groups.

We downloaded the immunophenoscore (IPS) from the cancer immunome atlas (TCIA) (https://tcia.at/home) to predict responses to ICB (21). IPS was calculated based on the expression of MHC molecules, immunomodulators, effector cells (ECs) and suppressor cells (SCs). It included four types of scores, ips_ctla4_pos_pd1_pos, ips_ctla4_pos_pd1_neg, ips_ctla4_neg_pd1_pos, and ips_ctla4_neg_pd1_neg, to better predict the efficacy of anti-CTLA-4 and anti-PD-1 antibodies. IMvigor210 cohort and two gynecologic tumor (TCGA-OC and TCGA-EC) cohorts were also used to validate the predictive value of Risk score for immunotherapy.

All statistical analyses were performed with R software (version 4.1.3). The Kaplan–Meier plotter was employed to depict survival curves. The Wilcoxon test was carried out to compare the difference between two groups, and the correlation coefficient was computed using the Spearman analyses. Two-tailed P< 0.050 was deemed statistical significance.

The TCGA cohort was divided into two clusters based on survival-related ICSs by using hierarchical agglomerative cluster ( Figure 2A ; Figures S1A, B ). The accuracy of the clustering was verified using principal component analysis (PCA) ( Figure 2B ). Furthermore, Kaplan-Meier survival curves showed significant survival difference between two ICSs clusters (P = 0.016; Figure 2C ). Differential correlation patterns of survival-related ICSs between two clusters were visualized as heatmap ( Figure S1C ). Enriched biological processes of ICSs clusters were summarized by Figures S1D and S1E , with specific data in the Data Files S4 .

To furtherly clarify the intrinsic substrates leading to different survival outcomes between two ICSs clusters, we performed a series of immune correlation analyses. CIBERSORT algorithm revealed that ICSs cluster A was characterized by high naive B cells, M0 macrophages, activated mast cells infiltration, which might be the cause of poor prognosis. By contrast, ICSs cluster B was marked by high CD8 T cells, memory activated CD4 T cells, follicular helper T cells, and M1 macrophages infiltration ( Figure 2D ). The heatmap of correlation coefficient was generated to visualize the cellular interaction of the tumor-infiltrating immune cell types ( Figure S1F ). MCPcounter and ESTIMATE algorithms similarly demonstrated a higher immune infiltration of ICSs cluster B ( Figure S1G ; Figure S1H ). In the subsequent analysis, we compared the differential expression levels of 27 common immune marker genes in two clusters (22, 23). The results exhibited that most of the immune marker genes were comprehensively elevated in the ICSs cluster B ( Figure 2E ). Furthermore, ssGSEA of specific gene sets demonstrated that the ICSs cluster B was highly active in multiple immune function pathways ( Figure 2F ), consistent with the results of the above immune correlation analyses. GSVA identified that the B cell receptor, T cell receptor, and JAK/STAT pathways were significantly activated in the ICSs cluster B ( Figure 2G ). The results were supplemented in the Data Files S5 .

The 986 DEGs (FDR< 0.050 and absolute fold-change > 2) were included in the univariate Cox regression model and 321 genes (P< 0.050) were found to be significant. Heatmap results displayed that the gene expression levels in ICSs cluster B were generally higher than that in ICSs cluster A ( Figure 3A ). All the 321 significant genes were incorporated into the LASSO and multivariate Cox regression model ( Figures 3B, C ). Finally, eight IRGs were contained in the ICSs prognostic model. The comprehensive Risk score was calculated as follows: Risk score = (0.13780 * CA9) + (0.29263 * FKBP10) + (-0.27821 * CKB) + (-0.40748 * GLIPR2) + (-0.42239 * ISG20) + (-0.42787 * S1PR4) + (-0.26784 * SDS) + (-0.20766 * VTCN1). Patients with CC were divided into high and low-risk groups using the cutoff value (1.27508) as the dividing line. The above results were showed in the Data Files S6 .

In the total population, training and validation sets, the distribution curves and survival scatter plots indicate that patients with high-risk scores have a poorer prognosis ( Figures S2A–F ), and Kaplan-Meier survival curves demonstrated that patients showed a significant difference between high and low-risk groups in survival rate (P< 0.001; Figures 3D–F ). In the total population, time-dependent receiver operating characteristic (ROC) curves showed that the ICSs prognostic model had a strong prognostic accuracy with the area under the curve (AUC) of 0.870 in 1 year, 0.785 in 3 years and 0.774 in 5 years ( Figure 3G ). The results of time-dependent ROC curves for the training and validation sets were shown in Figures 3H, I .

Univariate and multivariate Cox regression showed that Risk score and Stage were independent prognostic factors for CC ( Data Files S7 ). By combining the independent prognostic factors, we constructed a nomogram that serves as a clinically relevant quantitative method by which clinicians could predict the mortality of CC patients ( Figure 4A ). In addition, calibration plots indicated that the performance of the nomogram was similar to that of the ideal model ( Figure 4B ). The decision curve analysis (DCA) also revealed that the nomogram had high potential clinical utility ( Figure 4C ). We also compared the AUC of the Risk score and Stage, and the results showed that the Risk score had better predictive power ( Figure 4D ). Finally, we compared the distribution of Stage between high and low-risk groups, as shown in Figure 4E , with a higher proportion of Stage I-II patients in the low-risk group.

To verify our previous research results, IHC staining was performed on tissue samples collected from patients with CC to explore the expression of oncogenes (FKBP10) and tumor suppressor genes (S1PR4) with the largest weight coefficient in the ICSs prognostic model. Figures 5A and B show immunohistochemical images of FKBP10 and S1PR4 in different differentiation states. The immunohistochemical results showed that the expression level of oncogene FKBP10 was negatively correlated with the degree of differentiation, while the expression of tumor suppressor gene S1PR4 was the opposite. In addition, this study used immunohistochemical images of tissues at different stages to demonstrate that the expression level of FKBP10 gradually increased with increasing tumor stage, and that of S1PR4 gradually decreased with increasing tumor stage ( Figures 5C, D ). This evidence confirmed the expression of two key genes in CC tissues, and our results had a higher degree of confidence.

To further explore the correlation between the prognosis and TIME, we analyzed the immune infiltration landscape of CC. Alluvial diagram illustrated that most patients in ICSs cluster B have low Risk score and more alive status ( Figure 6A ; Figure S3A ). In addition, we assessed the differential correlation pattern of prognosis-related ICSs in the high and low-risk groups, and the result was shown in Figure 6B . Further analyses showed that low-risk group had higher immune score, stromal score and estimate score ( Figure 6C ). CIBERSORT and MCPcounter algorithm displayed that the immune cell infiltrating types of the low-risk group were similar to ICSs cluster B, and the high-risk group was consistent with ICSs cluster A ( Figures 6D, E ). The association heatmap visualized a negative correlation between Risk score and multiple tumor-infiltrating immune cells ( Figure S3B ). The distribution of immune subtypes showed that the low-risk group was mainly distributed in IS4 and C2 types, which reflected that the low-risk group had a more favorable anti-tumor immune response ( Figures 6F, G ). GSVA was used to further explore the potential role of Risk score in biological processes, which identified that B cell receptor, T cell receptor and ATP-binding cassette (ABC) transporters pathways were significantly activated in the low-risk group ( Figure 6H ). Detailed results of GSVA were listed in the Data Files S8 .

To further elucidate the effects of Risk score in the context of immunotherapy, we explored the associations between Risk score and several well-known immune marker genes. It was shown that immune marker genes were expressed at higher levels in the low-risk group ( Figure 7A ). We analyzed 29 immune function-related characteristics and found that the low-risk group had a more favorable immune activation phenotype, suggesting that they may have a more intense immune response ( Figure 7B ). Considering the great clinical significance of TMB for immunotherapy, we sought to explore the intrinsic correlation between TMB and Risk score. The “maftools” R package was carried out for assessing the distribution of somatic mutation in the high and low-risk groups, and depicted the top 20 driver genes with the highest alternative frequencies in Figures 7C, D . The heatmap of correlation coefficient demonstrated the interrelationship of the top 20 mutated genes in CC patients ( Figure S3C ). We found no difference in the prognosis of patients in the high and low TMB group ( Figure S3D ), as well as no significant correlation between TMB and Risk score ( Figures S3E , F ). In addition, we found that the expression levels of MHC molecules, and chemokines, which were responsible for the movement of immune cells, were comprehensively elevated in the low-risk group ( Figures 7E, F ). CYT and STING were relatively higher in the low-risk group (P< 0.001; Figures 7G, H ), and Risk score was significantly negatively associated with relatively higher ( Figures 7I, J ). Taken together, these results suggested that Risk score based on eight IRGs may be potential predictors of CC immunotherapy efficacy.

To assess the ability of the Risk score as a biomarker for predicting clinical response to ICB treatment, we assessed the immunogenicity of two prognostic subgroups by IPS analyses. The low-risk group had higher ips_ctla4_pos_pd1_pos, ips_ctla4_pos_pd1_neg, ips_ctla4_neg_pd1_pos, and ips_ctla4_neg_pd1_neg scores in the TCGA-CC cohort ( Figures 8A–D ). In the TCGA-CC and IMvigor210 cohort, Risk score was negatively correlated with MHC scores, EC scores and IPS scores. Regarding the SC score, opposite results were obtained ( Figures 8E, F ). In the subsequent analyses, IMvigor210 cohort were assigned high and low-risk score. As shown in Figures 8G, H , while there was no statistical difference in Risk scores between patients with complete response (CR)/partial response (PR) and those with stable disease (SD)/progressive disease (PD), clinical response rates with anti-PD-L1 therapy were higher in patients with low-risk group than those with high-risk group (P = 0.020). Noteworthy, patients with low-risk score had a significantly better prognosis than those with high-risk score in the IMvigor210 and two gynecologic tumors ( Figures 8I–K ). In addition, IPS analyses indicated consistent results for the TCGA-EC cohort with the TCGA-CC cohort ( Figures S4A–D ), while only the ips_ctla4_pos_pd1_neg score was higher in the low-risk group of the TCGA-OC cohort ( Figures S4E–H ). These results indicated that patients in the low-risk group may have a better response to immunotherapy.