HPV-driven heterogeneity in cervical cancer: study on the role of epithelial cells and myofibroblasts in the tumor progression based on single-cell RNA sequencing analysis

Data acquisition

The scRNA-seq data including five adenocarcinoma of the cervix (ADC) samples in the GSE197461 cohort were obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/) using the 10x Genomics and Illumina NovaSeq 6000 platforms for sequencing (Qu et al., 2023). Specifically, the ADC_1, ADC_2 and ADC_3 are human papillomavirus positive (HPV+) groups, while the ADC_4 and ADC_5 are the HPV− groups.

Analysis of scRNA-seq landscape of the ADC samples

We performed the scRNA-seq analysis by using the Seurat R package (Zulibiya et al., 2023). The Read10X function was used to read the expression matrix of each sample, and the cells with 200~7,000 genes expressed in at least three cells and <10% proportion of mitochondrial gene were retained. Then, the SCTransform function was used to perform data normalization and the harmony R package was applied to remove the batch effects between samples after principal component analysis (PCA) dimension reduction. The RunUMAP function further was used for the Uniform Manifold Approximation and Projection (UMAP) dimension reduction and visualization. Cell clustering (setting dims = 1:30, resolution = 0.1) was performed employing the FindNeighbors and FindClusters functions to subdivide epithelial cells and myofibroblasts at the dims = 1:20 and resolution = 0.2. Cell markers and annotation data were obtained from the CellMarker2.0 database.

Gene set enrichment analysis

To elucidate the enrichment difference in the functions of epithelial cells and myofibroblasts in HPV+ and HPV− groups, the FindMarkers function of Seurat R package was used to identify significantly differential genes between the two groups based on significant log₂ fold change value (Reimand et al., 2019). Subsequently, the gseGO function in the ClusterProlifer R package was used for the biological process (BP) analysis of gene ontology (GO) term and the gseKEGG function was used for the Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. Visualization of the pathways was achieved applying the GseaVis R package (Reimand et al., 2019).

Analysis of differentially expressed genes

The differentially expressed genes (DEGs) analysis among each cell subgroups was performed by using the FindAllMarkers function (setting only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25) (Song et al., 2023). After that, these DEGs were uploaded into the Database for Annotation, Visualization and Integrated Discovery database (DAVID, https://david.ncifcrf.gov/) for the enriched BPs (Xiang et al., 2021).

Copy number variation analysis

The inferCNV R package was used for the copy number variation (CNV) analysis for epithelial cells and myofibroblasts in HPV+ and HPV− groups, in which the B cells in the healthy group served as a reference (setting cluster_by_groups = TRUE, analysis_mode = “subclusters”, HMM_type = “i3”, denoise = TRUE, HMM_report_by = “subcluster”, HMM=TRUE) (Lu et al., 2023).

Cell culture and intervention

The human normal cervical endometrial cell line HUCEC (CP-H059) and the CC cell line Hela (CL-0101) were both obtained from Procell Lifescience Co., Ltd. in Wuhan, China. Cells were cultured either in Roswell Park Memorial Institute-1640 medium (PM150110; Procell Lifescience Co., Ltd., Wuhan, China) or in Minimal Essential Medium (PM150410; Procell Lifescience Co., Ltd., Wuhan, China) at 37 °C with 5% CO₂. All the media were additionally supplemented with 10% fetal bovine serum (FBS, 164210; Procell Lifescience Co., Ltd., Wuhan, China) and 1% antibiotics (PB180120; Procell Lifescience Co., Ltd., Wuhan, China).

The small interfering RNAs for TFF3 (target sequence: 5′-AAGCAGAAAAAATACATTTCAGG-3′) and the corresponding scramble negative control (5′-GGATAGCAGAACAAAATTTCAAA-3′) were ordered from GenePharma (Shanghai, China) and transfected into Hela cells using Lipofectamine 2000 transfection kit (11668-027, Invitrogen, Carlsbad, CA, USA) following the manuals.

Cell migration and invasion assays

A total of 1 × 10⁵ Hela cells were plated in 6-well dishes in a culture medium containing 10% FBS to complete confluence. Subsequently, the cell monolayers were wounded with sterile 10 μL pipette tips, and the debris was eliminated by rinsing in PBS. Following this, the wounded monolayers were incubated in serum-free medium for 48 h (h). The wound healing rate was observed using an inverted optical microscope (DP27, Olympus Corp., Tokyo, Japan) and analyzed by ImageJ (version 1.48, National Institute of Health, Bethesda, MD, USA).

For the invasion assay, transwell chamber with 8 μm pore size culture inserts (3422; Corning, Inc., Corning, NY, USA) was used. In detail, a total of 1 × 10⁶ cells/mL in 200 μL serum-free culture medium were placed in the upper chamber coated with Matrigel (354230, Corning, Inc., Corning, NY, USA), while the corresponding lower chamber was supplemented with 700 μL complete medium containing 10% serum. After 48 h, the cells invaded to the lower chamber were immobilized with 4% fixation solution (P1110; Solarbio, Beijing, China) and dyed with 0.1% crystalline violet staining (G1063, Solarbio, Beijing, China) for 20 min (min). The number of invaded cells at 48 h was calculated under the same optical microscope and processed with ImageJ.

Real-time quantitative PCR

TRIzol reagent (15596-026; Invitrogen, Waltham, MA, USA) was applied to isolate and prepare total cellular RNA, which were reverse-transcribed into corresponding complementary DNA using a related assay kit (K16325; Solarbio, Beijing, China) according to the protocol. RT-qPCR was then performed using SYBR Green PCR Mastermix (SR1110; Solarbio, Beijing, China) at the following parameters: at 95 °C for 3 min, and 40 cycles at 95 °C for 20 s (s), at 60 °C for 30 s and at 72 °C for 60 s. GAPDH was an endogenous control, and the relative mRNA level was calculated with the 2^−ΔΔCT method (Livak & Schmittgen, 2001). The primers (5′–3′) involved were listed below:

APOC1: Forward: AGGACAAGGCTCGGGAACTCAT; Reverse: GATGTCACCCTTCAGGTCCTCA.

CEACAM6: Forward: GCCTCAATAGGACCACAGTCAC; Reverse: AGGGCTGCTATATCAGAGCGAC.

FOXP3: Forward: GGCACAATGTCTCCTCCAGAGA; Reverse: CAGATGAAGCCTTGGTCAGTGC.

SFRP4: Forward: CTATGACCGTGGCGTGTGCATT; Reverse: GCTTAGGCGTTTACAGTCAACATC.

TFF3: Forward: TCCAGCTCTGCTGAGGAGTACG; Reverse: ATCCTGGAGTCAAAGCAGCAGC.

GZMA: Forward: CCACACGCGAAGGTGACCTTAA; Reverse: CCTGCAACTTGGCACATGGTTC.

GAPDH: Forward: GTCTCCTCTGACTTCAACAGCG; Reverse: ACCACCCTGTTGCTGTAGCCAA.

Statistical analysis

All calculations were performed using the R software (version 4.3.1) or GrpahPad Prism software (version 8.0.2). For bioinformatics analysis, the data between groups were compared with Wilcoxon’s test, while student’s t test was applied for cell validation test. A p < 0.05 was considered as statistically significant.

Epithelial cell and myofibroblasts had high infiltration levels in HPV+ group

After quality control, normalization and UMAP dimensionality reduction, a total of seven main cell subgroups were obtained from 12,431 cells, including epithelial cells, myofibroblasts, plasma B cells, cytotoxic NK/T cells, regulatory T cells, B cells and myeloid cells (Fig. 1A). Cell clusters (Figs. S1A and S1B), gene numbers (Fig. S1C) and mitochondrial ratio (Fig. S1D) in different samples were visualized. The cluster-specific marker genes were shown in Figs. 1B and 1C, from which it could be observed that NKG7, KRT8, AGR2 and CLDN4 were significantly expressed in cytotoxic NK/T cells and epithelial cells. Immune cells such as regulatory T cells, cytotoxic NK/T cells, myeloid cells, plasma B cells and B cells had higher proportion in the HPV− group, while the epithelial cell and myofibroblasts had higher proportion in the HPV+ group (Figs. 1D and 1E). This indicated that the patients in the HPV− group with a high immune infiltration had a stronger immunity against cancer, while the immune ability in the HPV+ groups may be inhibited, and that the epithelial cells and myofibroblasts may play an important role in cervical cancer progression.

Epithelial cells in the HPV+ group contributed to antigen presentation with a lower cell proliferation activity

Furthermore, we explored the differences of gene expression patterns of epithelial cells between the HPV+/− groups. The FindMarkers was used to identify significant genes, and GO enrichment analysis found that cytoplasmic translation, catecholamine secretion and transport process were enhanced in the HPV− group, while the HPV+ group has significantly enriched processes of nuclear division, chromosome division, positive regulation of cell cycle and detection of chemical stimulation (Fig. 2A) as well as remarkably upregulated HES1, SPHK1 and UBE2C in the positive regulation of cell cycle pathway (Fig. 2B). KEGG analysis revealed that the activities of ribosome, oxidative phosphorylation, antigen processing and presentation pathways in HPV− group were significantly enhanced, while the activity of cell cycle, cAMP signaling pathway, DNA replication, P53 pathway and steroid hormone pathway were noticeably inhibited (Fig. 2C). In HPV− group, the expressions of CDKN1A, SERPINB5 and TP53I3 in the P53 signaling pathway were downregulated (Fig. 2D). Activation of antigen-presenting pathways may be associated with immunotherapeutic response (Lee et al., 2020), and dysregulation of P53 signaling pathway is the focus of many targeted therapies (Liebl & Hofmann, 2021). To this end, we also observed that the marker genes of antigen processing and presentation pathway such as PSEM1, HSP90AB1 and HSPA1B were significantly upregulated in HPV− group (Fig. 2E), while the marker genes of P53 signaling pathway such as CDKN (1A and 2A), GADD45B and PERP were significantly upregulated in the HPV+ group (Fig. 2F), suggesting that epithelial cells in the HPV− group were involved in antigen processing and presentation with a lower activity of cell proliferation, which resulted in the depletion of epithelial cells.

Amplification of chr19p and chr20 may promote the epithelial cell proliferation in the HPV+ group

Analysis on the CNV of epithelial cells revealed amplifications in the chr1q and chr3q and deletions on the chr3p and chr6p between HPV+ and HPV− groups, and the difference was that the HPV− group had an amplification on the chr8q and a deletion on the chr17p (Fig. 3A). Specially, in the HPV− group, most negative regulation of epithelial cell proliferation genes such as RIDA, TP53INP1, STK3 and PEX2 were enriched on the chr8q amplification region (Fig. 3B), and the deletion region of chr17p contained genes associated with the MAPK, JNK cascade-related cell proliferation, for example, MAP2K3, PSMB6, GPS2 (Fig. 3C). These may help explain a low cell proliferation activity of epithelial cells in the HPV− group. In addition, the amplification region of chr19p and chr20 in the HPV+ group contained genes such as AURKA, BLCAP, CDC25B and MBD3, which were correlated with the chromatin remodeling, cell cycle, cell migration, and viral resistance response (Figs. 3D and 3E) that contributed the increase in the proportion of epithelial cells.

Activated cell cycle and proliferation pathways enhanced the heterogeneity of epithelial cells in HPV+ group

The heterogeneity of epithelial cells in the HPV+ and HPV− groups were analyzed and the epithelial cells were further subdivided into three subgroups (epithelial cells 1, 2 and 3) (Fig. 4A) and their marker genes were showed in the (Fig. 4B). Interestingly, the epithelial cells 1 had higher proportion in the HPV− group, while the epithelial cells 2 and epithelial cells 3 had higher proportion in the HPV+ group (Fig. 4C). We analyzed the functions of significantly high-expressed genes in the three epithelial cell subgroups and found that cell cycle, cell migration, angiogenesis, viral resistance response and hypoxia pathways were significantly activated in the epithelial cells 1 (p < 0.05, Fig. 4D), while pathways such as oxidative stress, glycolysis and gluconeogenesis were mainly enriched in the epithelial cells 2 (p < 0.05, Fig. 4E) and insulin response, cold and heat response, progesterone response and corticosterone response pathways were enriched in the epithelial cells 3 (p < 0.05, Fig. 4F). These results indicated that the epithelial cells in the HPV+ group had higher cell proliferation and migration, resistance to viruses and angiogenesis activity. Correspondingly, enhanced oxidative stress and glycolysis, insulin, progesterone, and corticosterone activity all enhanced the tumor immunity in the HPV− group. This indicated that different epithelial cell subpopulations may play different roles in the progression and malignant transformation of CC by participating in specific biological processes.

Activated tumor immunity increased the depletion of myofibroblasts in the HPV− group

Analysis on the gene expression pattern in the myofibroblasts in the HPV+ and HPV− groups identified significantly expressed genes. The GO enrichment analysis showed that in the HPV− group, immunoglobulin production, T- and B-cell-mediated immunity, antigen processing and presentation pathways were mainly activated, but cAMP cellular response, epidermal morphogenesis and adaptability of striated muscle and skeletal muscle pathways were significantly inhibited (Fig. 5A), HLA-F, DRB1, CTSH and HLA-C, A and E involved in T cell-mediated cytotoxicity were upregulated (p < 0.05, Fig. 5B). KEGG enrichment analysis showed that pathways including antigen processing and presentation, autoimmune thyroid disease and systemic lupus erythematosus, human papillomavirus infection were significantly activated in the HPV− group, while steroid hormone production, chemical cancer generation pathways were remarkably suppressed (p < 0.05, Fig. 5C). The genes (CYP1B1, EPHX1, GSTM5 and PTGS2) involved in the chemical carcinogenesis-DNA adducts in were noticeably upregulated the HPV+ group (p < 0.05, Fig. 5D). The bubble plot showed that the genes of the antigen processing and presentation (HLA-A, B, C and E) were significantly expressed in the HPV− group (Fig. 5E) and the genes of chemical carcinogenesis-DNA adducts (GSTM3 and 2, EPHX1) were significantly expressed in the HPV+ group (Fig. 5F), indicating that tumor immune activation may accelerate the depletion of myofibroblasts in the patients in the HPV− group.

Activation of anti-cancer immunity in myofibroblasts 1 in the HPV− group

Analysis on the CNV of myofibroblasts revealed an extensive deletion of on the chr6p in the two HPV groups, indicating that there was no obvious difference of CNV in myofibroblasts between HPV+ and HPV− groups. The myofibroblasts were further subdivided into the subgroups of myofibroblasts 1 and 2 (Fig. 6A), with the HPV+ group having a higher proportion of myofibroblasts 2 (Fig. 6B). Pathway enrichment analysis revealed that cell migration, positive regulation of interleukin-8 production, immune system process and positive regulation of T cell-mediated cytotoxicity were remarkably enriched in the myofibroblasts 1 (Fig. 6C), while signal transduction, cell adhesion, cell growth and proliferation, angiogenesis pathways were mainly enriched in the myofibroblasts 2 (Fig. 6D). The genes (HLA-E, A and B, B2M) involved in positive regulation of T cell-mediated cytotoxicity pathway were significantly expressed in the myofibroblasts 1 (p < 0.05, Fig. 6E), while the marker genes of positive regulation of cell division pathway including the PGF, PTN, MDK and PGF7 were significantly expressed in the myofibroblasts 2 (p < 0.05, Fig. 6F). These results suggested that the activated cell growth and proliferation, angiogenesis pathways could increase the number of myofibroblasts 2 in the HPV+ group, while activated immune response pathway of myofibroblasts 1 contributed to the anti-cancer immunity and promoted the depletion of myofibroblasts.

In-vitro cell validation

Due to the importance of method validation in biomarker analysis (Seyfinejad & Jouyban, 2022), the reliability of our identified marker genes based on HPV+ and HPV- samples were validated. To further validate the reliability of our identified marker genes based on HPV+ and HPV− samples, we assessed the levels of markers for specific cellular subpopulations including APOC1, CEACAM6, FOXP3, SFRP4, TFF3, and GZMA in CC cell line Hela and normal cervical endometrial cell line HUCEC. As shown in Figs. 7A–7F, we observed that except GZMA, the mRNA levels of all the rest markers were expressed at higher levels in Hela cells than in HUCEC cells (p < 0.01). Also, the results from wound healing assay (Figs. 7G and 7H) and Transwell assay (Figs. 7I and 7J) showed that the silencing of TFF3 visibly suppressed the migration and invasion of Hela cells in vitro (p < 0.001; Figs. 7G–7J).