Single-cell sequencing data were obtained by accessing the ArrayExpress (https://www.ebi.ac.uk/biostudies/arrayexpress) database with CC number E-MTAB-12305. The research data were derived from public databases, so no ethical approval was required.

The raw data were imported into R software (version 4.3.2) using Seurat R package for analysis. The DoubletFinder function (28–30) was used to perform strict quality control on the scRNA-seq data to filter low-quality cells and remove potential doublet cells. The quality control standards were as follows: 300<nFeature<6000; 500<nCount<50000; 0<pMT<25; 0<pHB<5.

The data were normalized using the NormalizeData function (31–35). The FindVariable Features function was used to find highly variable genes (36–40). The ScaleData R package was used for standardization processing of the data. Principal component analysis (PCA) (41–43) was performed on the data, and the top 30 principal components were selected. The harmony R package (44, 45) was used to remove batch effects.

The cells were clustered using the FindClusters function and the FindNeighbors function (28, 46), and the cells were annotated according to the CellMarker database (http://xteam.xbio.top/CellMarker/) and cell-typical markers. The FindAllMarkers function was used to find differentially expressed genes in different cell subpopulations.

Taking ECs as reference, the inferCNV algorithm was used to calculate the copy number variation of EPCs.

Gene ontology (GO)enrichment analysi (47–52) of differentially expressed genes was performed using ClusterProfiler R package (51, 53, 54). Gene set enrichment analysis (GSEA)was performed by downloading GSEA software from the GSEA website (http://software.broadinstitute.org/gsea/msigdb).

The getlineage function in the Slingshot R package is used to infer cell differentiation lineages (55–57). The getCurves function is used to calculate the cell expression levels of different lineages within the fitting time.

CellChat R package (58–60) was used to infer the interactions between cells. The identifyCommunicationPatterns algorithm is used to count the number of communication patterns, and the netVisual_diffInteraction algorithm is used to calculate the difference in the strength of communication between cells. Circle plots, violin plots, layer plots, and heat maps are used to visualize the communication between cells.

The activity of transcription factors in EPCs subsets was assessed using the SCENIC package in Python (version 3.7). The GRNBoost2 algorithm was used to infer co-expression modules between transcription factors and target genes. The RcisTarget algorithm was used to analyze the genes in each co-expression module to help identify enriched motifs. The AUCell score was used to evaluate the activity of transcription factors in EPCs.

HeLa cell lines and SiHa cell lines were propagated in MEM medium. Cells were seeded at 50% density in 6-well plates and then transfected with FOXM1-specific knockdown constructs (si-FOXM1–1 and si-FOXM1-2) and negative control constructs (si-NC). Transfection was performed using Lipofectamine 3000RNAiMAX (Invitrogen, USA) according to the manufacturer’s instructions. siRNA sequences: si-1: AAGAAGAAAUCCUGGUUAA; si-2: ACUAUCAACAAUAGCCUAU. qRT-PCR primers: F: AAACCTGCAGCTAGGGATGT; R: AGCCCAGTCCATCAGAACTC.

After being plated in 96-well plates, the cells were cultured for a full day. After adding 10 μL of CCK-8 labeling reagent to each well, the wells were left in the dark for two hours. Using an enzyme marker, absorbance at 450 nm was used to measure cell viability.

After the cells were starved for twenty-four hours, the cell suspension was combined with Matrigel and introduced onto the Costar plate’s upper chamber. The lower chamber was filled with the serum-containing media. The cells were fixed with 4% paraformaldehyde and stained with crystal violet following a 48-hour incubation period.

The cell monolayer in each well of a six-well plate was consistently scratched using a sterile 200 μL pipette tip, and scratch photos were recorded after 0 and 48 hours of incubation. Scratch width was measured using Image-J software.

After adding the EdU working solution, the cells were incubated for two hours. Following a PBS wash, the cells were fixed using 4% paraformaldehyde solution, permeabilized and quenched using a solution containing 0.5% Tr and 2 mg/ml glycine, and stained using 1X Apollo solution and Hoechst staining reaction solution.

All data were processed using R or Python. P< 0.05 was considered statistically significant.

To explore the evolution of transcriptome features during CC development, we performed scRNA-seq analysis on 9 samples from normal cervix (NAT), HSIL and cervical cancer (Tumor), and finally identified 85,591 high-quality cells, as shown in Figure 1A . Using known cell type marker genes ( Figure 1B ), we finally identified 10 types of cells in these samples, including T cells and NK cells, ECs, fibroblasts, smooth muscle cells (SMCs), EPCs, B cells, plasma cells, mast cells (MCs), neutrophils and myeloid cells As shown in Figure 1C , EPCs are the main cell type present in tumor samples, accounting for 65.3%, and are mainly in the G2M and S phases of the cell cycle, with active cell replication and vigorous proliferation ( Figure 1D ). Figure 1E confirms this result. In addition, we also looked at the number of molecules and total number of genes detected in different cell types, and the results showed that EPCs had higher nCount-RNA and nFeature-RNA expression ( Figure 1F ). This suggests that EPCs have higher cell activity and rich gene expression diversity.

As the precursor of tumor cells, the molecular characteristics of EPCs during cancer evolution are crucial to reveal the origin and development of tumors. Therefore, we conducted an in-depth exploration of the characteristics of EPCs during CC progression. First, we used inferCNV to analyze EPCs to view their copy number variations ( Figure 2A ). Next, we re-clustered EPCs and identified 6 EPCs subtypes: C0 SPRR1B+ EPCs, C1 IGFBP7+ EPCs, C2 MUC5B+ EPCs, C3 TOP2A+ EPCs, C4 CFAP126+ EPCs, and C5 PTPRC+ EPCs ( Figure 2B ). Figure 2C shows the expression of the top 5 marker genes of the 6 EPCs subtypes and different groups. The expression levels of the named genes of the 6 EPCs subtypes are shown in bar graphs and UMAP graphs ( Figures 2D, E ). It is worth mentioning that we found that C3 TOP2A+ EPCs had the highest expression of G2M and S phase cell cycle phase scores ( Figure 2F ), which means that this EPCs subpopulation has active proliferation and may play an important role in the rapid progression of tumors. In-depth analysis found that C3 TOP2A+ EPCs tend to be distributed in tumor samples ( Figure 2G ), suggesting that this EPCs subpopulation may be tumor-related EPCs.

To verify the above results, we used Slingshot to infer the differentiation trajectory of EPCs. As shown in Figure 3A , in the differentiation lineage 1 of EPCs, C3 TOP2A+ EPCs are located at the end of the differentiation trajectory, which is consistent with the results of Figure 2G , that is, the PCs subpopulation is tumor-associated EPCs and is of key significance in the evolution of CC. While in lineage 2, C4 CFAP126+ EPCs are located at the end. Figures 3B, C show the differences in the expression levels of the named genes of the six EPCs subtypes over time, among which TOP2A is mainly expressed at the end of lineage 1.

To further explore the biological characteristics of C3 TOP2A+ EPCs, we performed differential gene analysis and enrichment analysis on different subtypes of EPCs ( Figures 4A, B ). As shown in the figure, C3 TOP2A+ EPCs are mainly enriched in chromosome segregation, mitotic nuclear division, nuclear chromosome segregation, mitotic sister chromatid segregation, sister chromatid segregation and nuclear division-related pathways, which are closely related to chromosome division and mitosis. In addition, GSEA further confirmed the above findings ( Figure 4C ). It is worth mentioning that with the progression of CC, the metabolic activity of EPCs gradually increased, especially in pathways such as riboflavin metabolism, pyruvate metabolism, oxidative phosphorylation, and glycolysis/gluconeogenesis. Consistent with this, the metabolic activity of these pathways was also significantly enhanced in C3 TOP2A+ EPCs ( Figure 4D ).

To reveal the crosstalk relationship between C3 TOP2A+ EPCs and other cells, we first analyzed the number and strength of interactions between EPCs and other cell types using cellular communication ( Figure 5A ). Further research found that when C3 TOP2A+ EPCs serve as signal senders, they have a strong interaction with immune cells, such as T cells and NK cells, B cells and myeloid cells, suggesting that they are closely related to the shaping of the tumor immune microenvironment ( Figure 5B ). When C3 TOP2A+ EPCs serve as signal receptors, they interact strongly with stromal cells, such as fibroblasts and SMCs, affecting the matrix shaping of the tumor microenvironment ( Figure 5C ). This suggests that C3 TOP2A+ EPCs have significant plasticity.

Analysis found that C3 TOP2A+ EPCs interacted with other cells mainly through the LAMININ signaling pathway, especially the LAMC1 - (ITGA3+ITGB1) ligand receptor pair ( Figures 5D-F ). In this pathway, C3 TOP2A+ EPCs mainly interact with immune cells such as T-NK cells and myeloid cells through paracrine effects. Finally, we used centrality scores to visualize the role of C3 TOP2A+ EPCs in this pathway ( Figure 5G ).

Transcription factors can act directly on the genome, regulate gene transcription by binding to specific nucleotide sequences upstream of genes, and affect the biological function of cells. Therefore, we used the scenic algorithm to analyze the gene regulatory network of EPCs. Different transcription factors can jointly regulate the expression of certain genes. Therefore, we divided the transcription factors of EPCs into 5 modules based on the connection specificity index: M1, M2, M3, M4 and M5 ( Figure 6A ). As shown in Figures 6B, C , the transcription factors in the M4 and M5 modules may play a major regulatory role in the biological characteristics of C3 TOP2A+ EPCs. Further analysis of the top 5 transcription factors of different EPCs subtypes revealed that the transcription factors of C3 TOP2A+ EPCs were FOXM1, HOXA13, MYBL2, E2F8 and NFYB ( Figures 6D, E ). This is consistent with previous studies showing that E2F8 can maintain tumor cell proliferation and promote tumor cell migration (61, 62). Among them, FOXM1, as an important oxidative stress response regulator and proliferation-related factor, has a significant impact on the progression and outcome of tumors.

To verify the key role of FOXM1 in cervical cancer, we knocked down FOXM1 in CC cell lines ( Figures 7A, B ). CCK-8 showed that after FOXM1 knockdown, the cell viability of cervical cancer cells decreased significantly ( Figure 7C ), and cell proliferation slowed down ( Figure 7D ). The EDU experimental results are consistent with this ( Figure 7E ). The results of the scratch experiment and transwell experiment showed that the reduced expression of FOXM1 significantly inhibited the migration and invasion ability of CC cells ( Figures 7F, G ). Based on the above results, FOXM1 is a key target for inhibiting the progression and invasion of CC and a potential target for future clinical drug development.