The clinicopathological association and prognosis analysis of HPV-ind CCs

HPV-ind CCs accounted for 2.48% (25/1010) in our cohort (Fig. 1A). Of all HPV-ind CCs, 11/25 (44%) were squamous cell carcinoma (SCC), 12/25 (48%) were classic adenocarcinoma (ADC), 1/25 (4%) was neuroendocrine tumor (NET) and 1/25 (4%) was peripheral primitive neuroectodermal tumor (PNET). The mean age of patients was 52.8 years old (52.8 ± 14.3). Patient characteristics are detailed in Table 1. HPV-ind CCs patients were older (52.8 vs. 47.6-year-old, P = 0.009) and more likely to be in the post-menopausal stage (60% vs. 36%, P = 0.027). They were more accompanied with parametrial involvement (16% vs. 5%, P = 0.049). ADC was more commonly observed in HPV-ind CCs than HPV-asso CCs patients (48% vs. 15%, P < 0.001). The TCGA dataset included 307 samples, of which 18 were HPV-ind and 289 were HPV-asso (Table 1). The clinical characteristics were basically consistent with those of our cohort except for menopause status, and the proportion of deaths.

Table 1 Patients’ characteristics bewteen HPV-independent and HPV-associated cervical cancers.

In the median follow up of 54 months (range: 1–75 months), HPV status was significantly associated with relapse (HPV-ind: 48% vs. HPV-asso: 18%, P < 0.001). Univariate analysis revealed a striking association between HPV status and patient survival, concordant with previous studies [7, 8, 22], HPV-ind CCs showed a worse outcome than HPV-asso CCs (5-year PFS: 52% vs. 82%, P = 0.001, 5-year OS: 68% vs. 88%, P < 0.001) (Fig. 1B; Table 2). Our multivariate analyses revealed that HPV-ind in cervical cancer was an independent predictor for poorer PFS (HR = 0.53, 95% CI: 0.28–0.99, P = 0.045) (Table 2). Similarly, poorer OS and PFS were observed in patients with HPV-ind CCs from the TCGA cohort, although statistical significance was not achieved. (Fig. 1C, D). Collectively, these results revealed that HPV-ind CCs formed a rare but more deadly type of CCs distinct from HPV-asso CCs.

Table 2 HPV status in cervical cancer predicted independently better PFS in univariate and multivariate analyses.The PI3K pathway was highly activated in HPV-ind CCsGenomic alterations

WES was performed to evaluate the genomic variations in 21 HPV-ind CCs (Fig. 2A). It showed that more PI3K pathway-related mutation were found: PTEN (62%), PIK3CA (52%), and AKT2 (33%). We then analyzed the PIK3CA status of 25 HPV-ind CCs and 972 HPV-asso CCs, which had the detailed mutation information in exon9, and exon20 by Sanger sequencing in our previous study. Compared to HPV-asso CCs, HPV-ind CCs presented higher PIK3CA mutations (36% vs. 13%, P = 0.006, Fig. 2B). Other high-frequency mutations appeared in FGFR2 (48%), FBXW7 (43%) and TP53 (43%) in HPV-ind CCs. Similarly, single nucleotide variation (SNV) from the TCGA dataset was analyzed for validation (Fig. 2C, D). It demonstrated that TTN (32%), PIK3CA (27%), KMT2C (19%), MUC4 (18%), MUC16(17%) were the top 5 most frequent mutated genes in 274 HPV-asso CCs, while among 14 patients with HPV- CESC, PIK3CA (64%), SYNE1 (50%), TP53 (50%), PTEN (36%), TTN (36%) sited the 5 most frequent mutations, depicting a different mutation signature. PIK3CA, EP300, FBXW7, ARID1A, PTEN were identified as significantly mutated genes (SMGs) in HPV-asso CCs, which is consistent with previous findings of CCs overall [10], however, they were NDUFS1, PIK3CA, PTEN and TP53 in HPV-ind CCs. (Table S3, S4). Concordant with our cohort, HPV-ind CCs in the TCGA database had a higher mutation rate of PIK3CA and TP53 (Fig. 2E, PIK3CA 9/14 64% vs. 74/274 27%, p = 0.005; TP53 7/14 50% vs. 15/274 5.5%, P < 0.001 Fisher’s Exact Test). Most mutations of PIK3CA located in the activating helical domain E542K and E545K (Fig. S1A), whereas mutations in TP53 were scattered with no hotspots as reported in both HPV-asso and HPV-ind CESCs (Fig. S1B). No survival differences were observed between PIK3CA mutant and wild-type samples in either HPV-ind CCs or HPV-asso CCs, so as between TP53 mutation status (Fig. S1C–H).

Fig. 2: Mutation alterations of CCs.

A Somatic mutations in HPV-ind CCs from the SHCC cohort. The upper panel shows the number of mutational events for each sample, and the bottom panel shows the transition/transversions type. B Sanger sequencing of PIK3CA mutation in SHCC cohorts. Somatic mutations in HPV-ind CCs (C) and HPV-asso CCs (D) from the TCGA cohort. The middle panel shows the clinical features of each sample. E Stacked bar chart of PI3KCA (left) and TP53 (right) mutations in HPV-asso and HPV-ind CCs from the TCGA cohort. F Comparison of TMB between HPV-ind and HPV-asso CCs from both SHCC and TCGA cohort.

In addition, HPV-ind CCs demonstrated a higher tumor mutation burden (TMB) compared to HPV-asso CCs from the analysis of TCGA cohort (11.11 vs 3.84, p < 0.05, Fig. 2F). Due to lack of genomic information from WES in HPV-asso CCs patients, we analyzed the TMB only in SHCC HPV-ind CCs. It revealed that the TMB of HPV-ind CCs in our cohort was 16.20, which was significantly higher than that of HPV-asso CCs from the TCGA cohort.

Copy number variants (CNVs) were processed and analyzed by GISTIC2.0 in TCGA dataset. An average of 288.5 and 302.2 copy number variants were found for HPV-asso and HPV-ind CESCs, respectively (P = 0.91, Wilcoxon test), but a significant variance of 26.2 and 2.5 CNVs each for HPV-asso and HPV-ind CESCs were achieved after including CNVs restricted to q < 0.1 (Fig. S2A, B, P < 0.001, Wilcoxon test). It revealed 46 amplifications and 36 deletions in HPV-asso CESCs while only 1 amplification and 2 deletions in HPV-ind ones (Fig. S2C, D). The CNVs in HPV-ind CESC included amplification of 3q26.31 (GHSR, FNDC3B, 44.4%) and deletion of 5q11.2 (PDE4D, PART1, 27.8%) and 22q13.32 (MAPK11, TYMP, 27.8%), which displayed different variations with HPV-asso CCs. It is noteworthy that though 3q26.31 appeared in both HPV status, the distribution of its amplification seemed to be complementary to those without PIK3CA mutation in HPV-ind CESCs (Fig. S2E). As GHSR rarely expressed in CESCs, FNDC3B might play an important role in PIK3CA wild type HPV-ind CESCs.

Transcriptome analysis

We further explored the correlation between PIK3CA mutation mediated PI3K pathway activation and HPV status in cervical cancers. RNA-Seq was performed in all 25 HPV-ind CCs and an addition 28 HPV-asso CCs. It was shown that gene expression profiles differed between HPV-asso and HPV-ind CCs. For example, HDAC9, FGFR were highly expressed in HPV-ind CCs (Fig. 3A), which was also confirmed by the TCGA database (Fig. S3A). Gene Set Enrichment Analysis (GSEA) was then performed to display the different patterns of activated pathways between these two types (Figs. 3B; S3B). The PI3K/AKT signaling pathway was significantly enriched in HPV-ind CCs in both cohorts (Fig. 3C), while the P53 signaling pathway was highly activated in HPV-asso CCs (Fig. 3D; Fig. S3C), suggesting that the PI3K/AKT signaling, instead of TP53 signaling, may strongly engage in the oncogenic effects of HPV-ind CCs. Concurrently, the fibroblast growth factors (FGFR) family, including FGFR1 and FGFR4, was highly enriched in HPV-ind CCs (Fig. 3E; Fig. S3D). Other enrichment pathways between these two CCs included the abnormal activation of G-protein coupled receptor (GPCR) and the calcium signaling pathway in HPV-ind CCs, as well as hyperactivation of the senescence-associated secretory phenotype (SASP), G2M checkpoint and E2F related signaling in HPV-asso CCs. To undercover the downstream effects of activated PI3K/AKT signaling, a gene network of FGFR-PI3K-AKT signaling pathway was constructed (Fig. 3F). In addition to the receptors FGFR1, FGFR4 and ERBB4, the main effector AKT1 and its downstream effectors, MTOR, GSK3B, p21, WEE1, BAD and CREB5, which mediate cancer cell growth and progression, also exhibited marked changes in HPV-ind CCs. Collectively, the aberrant activation of the PI3K/AKT pathway may act as a potential driver of cervical carcinogenesis in the absence of HPV infection, thus indicating the potential role of PI3K inhibitors. Furthermore, the activation of the FGFR pathway suggested a combined therapeutic role of PI3K and TKI inhibitors in HPV-ind CCs.

Fig. 3: Transcriptome analysis revealed high activation of the PI3K/AKT pathway in HPV-ind CCs.

A Heatmap of significantly differential expression in cervical cancers with several of the most variable mRNAs annotated on the right. B Bar plot shows cancer-related pathways enriched by GSEA between HPV-ind and HPV-asso CCs in SHCC cohort. GSEA plot displays that PI3K/AKT signaling was enriched in HPV-ind CCs in both cohorts (C) while TP53 signaling pathway was enriched in HPV-asso CCs (D). FGFR pathways were also enriched in HPV-ind CCs (E). F Relative levels of transcript of main effectors in the PI3K/AKT signaling pathway.

Immune infiltration was suppressed in HPV-ind CCs

The tumor microenvironment (TME) was analyzed between HPV-ind and HPV-asso CCs in two cohorts. Sixty-four cell types in these two groups were quantified by the xCell algorithm, and the differences are shown in Fig. S4A and Fig. S4B. HPV-ind CCs contain more matrix cells, such as the pericytes, fibroblasts and preadipocytes. Additionally, compared to HPV-asso CCs, common immune cells such as lymphoid progenitor (CLP), pro B-cells and CD8+ naïve T-cells were less infiltrated in HPV-ind individuals. Next, immune-related pathways were further analyzed by ssGSEA (Fig. S4C, D). As expected, antigen processing- and presentation- related immune-promoted pathway were significantly downregulated, while members of the immune-inhibited pathway (TGF-β family members and their receptors) were upregulated in HPV-ind CCs. Furthermore, BCR, TCR signaling pathway, TNF family members receptors, chemokines and interferons were downregulated in HPV-ind CCs from TCGA cohort. (Fig. S4D). We next merged the two cohorts in calculating the immune cell infiltration by the CIBERSORT algorithm to acquire a higher confidence. Apparently, the HPV-ind CCs had less immune cell infiltration than the HPV-asso ones either in terms of innate immune cells or adaptive immune cells, including plasma cells, CD8 + T cells, activated CD4+ memory T cells, activated NK cells, dendritic cells, and activated macrophages, accompanied by enrichment of resting immune cells, anti-immune M2 microphage, mast cells and neutrophils. (Fig. S4E–G). Hence, a relatively suppressed immune state was exhibited in HPV-ind CCs, which also explained the worse prognosis, and indicated that HPV-ind CC patients may hardly benefit from immunotherapies such as immune checkpoint inhibitors or cancer vaccines due to the barren TME.

PI3K-α inhibitor BYL719 has a superior effect on HPV-ind cervical cancer cells

Considering the significant role of aberrant activation of the PI3K/AKT signaling pathway, the inhibitory effects of PI3K α-selective inhibitor BYL719 were evaluated in HPV-ind (C-33A) and HPV-asso cervical cancer cell lines (SiHa and MS751). It revealed that C-33A cells were more sensitive to BYL719 than MS751 and SiHa cells (IC50: 19.94 nmol/L (C-33A) vs. 26.69 nmol/L (SiHa); 38.39nmol/L (MS751), Fig. 4A). A concentration of 25nmol/L BYL719 significantly inhibited C-33A cell proliferation but barely affected that of SiHa and MS751 (Fig. 4B, C). Additionally, BYL719 treatment contributed to higher apoptosis rate (23.7% (C-33A) vs. 11.67% (SiHa); 7.5% (MS751), p < 0.001, Fig. 4D, F) as well as more cells arrest in the G0/G1 phase in C-33A cells than that of SiHa and MS751 cells (Fig. 4E, G). Consistently, corresponding changes of apoptosis, proliferation and cell cycle-related proteins were observed in C-33A cells (Fig. 4H). Increased phosphorylation of downstream substrates of the PI3K/AKT signaling pathway was significantly repressed on the treatment of BYL719 in C-33A cells (Fig. 4I), supporting our hypothesis of the aberrant activation of the PI3K/AKT signaling pathway and the inhibitory effect of BYL719 in HPV-ind CCs.

Fig. 4: PI3K-α inhibitor BYL719 has a superior effect on HPV-ind CC cells.

A CCK-8 assay was analyzed with BYL719 treatment for 72 h. Cells were treated with 25uM BYL719 for 10 days and then stained with crystal violet (B). Colonies with over 50 cells were quantified (C). D Cells were treated with BYL719 for 72 h and then stained with Annexin-V and 7-AAD. Apoptosis proportion was analyzed by flow cytometry. E Cells were treated with BYL719 for 72 h and then stained with PI. Cell cycle distribution was analyzed by flow cytometry. Proportion of apoptosis cells (F) and cell cycle (G) was presented as histogram. (H–I) The expression level of apoptosis, cell cycle, damage repair marker (H) and downstream substrates of PI3K-Akt signaling pathway (I) after treatment of BYL719 for 72 h determined by western blotting. All data are presented as means ± SD of three independent experiments. *P < 0.05; **P < 0.01. Wilcox test.

BYL719 showed marked anti-tumor effect on HPV-ind cervical cancer patient derived models

Next, we evaluated BYL719 in two HPV-ind and two HPV-asso patient-derived organoids (PDOs) (Fig. 5A). Consistent with the in vitro experiments, HPV-ind PDOs had lower IC50 as compared to HPV-asso PDOs (0.59 ± 0.22 nmol/L vs. 3.54 ± 0.16 nmol/L, P = 0.013, Fig. 5B). Similarly, the antitumor effects of BYL719 were determined in two patient-derived tumor xenograft (PDX) models from HPV-ind and HPV-asso patients. The xenografts were treated with BYL719 for 28 days. It showed that tumor proliferation was significantly repressed with a higher tumor growth inhibition rate (TGI%) in HPV-ind PDXs than that with HPV-asso xenografts (Fig. 5C–F), however, an obvious body weight loss was also observed (Fig. 5G). Furthermore, tumors of HPV-ind PDXs treated with BYL719 displayed reduced staining of Ki67, and increased staining of cleaved caspase3 and γH2AX, indicating that BYL719 had superior effect on inhibiting tumor proliferation, promoting tumor apoptosis and DNA damage in HPV-ind PDXs than that with HPV-asso PDXs (Fig. 5H).

Fig. 5: BYL719 showed marked anti-tumor effect on HPV-ind CC patient derived models.

A Quantification of cell death in four HPV-ind and HPV-asso PDOs after 120 h of treatment with different concentrations of BYL719. B CCK-8 assay of HPV-ind and HPV-asso PDOs were analyzed with BYL719 treatment for 120 h. Tumor growth (C), tumor weight (D) and TGI% (E) curve of HPV-ind and HPV-asso PDX models treated with BYL719 p.o. for 28 days. F Comparison of the gross appearances of HPV-ind PDX and HPV-asso PDX between control and BYL719 treated groups (n = 2 samples for each group). G The body weight curve showed a decrease in HPV-ind PDX after 7 days of treatment with BYL719. H Representative images of Ki67, cleaved caspase3 and γH2AX staining of HPV-ind and HPV-asso PDX tissues. I GSEA for untreated PDXs showing enrichment of the PI3K/AKT signaling pathway in HPV-ind PDX. GSEA for PDXs treated with BYL719 p.o. for 28 days showing different patterns of changed tumor-related signaling pathways. Top 5 significant enriched pathways were shown for HPV-ind PDX (J). All 3 significant enriched pathways were shown for HPV-asso PDX (K). All data are presented as means ± SD of three independent experiments. *P < 0.05; **P < 0.01. Wilcox test.

To further understand the effect of BYL719 on TME in HPV-ind and HPV-asso PDXs, RNA-Seq was performed with tumor tissues. It identified that the PI3K-Akt signaling pathway was highly enriched in HPV-ind PDXs compared to HPV-asso PDXs in untreated groups, although PIK3CA mutation was not carried (Fig. 5I). While it was evidently inhibited under the treatment of BYL719, the treatment also influenced the other cancer-related pathways, such as JAK-STAT, NF-kB, Ras, TNF signaling pathways (Fig. 5J). However, only changes in calcium, TNF and TGF-beta signaling pathways were found after the exposure of BYL719 in HPV-asso PDXs (Fig. 5K). In summary, these results revealed that the PI3K-Akt signaling pathway was aberrant activated in HPV-ind CCs, and these patients may benefit from the PI3K α-selective inhibitor BYL719.

Molecular subgroups of HPV-ind CCs were identified by RNA-Seq

To investigate the heterogeneity among HPV-ind CCs, our RNA-Seq data was analyzed using unsupervised hierarchical clustering to classify 25 HPV-ind CCs into different subgroups. The optimal number (K = 2) of clusters was chosen by ConsensusClusterPlus package. (Fig. 6A and Fig. S5A). Survival analyses of the two clusters illustrated an obviously lower survival probability for subtype 1 in terms of both OS (HR: 5.9, CI: 1.3–27.2, P = 0.018, Fig. 6B) and PFS (HR: 8.4, CI:1.8-39.1, P = 0.0044, Fig. 6C). Subsequently, the same features were also validated in the TCGA cohort (Fig. S5B, C). Survival analyses showed a difference in survival probability between the two clusters, albeit with no statistical significance in the log-rank test (Fig. 6D, E). Principal component analysis (PCA) conducted managed to differentiate the clusters in both cohorts (Fig. S5D, E), and the clinical features of the two clusters in SHCC and TCGA cohorts were depicted in Table S5

Fig. 6: Molecular clustering of HPV-ind CCs.

A Heatmap with clinical features annotated above of 2 clusters from the SHCC cohort identified by unsupervised hierarchical clustering. Kaplan-Meier plot in two clusters of OS (B) and PFS (C) in the SHCC cohort, OS (D) and PFI (E) in the TCGA-CESC cohort. F Bar plot shows cancer-related pathways enriched by GSEA between two subtypes in merged cohort. The representative GSEA plot shows that the fatty acid metabolism was enriched in subtype 1 (G) while the PD-1 signaling and signaling by NOTCH were enriched in subtype 2 (H). I Bar plot displays different enrichments in TILs and immune-related pathways between two subtypes in merged cohorts. NS no significance. *P < 0.05; **P < 0.01. Wilcox test.

Further independent analysis of both cohorts showed that subtype 1, named metabolic subtype, was more active in metabolism and epithelial-mesenchymal transition, while immune-related pathways were distinctly enriched in subtype 2, named immune subtype (Fig. S5G, H). The results were consistent with those of separate analyses when the samples from two cohorts were combined (Figs. 6F and S5I). In detail, the metabolic subtype was enriched in fatty acid metabolism while the immune subtype was enriched with targets of novel antitumor targeted therapies, such as PD-1 signaling and Notch signaling (Fig. 6G, H), In addition, PIK3CA was found highly mutated in immune subtype in the SHCC cohort (15% (subtype1) vs. 58% (subtype2), P = 0.041, Fisher exact test), but there were no differences in PI3K-AKT signaling activation between these two subtypes (Fig. S5F), suggesting an alternative PI3K-AKT signaling activating mechanism beyond the PIK3CA mutation in the metabolic subtype. Furthermore, the immune subtype had more tumor-infiltrating lymphocytes (TILs) especially CD8+ and activated CD4 + T cells, and higher immune score in most immune-promoting pathways, including antigen processing and presentation, chemokines, interferons, and TNF-related pathways, alongside lower scores for the immune-suppressive TGF-β pathway. (Fig. 6I).

Collectively, we identified two subtypes with distinct biological characteristics and outcomes in HPV-ind CCs, which indicated the heterogeneity of HPV-ind CCs and personalized therapies should be suggested in different subtypes.