Patients with treatment-naive cervical squamous cell carcinoma, who were admitted to the Affiliated Cancer Hospital of Xinjiang Medical University between November 2021 and August 2022, were included in the study. Relevant clinical data were systematically collected for analysis.

Inclusion criteria for patients were (1): histopathologically confirmed cervical squamous cell carcinoma with no prior treatment for cervical cancer at the time of enrollment; (2) Karnofsky Performance Status (KPS) score ≥80; (3) age ≥ 18 years; (4) complete serological and clinical data. Exclusion criteria were: (1) pathologically confirmed types other than squamous cell carcinoma; (2) history of any other malignancy; (3) concurrent chronic infections, infectious diseases, immune system disorders, etc.; (4) pregnant females. We enrolled a total of fifty eligible patients, with ages ranging from 30 to 82 years and a median age of 55 years. Among them, 15 were elderly individuals aged 60 and above, and 35 were 60 years or younger. Histological differentiation revealed 17 cases of low differentiation and 33 cases of moderate differentiation. Clinical staging included 9 cases of stage I, 23 cases of stage II, 16 cases of stage III, and 2 cases of stage IV. HPV16-positive cases were observed in 40 instances, with 4 other HPV-positive patients and 6 HPV-negative cases. Lymph node metastasis was present in 14 cases out of 50. Surgical interventions were performed in 31 cases, while 19 cases did not undergo surgery while 19 cases did not undergo surgery at the initial treatment. Among the patients, 17 were premenopausal, and 33 were postmenopausal. In terms of tumor infiltration depth, 7 cases had superficial infiltration (1/3), 8 had intermediate infiltration (1/3), and 16 had deep infiltration (1/3). Additionally, 20 cases presented with intravascular cancer emboli, while 13 cases had no vascular involvement. In 36 cases, the Squamous cell carcinoma antigen(SCC) was ≥1.5, and in 14 cases, it was <1.5. All patients were staged using the 2018 Federation International of Gynecology and Obstetrics(FIGO) classification. This research obtained approval from the Ethics Committee of the Affiliated Cancer Hospital of Xinjiang Medical University, and all participants provided informed consent before their inclusion.

Ten milliliters of peripheral blood obtained from patients was anticoagulated with ethylene diamine tetraacetic acid (EDTA), followed by isolation of peripheral blood mononuclear cells (PBMCs) using density gradient centrifugation. The total cell count was determined using Trypan Blue staining. Fresh tissues (0.5 cm * 0.5 cm), acquired from treatment-naive cervical cancer patients. The samples were then dissociated using C tubes (miltenyi-biotec, catalog NO. 130-093-237) and the Human Tumor Dissociation Kit (miltenyi biotec, catalog No. 130-095-929) as per manufacturer instructions. The dissociated cells were then passed through a mesh with a pore size of 50 μm to create a solution consisting of single cells.

The isolated peripheral blood and tissue mononuclear lymphocytes were counted, and 6 × 10^5 cells were taken. After incubating with Live/Dead BV510 dye for 20 minutes, surface staining with antibodies was performed for an additional 20 minutes. Subsequently, 200 μL of fluorescence-activated cell sorting (FACS) fixing buffer was added to fix the cells, and flow cytometry was employed for detection. The BD LSRFortessa flow cytometer was purchased from BD in the United States, Live/Dead BV510 dye was obtained from Life Technologies in the United States, and materials including Fetal Bovine Serum (FBS), RPMI 1640 cell culture medium, human lymphocyte separation solution, Human Tumor Dissociation Kit (Miltenyi), and Phosphate Buffered Saline (PBS) were all obtained from Sigma in the United States. Monoclonal antibodies, specifically anti-human CD3-AF700, anti-human CD8-APC-H7, anti-human CD4-FITC, anti-human 4-1BB-BV786, anti-human ICOS-BV711, and anti-human CD28-PerCP-CyTM5.5, were purchased from BD in the United States.

Single-cell transcriptome files of CC samples were obtained from the GEO database for datasets GSE168652 and GSE171894. Gene expression and matched survival data for 304 CC patients were downloaded using the TCGAbiolinks package, followed by normalization of the gene expression data. The 304 patients were divided into two groups according to median PD-1 and 4-1BB expression levels. This yielded a final cohort of 250 CC patients.

Using the ‘pRRophetic’ Package (25), pharmacological sensitivity analysis was performed. Ridge regression determined the Half-Maximal Inhibitory Concentrations (IC50) using data from the Genomics of Drug Sensitivity in Cancer (GDSC) Database (https://www.ancerrxgene.org/). The Cancer Immunome Atlas (TCIA) database (https://tcia.at/) (26) is where the Immunophenoscores (IPS) of CC were obtained.

The metastatic urothelial tumors cohort (27) (IMvigor210) and the advanced melanoma cohort (28) (BMS038), both treated with immunotherapy, were included in the analysis. The R package “IMvigor210CoreBiologies” preprocessed the IMvigor210 data. RNA-seq data underwent filtering and normalization using the “edgeR” (29) package, followed by transformation using the “voom” function within the “limma” (30) package in R. Quantification of PD-1 + 4-1BB was conducted in the IMvigor210 and BMS038 cohorts.

Statistical analysis was conducted using SPSS 26.0 software. The chi-square test or Fisher’s exact probability test compared count data, while the t-test was applied for quantitative data comparisons. The rank sum test was used for non-normally distributed data. Overall survival (OS), disease specific survival(DSS) and progression free interval(PFI) were estimated via the Kaplan-Meier method, with the Log-rank test, and Cox proportional risk regression model for multi-factor analysis. R (version 3.6.3) was used for statistical analysis and visualization, with ggplot2 (version 3.3.3) utilized for visualization. Differences with p-values below 0.05 were considered statistically significant.

To investigate the co-expression levels of PD-1/4-1BB, PD-1/ICOS, and PD-1/CD28 in CC tissues and blood, a paired analysis was conducted. Flow cytometry was utilized to detect the co-expression levels of PD-1/4-1BB, PD-1/ICOS, and PD-1/CD28 on CD8+ T cells from peripheral blood and matched cancer tissues of 50 cervical squamous cell carcinoma patients ( Figure 1A ). The co-expression levels of PD-1/4-1BB on CD8+ TILs were 2.59% (interquartile range, 1.05%-6.15%), PD-1/ICOS was 19.50% (interquartile range, 11.62%-34.37%), and PD-1/CD28 was 7.45% (interquartile range, 2.94%-15.60%). On peripheral blood CD8+ T cells, the co-expression levels of PD-1/4-1BB were 0.14% (interquartile range, 0.06%-0.27%), PD-1/ICOS was 0.58% (interquartile range, 0.34%-1.05%), and PD-1/CD28 was 2.93% (interquartile range, 1.75%-3.83%). Comparative analysis revealed that the co-expression level of PD-1/4-1BB on CD8+ TILs was higher than that in PBMCs ( Figure 1B ), with a median difference of 2.37% (interquartile range, 0.76%-6.07%), (P < 0.001). Similarly, the co-expression level of PD-1/ICOS on CD8+ TILs was higher than that in PBMCs ( Figure 1C ), with a median difference of 18.76% (interquartile range, 11.32%-34.11%) (P < 0.001). Moreover, the co-expression level of PD-1/CD28 on CD8+ TILs was also higher than that in PBMCs ( Figure 1D ), with a median difference of 4.52% (interquartile range, 0.10%-12.31%) (P < 0.001). In CD8+ PBMCs, the co-expression level of PD-1/4-1BB was lower than that of PD-1/ICOS and PD-1/CD28 ( Figure 1E ). On CD8+TILs, the co-expression level of PD-1/4-1BB was lower than that of PD-1/ICOS ( Figure 1F ), and the differences were statistically significant. Furthermore, 4-1BB expression was lower than that of ICOS and CD28 in CD8+PBMCs and CD8+TILs ( Supplementary Figure 2 ). Compared to PD-1- CD8+ TILs, 4-1BB exhibited a significant upregulation in PD-1+CD8+ TILs ( Figure 1G ).

Based on median values, the co-expression of PD-1/4-1BB, PD-1/ICOS, and PD-1/CD28 were categorized into high and low expression groups. Clinical characteristic analysis indicated that the co-expression level of PD-1/4-1BB on CD8+ TILs correlated with lymph node spread (P=0.019), FIGO stage (P=0.041), surgical intervention status (P=0.016), and SCC levels (P=0.039). PD-1/CD28 co-expression on CD8+ TILs correlated with SCC levels (P=0.03), whereas PD-1/ICOS co-expression showed no association with clinical features ( Table 1 ). Compared to PD-1/ICOS and PD-1/CD28, PD-1/4-1BB demonstrated a closer association with cervical cancer. Upon stratification by median values, the clinical feature analysis showed that the expression level of PD-1 in CD8+ TILs was correlated with lymph node spread (P=0.03) and tumor infiltration depth (P=0.002), while 4-1BB expression was related to surgical status (P=0.003) and FIGO stage (P=0.041) ( Supplementary Table 1 ). Therefore, our focus was directed towards the co-expression of PD-1/4-1BB.

To comprehensively explore the cellular composition and structure of cervical cancer, we obtained cervical cancer sample from the GSE168652 dataset in the GEO database for analysis. We constructed an atlas consisting of 15875 (tumor, 6517, normal, 9358) single cells that passed strict quality filtering. We first applied U-MAP dimensionality reduction and clustering to generate a 2D map containing 13 clusters. Using singleR annotation, we classified all cells into five distinct cell types: endothelial cells, fibroblasts, CD8+ T lymphocytes, macrophages, and epithelial cells ( Figure 2A ).

We employed U-MAP analysis to re-cluster CD8+ T cells into six clusters ( Figure 2B ). To identify the major cell subtypes, we annotated each cluster based on its marker gene expression (31). CD8+ T lymphocytes were segmented into six distinct cell clusters, including CD8_1 Exhaustion/cell-cycle (TK1), CD8_2 Exhaustion/heat shock protein (HSPA1A), CD8_3 Exhaustion (GIMAP6), CD8_4 Memory (CD55), CD8_5 Early activated cells (PFKFB3) and CD8_6 Memory/effector (LTB) ( Figure 2C ). Violin plots showed elevated expression of PD-1 in the CD8_1 cluster and 4-1BB in the CD8_6 cluster ( Figure 2D ).

We utilized CytoTRACE to assess the differentiation potential of distinct cell subtypes. The results showed that CD8_1, CD8_6, and CD8_2 were less differentiated, while CD8_4 and CD8_5 were highly differentiated ( Figure 3A ). In order to explore the evolutionary dynamics of the CD8+ T-cell lineage in cervical cancer, we performed pseudo-temporal cell trajectory analyses of six CD8+ T-cell subpopulations. Cells of the CD8_1 population were primarily situated at the final stage of the cell trajectory, and cells of the CD8_2 population were predominantly distributed in the initial stage, with four branching points in the differentiation process ( Figures 3B, C ). Variations in PD-1 and 4-1BB expression during cell differentiation showed low expression at the initiation stage and higher expression at the terminal stage ( Figure 3D ). PD-1 expression on CD8+ T cells positively correlates with 4-1BB expression level ( Supplementary Figure 3A ). The heatmap illustrated the gene expression changes from the initial state to cell fate 1 or 2, identifying six distinct transformation patterns( Figure 3E ), in which PD1 was also included. Enrichment analysis was performed on DEGs in various cell trajectory fates. GO analysis revealed that DEGs are mainly associated with positive regulation of response to external stimulus, positive regulation of leukocyte chemotaxis, and cytolytic granule. KEGG analysis indicated that DEGs are primarily related to Antigen processing and presentation ( Supplementary Figure 3B ). The dynamic expression patterns of PD-1 and 4-1BB are similar across different cell differentiation fates ( Figure 3F ).

To identify potential interactions between the six CD8+ T cells, we performed Cell Chat analysis for cell-cell signaling connections (19). Cell communication analysis revealed close interaction among the CD8_1 and CD8_6 cluster ( Figure 4A ). The CD8_1 cluster was mainly a signal sender and the CD8_6 cluster was mainly a signal receiver ( Figures 4B, D ). We examined particular ligand-receptor interactions between distinct cell subpopulations. The ligand-receptor pairs MIF-(CD74+CD44) and CD70-CD27 were significantly upregulated in the immune cell subsets CD8_1 and CD8_6 ( Figure 4C ). The pathway of macrophage migration inhibitory factor (MIF) is activated ( Figure 4E ). SCENIC analysis revealed specific activation of transcription factors NR3C1, E2F1, HOXB2 in the CD8_1 subset, and THAP1, FOXP3 in the CD8_6 subset ( Figure 4F ). The CD8_1 and CD8_6 subsets were defined as the PD-1/4-1BB high co-expression group, while the 2/4/5 subsets were defined as the PD-1/4-1BB low co-expression group. Marker genes from these two groups were selected and applied to the cervical cancer TCGA expression matrix. The expression ratio of the mean values of marker genes in the two groups was utilized to stratify them into high and low categories, followed by survival analysis. The outcomes demonstrate that the PD-1/4-1BB high co-expression group exhibits a more favorable prognosis ( Figure 4G ). Validation in the GSE171894 dataset consistently demonstrates a superior prognosis in the PD-1/4-1BB high co-expression group ( Figure 4H ).

To explore the immune-related characteristics of PD-1/4-1BB high and low co-expression groups, we conducted an analysis assessing their correlation with immune-infiltrating cells, immune scores, immune checkpoints, cytotoxic reactions, T cell functions, and TMB. Initially, the comparative analysis of immune-related cytotoxic reactions revealed that the high co-expression group of PD-1/4-1BB exhibited stronger cytotoxic responses, including elevated levels of Interferon Gamma (IFNG), Granzyme A (GZMA), Perforin-1 (PRF1), and Granzyme B (GZMB) ( Figure 5A ). The analysis of TMB demonstrated a higher TMB among those with high expression ( Figure 5B ). Later, compared the immune scores among the PD-1/4-1BB high and low co-expression groups, revealing that the Stroma score, Immune Score, and ESTIMATE Score were all higher in the high expression group of PD-1/4-1BB ( Figure 5C ). Analysis demonstrated that the high co-expression group of PD-1/4-1BB exhibited increased expression levels of CD8 T cells, CD4 memory activated, T cells follicular helper, Dendritic cells resting, and M1 macrophages compared with the low co-expression group ( Supplementary Figure 3C ), with significant differences. The ssGSEA score heatmaps revealed greater immune infiltration in the PD-1/4-1BB high co-expression group compared to the low co-expression group ( Figure 5D ). Immune checkpoint analysis indicated that most inhibitory and stimulatory immune checkpoints had higher expression in the high co-expression PD-1/4-1BB group ( Figure 5E ). These imply a stronger association among the high co-expression group of PD-1/4-1BB and T-cell functionality. Afterwards, we conducted additional research on the distribution patterns of the top 10 somatic mutations utilizing the TCGA dataset with the aid of the “maftools” package. The most prevalent mutations were identified in TTN and PIK3A ( Figures 5F, G ).

Based on the TCGA database, Kaplan-Meier analysis revealed that patients who exhibited an elevated level of co-expression of PD-1/4-1BB had improved DSS (P=0.03) ( Figure 6A ) and PFI (P=0.02) ( Figure 6B ) in comparison to the group with low co-expression, with statistically significant differences. However, no statistically significant variance emerged in terms of OS and DFI. Univariable analysis revealed that FIGO stage IV (HR 5.947, 95% CI 2.905-12.177; P < 0.001) and PD-1/4-1BB co-expression (HR 1.988, 95% CI 1.052-3.755; P = 0.034) were risk factors influencing DSS. Similarly, FIGO stage IV (HR 4.095, 95% CI 2.182-7.686; P < 0.001) and PD-1/4-1BB co-expression (HR 1.829, 95% CI 1.073-3.119; P = 0.027) were risk factors affecting PFI in univariable analysis. In multivariable analysis, FIGO stage IV (HR 5.763, 95% CI 2.820-11.815; P < 0.001) and PD-1/4-1BB co-expression (HR 1.953, 95% CI 1.032-3.696; P = 0.04) were independent prognostic factors for DSS. Moreover, FIGO stage IV (HR 3.994, 95% CI 2.126-7.504; P < 0.001) and PD-1/4-1BB co-expression (HR 1.860, 95% CI 1.090-3.176; P = 0.023) were identified as independent prognostic indicators of PFI. To explore the potential molecular mechanisms underlying the differences in DSS and PFI among the groups with high and low co-expression of PD-1/4-1BB, we conducted differential genes ( Figure 6C ) KEGG and GO enrichment analyses. GO analysis indicated enrichment of immune-related cellular functions across Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF). Differentially expressed genes in BP exhibited enrichment in lymphocyte-mediated immunity. CC analysis showed enrichment in T cell receptor complex, while MF analysis indicated enrichment in immune receptor activity ( Figure 6D ). KEGG pathway analysis demonstrated enrichment in pathways related to cytokine-cytokine receptor interaction, viral protein interaction with cytokines and cytokine receptor, hematopoietic cell lineage, and cell adhesion molecules ( Figure 6E ).

To investigate the efficacy of the high and low PD-1/4-1BB co-expression groups as indicators of how patients with CC will respond to different medicines, we calculated IC50 values for 89 drugs. Our results showed that individuals with high expression levels may demonstrate increased sensitivity to drugs such as Mitomycin C, Dasatinib, Ispinesib Mesylate, Gemcitabine, among others, while those with lower expression levels may exhibit enhanced responsiveness to Bleomycin, MS-275, and Etoposide( Figure 6F ). Furthermore, our analysis revealed noteworthy disparities in IC50 values for six commonly employed CC therapeutics, namely Paclitaxel, Gemcitabine, Bleomycin, 5-Fluorouracil, Etoposide, and Mitomycin C ( Figure 6G ).

Given the significance of checkpoint inhibitors in clinical practice, two approaches were employed to validate the predictive capability of PD-1/4-1BB co-expression in immunotherapeutic benefits. Utilizing the IPS, recognized as a reliable predictor of response to ICIs (26), we evaluated the effectiveness of immunotherapy. Next, we conducted a comparison of the IPS values among the cohorts with high and low expression. Elevated IPS scores are indicative of a more favorable response to ICI therapy. Our research findings reveal that in treatments involving CTLA4+/PD1−, CTLA4−/PD1+, and CTLA4+/PD1+, the high-expression cohort consistently exhibited notably higher IPS values compared to the low-expression cohort. This suggests that patients with high-expression demonstrate superior responses to anti-CTLA4, anti-PD-1, and combined anti-CTLA4 and anti-PD-1 therapies compared to those with low-expression. ( Figures 7A-D ). Compared to the low-expression group in the IMvigor210 cohort, a significantly greater percentage of patients in the high expression group achieved complete response/partial response (CR/PR) (P <0.01) ( Figure 7E ). Additionally, PD-1/4-1BB co-expression levels were considerably elevated in the CR/PR subgroup compared to the stable disease (SD)/progressive disease (PD) subgroup (P <0.01) ( Figure 7F ). Consistent with these findings, within the BMS038 cohort, a notably greater percentage of CR/PR patients was noted in the high expression group than in the low expression group (P< 0.05) ( Figure 7G ), with concomitantly higher co-expression levels of PD-1/4-1BB in the CR/PR subgroup compared to the SD/PD subgroup (P<0.01) ( Figure 7H ). Collectively, these findings reinforce the rationale for considering PD-1/4-1BB as a predictive marker of immunotherapy response. They indicate that individuals exhibiting high expression levels may stand to gain enhanced therapeutic outcomes from such interventions.