Integrated-analysis of single-cell RNA-seq and bulk RNA-seq reveals MMP-mediated expression patterns by distinct tumor microenvironment immune profiles in cervical cancer

doi:10.21203/rs.3.rs-3910841/v1

Download PDF

Article

Integrated-analysis of single-cell RNA-seq and bulk RNA-seq reveals MMP-mediated expression patterns by distinct tumor microenvironment immune profiles in cervical cancer

https://doi.org/10.21203/rs.3.rs-3910841/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

Few studies have analyzed the effect of matrix metalloproteinase (MMP) expression patterns on the tumor microenvironment (TME) during development of cervical cancer (CC).

Methods

We elucidated the landscape and score of MMP expression in CC using single-cell RNA sequencing (scRNA-seq) and RNA sequencing datasets. Further, we aimed the MMPscore to probe the infiltration of immune cells. Further, MMP expression was measured by quantitative Real-Time Polymerase Chain Reaction (qRT-PCR).

Results

We found MMPs were cell-type specific expressed in diverse types of CC cells, regulating the relative pathways of CC progression. Two distinct MMP expression patterns that associated infiltrated tumor microenvironment (TME) were identified. We discovered MMP expression patterns can predict the stage of tumor, subtype, stromal activity in the TME, genetic variation, and patient outcome. Patients with high MMPscore benefited from significantly better treatment and clinical outcomes.

Conclusion

These results indicate high MMPscore in diverse cell types may regulate immune response and improve the survival of patients with CC, which assist in developing more effective immunization strategies.

Biological sciences/Cancer/Cancer genomics

Biological sciences/Cancer/Cancer microenvironment

Biological sciences/Cancer/Gynaecological cancer/Cervical cancer

Biological sciences/Computational biology and bioinformatics/Computational models

Biological sciences/Computational biology and bioinformatics/Genome informatics

Biological sciences/Immunology/Gene regulation in immune cells

Biological sciences/Immunology/Inflammation

Biological sciences/Immunology/Tumour immunology

Matrix metalloproteinase

Cervical cancer

Tumor microenvironment

Extracellular matrix

Single-cell sequencing

Immunotherapy

Cervical cancer (CC) is the most common gynecological malignancy. Risk factors for CC include human papillomavirus (HPV) infection, long-term unhealthy lifestyle maintenance, eating habits, vaginal microenvironment, and regional differences ^1–4. Although the use of HPV vaccines has been promoted in recent years, the global incidence of CC remains high. According to the American Cancer Society, the estimated number of new CC cases in the United States will reach 14,100, and the estimated number of deaths will reach 4,280 in 2022, compared with 14,480 and 4,290, respectively, in 2021 ⁵. In 2020, the diagnosis rate of CC in Europe was 58.2/100,000, with a mortality rate of 26.0/100,000 ⁶. The citywide CC incidence in Guangzhou, China, displayed an annual increase of 2.1% from 2004 to 2018, with the most substantial increase in rural areas, with an average annual percentage change of 6.6% ⁷. In recent years, the proportion of young patients with CC has increased. The incidence of undiagnosed early-stage squamous cell carcinoma has consistently increased in young women aged < 50 years in Japanese metropolitan areas ⁸. Invasive and distant metastases in CC are associated with a substantially low five-year survival rate. Thus, developing an index to assess tumor aggressiveness is crucial to evaluating the prognosis of patients with CC and guiding their clinical treatment.

The tumor microenvironment (TME) is the cellular environment in which tumors or tumor stem cells exist, including tumor cells, adipocytes, fibroblasts, lymphocytes, dendritic cells, cancer-associated fibroblasts, and tumor vasculature, and is widely associated with tumorigenesis. These cells interact with circulatory and lymphatic systems to promote tumorigenesis and progression ^9–11. Single-cell transcriptomics has revealed heterogeneity in tumor- and tumor-derived endothelial cells in CC ¹². Other cancers exhibit heterogeneity in tumor cells and tumor-associated cells. Cancer-associated fibroblasts of different origins contribute to the heterogeneity of tumor cells and exert functional effects on tumors through various mechanisms ^{13, 14}. Tumor-associated neutrophils are associated with the heterogeneity of lung cancer cells ¹⁵. In breast cancer, tumor-associated macrophages exhibit marked heterogeneity and comprise anti-tumor M1-like tumor-associated macrophages (TAM) or pro-tumor M2-like TAM ^{16, 17}.

Two critical stages of tumor development are the degradation of the basement membrane and the invasion of tumor cells into the surrounding tissues. Metastasis of cancer cells is a complex multistep process involving changes in intercellular adhesion, degradation of the extracellular matrix (ECM) and basement membrane, detachment of tumor cells in situ, and extensive infiltration of proteolytic enzymes into lymphatics or blood vessels, which ultimately enhance the ability of tumor cells to invade and metastasize. Matrix metalloproteinases (MMPs) are a series of diverse protein enzymes involved in ECM degradation, mainly by degrading collagen IV and laminin. Various MMPs are produced by tumor and tumor-related cells. During tumor development, MMPs enable tumor cells to cross the matrix membrane barrier and increase their invasiveness, migration, and metastasis, thereby promoting tumor progression. MMP-9 secreted by neutrophils, mast cells, and macrophages degrade the main components of the basement membrane to promote tumor invasion ^18–20. Fibroblasts and tumor cells can secrete MMP-13, MMP-7, and MMP-14. MMP-13 promotes tumor angiogenesis ²¹, MMP-7 degrades heparin binding epidermal growth factor and E-cadherin in the basement membrane ^{22, 23}, and MMP-14 degrades CD-44 and electron-cadherin in the basement membrane ^{24, 25}. Together, these MMPs play a vital role in tumor invasion. MMP-10 is highly expressed in squamous cells and promotes the recruitment of infiltrating cells by remodeling the ECM. MMP-10 can also upregulate the expression of MMP-7, MMP-9, and MMP-13, which are critical for tumor progression ²⁶. MMPs are essential for neovascularization, the inflammatory response, and apoptosis, and play a major role in leukocyte infiltration and tissue inflammation. In the TME, MMP-1 and MMP-2 secreted by fibroblasts promote tumor growth ^{27, 28}. MMP-3 secreted by fibroblasts is an essential mediator of tumor angiogenesis and progression ²⁷. Collectively, MMPs play an important role in regulating the mechanisms of tumor metastasis.

However, the difference in the expression of MMPs secreted by heterogeneous tumor cells and tumor-related cells, as well as the relationship between the expression level of MMPs and the occurrence and progression of CC, remains unclear. Furthermore, its predictive effect on the prognosis of CC also remains unexplored. In this study, we used single-cell RNA sequencing (scRNA-seq) datasets, including patients with CC and healthy individuals, to analyze the MMP landscape at the single cell level. We discovered that MMPs displayed cell-type-specific expression patterns in CC cervix tissues. We analyzed differentially expressed MMPs and their association with the infiltration of immune cells. We quantified the MMP expression of individual tumor cells by constructing the MMPscore to assess the effect of heterogeneity in MMP expression. We believe our study would assist in evaluating whether the MMPscore can effectively predict the prognosis of patients with CC and its significance in guiding the clinical treatment of patients with CC.

2.1 Single-cell sequencing

Single-cell RNA -sequencing (scRNA-seq) datasets of CC cells and healthy donor cervix cells were acquired from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/, ID: GSE168652), The scRNA-seq data were acquired from one example of CC patient and one healthy donor, and RNA-seq was performed by a 10x Genomics Chromium platform. The UMAP algorithm was used to project the data for visualization, and an unsupervised clustering algorithm was used to classify the cells. Differentially expressed genes (DEGs) were identified using the FindMarkers function embedded in Seurat to precisely annotate cell types; we manually curated genetic markers for each cell type. Most markers that distinguished between different cell types were retrieved from the Cell Marker Database (https://www.labome.com/method/Cell-Markers.html), with the following cutoff thresholds: Benjamin-Hochberg's adjusted p value < 0.01 and Fold Change > 1.5. The DEGs were loaded into the cluster profile package for Gene Ontology (GO) enrichment analyses. Pathways with an adjusted p value of < 0.05 were considered significantly enriched. We then performed a gene set enrichment analysis (GSEA) using the gsea function embedded in the cluster profiler package to detect which gene set was significantly enriched in each specific cell cluster. Another set of scRNA-seq datasets from the ArrayExpress database with access to E-MTAB-11948 was employed to demonstrate the correlation between MMP regulatory factors and major functional pathways in CC tissue, the samples are from 3 CC patients and 3 healthy donors.

2.2 CC dataset origin and pretreatment

With the workflow of our study shown in graphical abstract, public gene-expression data and full clinical annotation were searched in GEO and the Cancer Genome Atlas (TCGA) database. Patients without survival information were removed from further evaluation. Cancer genome map from TCGA database (http://cancergenome.nih.gov/), including 11 patients with CC and 264 with normal uterine cervical lesions (n = 261) and normal ectocervical tissue (n = 3), we identified co-expressed genes in the two expression modes of MMPcluster-A and MMPcluster-B, constructed Venn diagrams, and analyzed the copy number variation, somatic mutation data, and survival information. With gene expression summary from GEO database (https://www.ncbi.nlm.nih.gov/geo/,ID: GSE192897), we retrieved RNA-seq data from 16 patients with CC, encompassing 11 patients with CC and 5 with normal uterine cervical lesions (n = 11) and normal ectocervical tissue (n = 5). The clinical information obtained included sex, age, tumor stage, grading, and survival analysis.

2.3 Unsupervised clustering of different MMP expression levels

Three differentially expressed MMPs between CC patients and normal participants were selected by copy number variation (CNV) analysis, and the remaining seven MMPs closely related to CC were selected based on document retrieval. Through heatmaps, Spearman analyses and correlation network diagram analyses, we further confirmed the correlation of 10 MMPs. We performed Spearman’s correlation analysis to obtain the Pearson correlation between 10 MMP phenotypes associated with CC and determined two different expression patterns by unsupervised clustering: MMPcluster-A and MMPcluster-B. The number of clusters and their stability were determined by the consensus clustering algorithm. Principal component analysis was used to verify the existence of these two expression patterns. We used the ConsensuClusterPlus package to per-form the above steps with 1000 times repetitions for guaranteeing the stability of classification ²⁹.

2.4 Enrichment Analysis (GSVA) MMPs modifier mode and functional analysis of different phenotypes

The differences in biological processes between the two expression modes of MMPs were analyzed by the “GSVA” package of R software for GSVA enrichment analysis. In nonparametric and unsupervised methods, GSVA was often utilized to estimate changes in the pathway and biological process activities in samples of expression datasets³⁰. We downloaded the gene set for “c2.cp.kegg.V6.2.Symbols” from the MSigDB database for running the GSVA package (https://www.gsea-msigdb.org/gsea/msigdb/). An adjusted p value less than 0.05 was thought statistically crucial. Functional annotation of MMP-related genes was performed by the clusterProfiler R software package with the cutoff value of FDR less than 0.05.

2.5 Estimation of TME cell infiltration

We used the single-sample gene-set enrich-ment analysis (ssGSEA) algorithm to quantify the relative abun-dance of each cell infiltration in the CC TME. In a study by Charoentong, a gene set labeling each TME-infiltrated immune cell type was obtained, which stored a variety of human immune cell subtypes, including activated dendritic cells, activated CD8 + T cells, natural killer T cells, regulatory T cells, macrophages and so on^{31, 32}. The relative abundance of each TME-infiltrated cell in each sample was represented by the enrichment fraction calculated by ssGSEA analysis.

2.6 Typing of MMP phenotype co-expression gene

First, the R software package “limma” was used to determine the DEGs between CC and normal cervical tissue. From the DEGs, we further searched for genes to express genes in MMPcluster-A and B modifier modes. Finally, unsupervised clustering was employed to reveal the three different MMP phenotypes, which were named MMP gene clusters A–C. Unsupervised clustering analysis was applied to identify the distinct MMP expression patterns and classify patients for further analyses.

2.7 Generation of MMP gene signature

To calculate module scores and the fraction of enrichment for MMP-related gene expression in single cells, AddModuleScores were performed in scRNA-Seq database. A scoring system to evaluate the MMP expression pattern of individual patients with CC was constructed to quantify the MMP expression pattern of a single tumor. The procedures for MMPscore establishment were as follows: first, the two expression modes of MMPcluster-A and -B were processed, and the overlapping genes were extracted. The patients were classified into several groups for deeper analyses by adopting unsupervised clustering methods for studying overlap DEGs. The consensus clustering algorithm was utilized for defining the number of gene clusters as well as their stability. The prognosis of each gene in the feature was analyzed by a one-factor Cox regression model. Further analysis of the genes with significant prognoses was performed. Principal component analysis (PCA) was then made to construct MMP-related gene signatures, with principal components 1 and 2 as signature scores. The advantage of this method was that the scores were concentrated on the set with the largest block of well correlated (or anticorrelated) genes in the set, while down-weighting contributions from genes didn't track with other set members. We then defined the MMPscore, whose method was similar to GGI^{33, 34}:

$$\text{M}\text{M}\text{P}\text{s}\text{c}\text{o}\text{r}\text{e}=\sum _{}^{}\left({PC1}_{i}+{PC2}_{i}\right)$$

where i is the expression of MMP genes.

2.8 Relevance of MMP gene traits to other relevant biological processes

We searched for genes associated with the biological processes of CC by a public database, including immune checkpoints and antigen processing and presentation. The analysis of MMP-related gene traits and TME immune cells further revealed a link between MMP gene traits and related biological pathways.

2.9 Genomic and clinical information on immune checkpoint inhibitors

First, by systematically searching the immune checkpoints to block the gene expression spectrum, two immunotherapy cohorts were included: PD-1 and CTLA4, as representatives of combined immunotherapy. Thus, MMP expression patterns can predict a patient’s response to immune checkpoint blockade therapy.

2.10 Patients and data collection

Patients were selected from CC screening of eligible women who were residents in Shanxi province. Non-typical squamous cells and above with uncertain significance were diagnosed using liquid-based cytology, while the control group consisted of normal cervical cells. Participants with abnormal Pap test results were examined via colposcopy and histopathology. Three cases of CC tissue and three cases of normal cervical tissue were included in the experiment.

The following participants were excluded: (1) pregnant women, (2) patients with a history of hysterectomy, (3) patients with a history of treatment for cervical and vaginal lesions; (4) patients with other malignant tumors; and (5) patients with blood and digestive system disorders. This study was approved by the Ethics Committee of the Second Hospital of Shanxi Medical University with Ethics Number: 2023 YX No.158. We obtained written informed consent from all participants. All patients with CC underwent colposcopy, histopathology, and HPV nucleic acid testing and typing. All normal patients underwent histopathological and HPV nucleic acid testing and typing, and samples were taken from the cervix of patients undergoing total hysterectomy due to uterine fibroids. There was no statistically significant difference in age between the two groups of patients (p > 0.05).

For CC tissue, colposcopy examination was performed by a gynecologist with more than two years of experience. To limit the deviation of routine pathological diagnosis of CC, forceps were used to remove two additional pieces of tissue (approximately 5 mm each). The tissue samples were washed with physiological saline and placed in a solution containing RNA preservation and tissue fixation. The samples used for quantitative Real-Time Polymerase Chain Reaction (qRT-PCR) were stored in a refrigerator (4 ℃) overnight and then transferred to a -20 ℃ refrigerator for storage. Hematoxylin-eosin staining (HE) samples kept at room temperature.

2.11 Hematoxylin-eosin staining(HE)

Samples were fixed using a paraformaldehyde fixing solution overnight at room temperature or 37 ℃ and dehydrated using a Histo-Tek VP1 dehydrator (Sakura, Japan). Tissue samples were cut, embedded, and placed in a freezer until the paraffin was completely set. The paraffin blocks were sectioned using a LEICA RM2235 paraffin microtome (LEICA, Germany) at a thickness of 3–5 µm. Slice harvesting: TKY-TKB slice baking machine (Taikang, Hubei, China) is used to spread, flatten, shape and dry slices. Samples were incubated at 60 ℃ for 2–3 h, then removed and cooled before dewaxing. Samples were passed through several solutions as follows: xylene I 10 mins, xylene II 10 mins, xylene III 10 mins, 100% alcohol I 5 mins, 100% alcohol II 5 mins, 95% alcohol 5 mins, 80% alcohol 5 mins, 70% alcohol 5 mins, distilled water for 5 mins, PBS for 5 mins three times. Samples were dyed using hematoxylin (40 s), rinsed with tap water 3–5 times, differentiated using ethanol hydrochloride (1 s), and rinsed with tap water 3–5 times. Thereafter, samples were dyed with eosin (40 s), rinsed with tap water 3–5 times, and rinsed with 70% alcohol for 10 s, 80% alcohol for 30 s to 1 min, 95% alcohol for 2 mins, 100% alcohol for 5 mins, and 100% alcohol II for 5 mins. The slide was made transparent by xylene three times (5 mins each time), and then the neutral gum was dropped in the center of the tissue on the slide. The cover glass was gradually covered placed over the sample. The mounted slides were then examined and photographed using the 3D Histech Digital Pathology System (Budapest, Hungary).

2.12 Quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNA was isolated from tissue samples using Trizol reagent (Ambion, USA), and total mRNA was reverse transcribed to cDNA using the cDNA Synthesis SuperMix (TransGen Biotech, Beijing). Real-time PCR amplification was performed using a Real-Time PCR system (ROCHE, Lightcycler 96, Switzerland). The following thermal profile was used: segment 1 – one cycle at 95 ℃ for 30 s, segment 2–40 cycles at 95 ℃ for 5 s, followed by 60 ℃ for 30 s. The mRNA expression level was normalized to that of the housekeeping gene GAPDH.

The fold-change in the expression of each target mRNA relative to GAPDH was calculated using the CT (2 − ΔΔCT) method. The experiments were performed at least three times, and the resultant data were statistically analyzed. The primer sequences are listed in Table 1.

Table 1: Primer sequences used in the study.

Genes	Primer nucleotide sequence (5′ to 3′)
MMP-2	F: CCTACACCAAGAACTTCCGTCTG	R: GTGCCAAGGTCAATGTCAGGAG
MMP-3	F: CCTTTCCTGGCATCCCGAAGTG	R: GCCTGGAGAATGTGAGTGGAGTC
MMP-7	F: GAGGATGAACGCTGGACGGATG	R: AGGATCAGAGGAATGTCCCATACCC
MMP-8	F: GGAACGCACTAACTTGACCTACAGG	R: AACACTCCAGAGTTCAAAGGCATCC
MMP-9	F: CTGGTCCTGGTGCTCCTGGTG	R: CTGCCTGTCGGTGAGATTGGTTC
MMP-12	F: TGGACCTGGATCTGGCATTGGAG	R: TCGTGAACAGCAGTGAGGAACAAG
MMP-13	F: GGAGATGAAGACCCCAACCCTAAAC	R: CGGAGACTGGTAATGGCATCAAGG
MMP-19	F: CGGAGACTGGTAATGGCATCAAGG	R: CAAAGGGCAGACACTCGGAACAAG
GAPDH	F: GCTCTCTGCTCCTCCTGTTC	R: ACGACCAAATCCGTTGACTC

2.13 Statistical analysis

The correlation coefficient between TME-infiltrated immune cells and the phenotypic expression of MMPs was calculated using Spearman and distance correlation analysis. Differences between three or more groups were compared using one-way analysis of variance and the Kruskal–Wallis tests ³⁵. The cutoff point for each dataset subgroup was determined using the SurvMiner R software package and was based on the correlation between the MMPscore and patient survival. The MMPscore was bisected using the “SURV-Cutpoint” function, and patients were divided into high and low MMPscore groups based on the log-rank statistic of the largest selection. A survival curve for prognostic analysis was generated using the Kaplan–Meier method, and the log-rank test was used to determine the significance of the difference. The value of the MMPscore in predicting the efficacy of different tumor immunotherapies was evaluated using the Tumor Immune Dysfunction and Exclusion (TIDE) database. The specificity and sensitivity of the MMPscore were evaluated using the receiver operating characteristic (ROC) curve, and the area under the curve (AUC) was quantified using the PROC R software package(https://tide.dfci.harvard.edu/). Copy number variation atlases of the three MMP phenotypes on 23 pairs of chromosomes were mapped using the R package of RCIRCOS. The p-values were analysis from scRNA-Seq datasets of 20,938 cells both bilateral, and a p value of < 0.05 was considered statistically significant. All data processing was performed using R3.6.1 software.

3.1 Single-cell transcription mapping and cell typing of the CC matrix

To understand the cellular diversity and molecular features of cervical tissues in patients with CC, we analyzed scRNA-seq datasets of 20,938 cells from patients with CC and healthy control individuals. Of these cells, 10,395 were obtained from the tumor samples that comprised seven cell types, including Myeloid cells, lymphocytes, endothelial cells, endometrial stromal cells, fibroblasts, smooth muscle cells, and CC cells, which were compared with those of the normal group and annotated using classical marker genes (Fig. 1a, and S1a). The heatmap showed marker genes (Figure S1b and S1c). To understand the proportion of MMP molecules in each cell type, we used the bar graph to reveal the percentage of MMPs gene expression in each cell type. We found that MMP-1, MMP-7, MMP-13 and MMP-14 had a high proportion of diverse cell types (Fig. 1b). We also made conclusion that MMP-3, MMP-7, MMP-13 and MMP-14 were up-regulated in CC cells in individuals and MMP-9, MMP-14 and MMP-19 were up-regulated in Myeloid cells by employing the heatmap(Fig. 1c). To elucidate the mechanism of MMP regulation in diverse cell types in CC, we identified DEGs in individuals between MMPs high-expression and low-expression groups, finally it was found MMP-1, MMP-12, MMP-13 and MMP-7 were upregulated, but MMP-14, MMP-19, MMP-2, and MMP-21 were downregulated in diverse cell types (Fig. 1d), suggesting biochemically redundant members of the MMP family may have intricate interplay in tumor progression. To further explain the relationship between the expression of five MMPs in diverse cell types, we performed the re-clustering and Uniform Manifold Approximation and Projection (UMAP) dimensionality reduction profiles of MMPs in the CC and control groups. We found that the expressions of MMP-14, MMP-19 and MMP-2 were upregulated in normal individuals compared to CC (Fig. 1f–h), during which fibroblasts and myeloid cells displayed abundant MMP expression, indicating that the MMPs gene set could participate in inhibits tumor progression of CC. To identify whether MMP molecules can be linked to an favorable prognosis in CC, we then put MMPs as a gene module to calculate the score using AddModuleScore methods. The violin plot illustrated the fibroblasts, CC, and Mye cell had a higher score than the other cell types, consistently showing that the MMPscore was associated with immune responses in CC (Fig. 1e).

To investigate the association between the classical CC pathways and MMPscore more accurately, we analyzed another scRNA-seq datasets from 3 normal and 3 CC (Fig. 1i, Fig. 1j), which encompassed 50,014 cells from patients with CC and healthy control individuals that consisted of seven cell types, including macrophage, lymphocytes, Neutrophils, endothelial cells, fibroblasts, smooth muscle cells, and epithelial cells (CC cells) (Fig. 1l). The CC patient group had more abundant epithelial cells, macrophage and lymphocytes (Figure S1d). To understand the proportion of MMP molecules in each cell type, the bar graph showed the percentage of MMP-1, MMP-2, MMP-7, MMP-13 and MMP-14 had a high proportion of diverse cell types (Fig. 1n). By identifying DEGs between MMP high-expression and low-expression groups, MMP-1, MMP-3, MMP-12, MMP-13 and MMP-7 were upregulated, but MMP-2, MMP-14 and MMP-19 were downregulated in diverse cell types (Figure S1e-g). We also put MMPs as a gene module to calculate the score using AddModuleScore methods and to define scores separately. Then we found that fibroblasts, CC, and Mye cell also had a higher score, suggesting that MMP-related genes played function more in these cell types (Fig. 1k). A higher score for CC patients than for normal may indicate a potential association of MMP-related genes and CC. To probe the association between the MMP-regulation and progression of CC, we used the functional enrichment analyses based on the GSEA database, and we performed the correlation between MMPscore by AUCell and classical pathways in CC to explore the influence of MMP regulators on CC pathways. Notably, we found MMPscore was upregulating TGF-β pathway (R > 0.6) and epithelial cell proliferation (R > 0.8), but downregulating PD-1 (R < -0.8), HPV copy number (R < -0.4), Wnt pathway (R < -0.8) and CC cell proliferation (R < -0.2) in CC patients (Fig. 1m). To explore whether MMPs were associated with Wnt pathway in regulation of cell proliferation, with defense to virus by host and epithelial-mesenchymal transition (EMT). A correlation analysis combined with the MetaCell algorithm ( K = 30 ) revealed that MMP-1 (R = -0.79), MMP-7 (R = -0.69) and MMP-13 (R = -0.73)were negatively associated with EMT in CC. MMP-1 (R = -0.8), and MMP-7 (R = -0.81) and MMP-13 (R = -0.76)were negatively associated with Wnt pathway in regulation of cell proliferation. Scatter plot showed that MMP-7 (R= -0.57) and MMP-13 (R = -0.55) were negatively associated with positive regulation of defense response to virus by host in CC (Fig. 1m and Figure S1h-k). These results indicated that cell-type-specific up-regulation of MMP molecules was in normal cervical epithelial cell proliferation and inhibition to HPV infection, Wnt signaling, CC cell proliferation, EMT and PD1 in tumor proteins which controlled CC progression.

3.2 MMP genetic variation and expression landscape in CC

In this study, we first summarized the prevalence of MMP CNV and somatic mutations in CC. Investigating CNV alteration frequency showed prevalent CNV alteration in 21 regulators of three MMPs, mostly focusing on copy number amplification (Fig. 2a2). The chromosomal locations of MMP CNV alterations were shown in Fig. 2a1. By investigating MMP mRNA expression levels in normal and CC samples to determine whether the above genetic variations influenced MMP expression in patients with CC, we found that CNV change may be the main factor leading to abnormal MMP expression. Compared with normal cervical tissues, MMP-9 with amplified CNV demonstrated markedly higher expression in CC tissues (Figure S2e). There were significant differences in the expression of ten MMPs between patients with CC and normal controls (Fig. 2b). The correlation network diagram in Fig. 2c depicted MMP interactions and further confirmed the ubiquitous correlation of the ten MMPs. We then examined the Pearson correlation between ten MMPs associated with CC by Spearman’s correlation analyses, which observed a positive correlation among MMPs (Fig. 2d). In TCGA cohort, we used the R package “ConsensusClusterPlus” to categorize a series of patients with different MMP expression patterns according to the expression of the ten MMPs. With final identification of two distinct modification patterns by unsupervised clustering (Fig. 2e), we named these patterns MMPclusterA and -B that were supported by the principal component analysis results (Fig. 2f, Figure S2f-h, Figure S3a-b). Prognostic analysis of the two MMP expression subtypes revealed a prominent survival advantage within the MMPcluster-B expression pattern (p < 0.05) (Fig. 2g). We further explored the different clinicopathological features and prognoses between the two groups (Fig. 2h). The assessment of these data showed ten highly expressed MMPs in MMPcluster-A (Fig. 2i). In addition, we also examined the different expression of four MMP types, and a significant correlation with overall survival prediction was seen in patients with CC (Figure S2a-d).

3.3 Immune cell infiltration characteristics in the TME under different MMP expression patterns

The biological behaviors of various MMP expression patterns were analyzed by Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment. MMPcluster-A was significantly enriched in the matrix and carcinogenic activation pathways, like ECM receptor interactions, small cell lung cancer signaling pathways, cell adhesion, and MAPK signaling pathways, exhibiting significant aggregation during cell matrix catabolism. MMPcluster-B was mainly enriched in gene mutations, encompassing activation of homologous recombinant signaling pathways and RNA degradation signaling pathways (Fig. 2k). Furthermore, MMP expression patterns greatly differed in terms of cell composition, molecular function, and biological processes. MMPcluster-A presented a generally enriched result, especially in collagen catabolism, metabolism, ECM catabolism, and constituent components, unlike MMPcluster-B (Fig. 2l). Subsequent analyses indicated a higher degree of immune cell infiltration in the TME of MMPcluster-A, including activated T cells, CD56 + NK cells, MDSC, and plasmacytoid dendritic cells, while infiltration of CD56- NK cells, dendritic cells, and monocytes exhibited no clear difference between MMPcluster-A and -B (Fig. 2j). Combination of KEGG, GO, and TME immune cell infiltration analysis showed that MMP expression patterns had obvious immune cell infiltration characteristics. Despite the significant immune cell infiltration, MMPcluster-A did not show a matching survival advantage (Fig. 2g).

3.4 Generation of MMP gene signatures and functional annotations

To further investigate the potential biological behavior of each MMP expression pattern, we identified two expression patterns (MMPcluster-A and MMPcluster-B) by unsupervised clustering. To prepare for the establishment of the MMPscore and visualize the heatmap of DEGs between MMPclusterA and -B, we obtained a total of 350 genes using DEG analyses (Fig. 3a). Furthermore, through unsupervised clustering, DEGs were once again divided into three gene clusters: MMP gene cluster A-C (Fig. 3b), representing the genetic differences in MMPs. To quantify the MMP landscape and facilitate the identification of key genes, PCA was used to compute the aggregate score of feature genes from three different MMP genome phenotypes, respectively. We obtained the sum of scores and defined them as the MMPscore. All TCGA patients were stratified into two groups with high or low MMPscore. As indicated from the prognostic analyses, the prognosis of the high score was better than that of the low score, 27 of 253 patients with CC clustered in gene cluster B, which was shown to be associated with a better prognosis. In contrast, 94 patients in gene cluster A exhibited poorer prognosis. An intermediate prognosis was observed in gene cluster C with 132 clustered patients (Fig. 3c). In addition, we found that younger patients and patients with advanced diseases were significantly associated with lower MMPscore, implying that these patients were characterized by a poorer clinical outcome in the MMPcluster-A modification pattern (p < 0.05) (Fig. 3d), which suggested the two distinct MMP expression patterns in CC. We observed that MMP gene clusters A and B were consistent with the above expression patterns, whereas tumors in MMP gene cluster C exhibited poorer differentiation. Patients with clinically advanced diseases were characteristic of MMP gene cluster A, consistent with the above results (Fig. 3e). By performing ontology enrichment analyses of different MMP gene clusters, we found that they were significantly enriched in ECM, collagen, and extracellular structural tissues. Considering the individual heterogeneity of MMP expression (Fig. 3f), we quantified MMP expression in individual tumor cells by calculating the MMPscore. Compared with the other clusters, MMPcluster-A showed a significantly lower MMPscore, while MMPcluster-B displayed a high median score (p < 0.05) (Figure S3c). The Kruskal-Wallis test implied crucial differences in MMPscore between MMP gene clusters (p < 0.05). Gene cluster A had the lowest median score, while gene cluster B had the highest median score (Fig. 3d). Consistent with the above analysis, this suggested that a high MMPscore may be strongly associated with immune activation, whereas a low MMPscore may be linked to stromal activation. Next, we sought to further determine the importance of the MMPscore in predicting patient prognosis. Patients were divided into low and high MMPscore groups. Patients with high MMPscore showed a significant survival benefit (p < 0.001), approximately twice that of patients with low MMPscore (Fig. 3g). Besides, we found extremely different immune cell infiltration in the three gene clusters (Fig. 3h). Further single-sample GSEA revealed that different MMPscore were significantly associated with high and low levels of immune infiltration in tumor tissues (Fig. 3i, Fig. 4a). The high MMPscore group was dominated by immune cell activation that was closely correlated with prognosis (Fig. 3i). The above results strongly displayed that a high MMPscore was associated with increased immune activation. The MMPscore allowed for a better assessment of the MMP expression patterns of individual tumors and further assessed the TME infiltration characteristics of tumors to distinguish between true and false TME immune infiltration.

3.5 Role of MMP expression patterns in immunotherapy

To better characterize the MMP immune profile and test the correlation between immune cells and MMPscore, we examined the specific correlation between each TME-infiltrating cell type and high and low MMP expression, which showed a tight correlation (Fig. 3j). Our study implied that TME immune cell infiltration was significantly increased in tumors with high MMPscore, which meant a more significant positive correlation with T cell subsets, mast cells and B cells (p < 0.05) (Fig. 3k-m, Figure S3e-m). The investigation about the correlation between MMPscore and adhesion molecules, as well as between HLA molecules and interleukins, showed that MMPscore correlated significantly with immune checkpoints (Fig. 4a-c) with the most significant correlation of CD44 and TNFSF9 (Fig. 4d-f, Figure S3n-p). MMPscore can predict the strength and weakness of immune function with immunotherapy performed and evaluated for immune checkpoints. In addition, we found that different types of HLA also correlated with MMPscore (Fig. 4b), with HLA-E and HLA-C the most prominent (Fig. 4e, Figure S3p). For immune regulation, we found a correlation between MMPscore and interleukins (Fig. 4c), with TLSP and IL-33 of more value (Fig. 4f, Figure S3o). MMPscore also had predictive value in evaluating immune escape. The mRNA transcriptome showed differences with various MMPscore and correlated importantly with immune-related biological pathways (Fig. 4g). The expression of MMP may be involved in immune pathways and activation or inhibition. CD47 increased immune escape, and a high MMPscore was associated with low immune escape in a study examining the effect of MMPscore on immune checkpoint blockade therapy (Figure S4a). These results showed that tumors with high MMP expression showed significantly correlated immune activation pathways.

3.6 Characterization of MMP expression in TCGA molecular subtypes and tumor somatic cell mutations

TCGA constructed a comprehensive molecular landscape for CC and investigated the predictive ability of MMPscore in CC prognosis in patients with different tumor loads, revealing better survival prognosis in the high MMPscore group than in the low MMPscore group (Fig. 4h), which was not significantly related to high or low tumor loads (Fig. 4i). We counted 179 MMP mutations in 208 CC samples with a mutation frequency of 86.06%, of which TTN showed the highest mutation frequency followed by PIK3CA, and other genes also displayed varying degrees of mutation (Fig. 4j). The somatic mutation distribution differed between high and low MMPscore in the cohort (Fig. 4j, k): the low MMPscore group exhibited a wider range of tumor mutations, compared with the high MMPscore group. The rates of the 2nd and 14th most significant mutated genes were 27% and 7%, and 9% and 28%, respectively. The above analyses indicated that a high MMPscore was associated with long-term survival and durable clinical benefit. The combination of MMPscore and tumor mutation load could more accurately determine prognosis in different patient cohorts, and the predictive advantage of ROC curve assessment was reflected in the survival prognosis of patients with CC (Figure S3d).

3.7 Clinical features of MMP expression patterns

Consistent with the above findings (Fig. 5a), gene clusters B and C both showed high MMP expression and better survival outcomes, while gene cluster A, containing both high and low MMP expression cases, exhibited poorer survival outcomes, approximately 50% of patients with a prognosis of death. The above results again suggested that MMP expression profiles played a non-negligible role in shaping different TMEs. Next, we used MMPscore to systematically evaluate CC in terms of its clinical characteristics, including age, weight, smoking, clinical stage, and prognosis (Fig. 5b-k, Figure S4c-l). Once again, the high MMPscore was dominated by survival outcomes in the Fustat survival prognosis (Fig. 5b, d, h). In addition, high MMPscore also correlated significantly with TNM stage, especially in patients with Nx, T2, and G3 stages of better survival prognosis (Fig. 5e-g, 5i-k), which showed better assessment in tumor infiltration and metastasis. We also found that the MMPscore could indicate sensitivity to more than ten drugs, including nilotinib and Adriamycin. It can be inferred that high MMPscore can make drug response more sensitive (p < 0.05) (Fig. 5l-o, Figure S2i-p).

3.8 MMP expression is generally increased in CC tissue

By verifying the expression of the MMP gene set in CC tissues using immunohistochemistry, we found that MMP expression in normal cervical tissue showed low or no expression, demonstrating the effectiveness of the MMPscore (Fig. 6k-n). To confirm whether the expression of the MMP gene set was generalized at the molecular level in tissues of patients with CC, we conducted HE and qRT-PCR experiments on cervical tissues from three cases of CC and three healthy individuals. As listed in Table s1, the six subjects enrolled in this study were divided into the CC group (three patients with CC receiving surgery or other treatments) and the non-CC group (three subjects receiving abdominal hysterectomy for uterine leiomyoma or prolapse of uterus). The age distribution, the parity, and times of pregnancy showed no significant group differences. HE staining revealed cell polarity disorder, increased mitotic Figures, and abnormal cells breaking through the basal layer in CC tissues, which was consistent with the diagnosis of CC. HE staining of cervical tissue in the control group conformed to normal cervical pathological characteristics (Fig. 6a-b). We observed that the expression of MMP-2, MMP-3, MMP-7, MMP-9, MMP-12, MMP-13, MMP-14, and MMP-19 genes in CC tissue was generally higher than that in normal cervical tissue (Fig. 6c-j). qRT-PCR analysis displayed that the mRNA expression of MMP-2, MMP-3, MMP-7, and MMP-9 in the cervical tissue of patients with cancer was significantly greater than that in the normal group (p = 0.005, p < 0.001, p < 0.001, and p = 0.038, respectively). mRNA expression levels of MMP-12, MMP-13, MMP-14, and MMP-19 in the cervical tissue of patients with cancer were greater than that observed in the normal group (p = 0.045, p = 0.025, p = 0.003, and p = 0.025, respectively).

MMPs regulate cancer progression by regulating angiogenesis, invasion, and immune escape, but their potential for CC diagnosis and prognosis has not been studied considering the overall MMPs expression pattern ³⁶. MMPs can bind to the cancer cell surface, act on all stages of cancer, and are expressed in early tumor cells to facilitate ECM remodeling and release membrane-bound growth factors that prepare the microenvironment for tumorigenesis ^{37, 38}. MMP-9 plays a role in a wide variety of cancers due to its effects on immune cell infiltration, and it has shown potential to determine prognosis ^39–42. In addition, MMP-1 and MMP-2 have also been associated with CC prognosis ³⁹. MMP-7 can predict a more aggressive colon cancer phenotype and is inversely correlated with patient survival ⁴³. MMP-11 can be considered a potential tumor marker and therapeutic target for advanced prostate cancer ⁴⁴. The mechanism underlying MMP effects in CC still requires further study. Recently, m6A methylation has shown a strong evaluation ability in terms of clinical pathological characteristics, prognosis, and immune phenotype, which revealed MMP expression characteristics in the CC TME in this study ⁴⁴.

Accumulating evidence indicates that the MMP expression profile plays an indispensable role in inflammation ⁴⁵, immunity, and inhibiting tumor progression, especially in the development and progression of digestive tract tumors, but there is no related study concerning CC. Furthermore, most current studies focus on a single TME cell type or a single protease, and the overall TME infiltration characteristics mediated by the combined effects of multiple MMPs have not been comprehensively recognized. We used the scRNA-seq dataset composed of CC patients and healthy individuals to find that MMP exhibited a cell type specific expression pattern in cervical tissue of cancer at the single cell level. This study revealed two MMP expression patterns based on 10 MMPs, which have shown distinctly different TME cell infiltration characteristics. MMPcluster-A is a type of immune rejection characterized by innate immune cell infiltration and stromal activation, as well as an immunoinflammatory phenotype characterized by adaptive immune cell infiltration and immune activation. In contrast, MMPcluster-B was activated by mutations. The immune rejection phenotype comprises non-inflamed tumors, while immune inflammatory phenotypes, referred to as hot tumors, manifest as massive immune cell infiltration in the TME ^{46, 48}. Although the immune rejection phenotype included a large number of immune cells, they remained in the stroma surrounding the tumor cell nests, rather than penetrating their parenchyma. The stroma can be confined to the tumor envelope or penetrate the tumor itself, making immune cells appear truly inside the tumor ^{49, 51}. More importantly, we found that MMPcluster-A exhibited a distinct stromal activation status, combined with TME cell infiltration features in each cluster. Patient prognosis opposed expectations; therefore, we speculated that stromal activation in MMPcluster-A inhibited the antitumor effect of immune cells. Significant prognostic differences existed between the two clusters, confirming the reliability of our immunophenotype classification for different MMP expression patterns. Therefore, fully exploring the TME cellular infiltration characterization induced by different MMP expression patterns demonstrated that MMPcluster-A could further lead to poor prognosis through the function of suppressed immune cells. In addition, comparing MMP genetic and expression alterations between CC tissues and normal tissues showed a certain heterogeneity, indicating that MMP expression imbalance may play an important role in CC occurrence and progression. Our seminal exploration of the overall MMP expression pattern role in TME infiltration in CC will contribute to a deeper understanding of the mechanism of TME antitumor immune response and a more effective strategy for guiding immunotherapy.

In this study, MMP gene analysis identified three genomic subtypes that were also significantly associated with matrix activation and immune response, similar to the clustering results for MMP expression. This demonstrates that MMP expression has important implications for shaping different TMEs. Therefore, a comprehensive evaluation of MMP expression patterns will enhance our understanding of the characterization of TME cellular infiltration. However, previous analysis was mainly based on patient population and could not accurately predict the expression pattern in individual patients. Considering the individual heterogeneity of MMP expression, its pattern must be urgently quantified in single tumors. In this study, this deficiency was well compensated by constructing an MMP scoring system, evaluating the MMP expression pattern in CC patients, and visualizing property changes in individual patients. The expression pattern dominated by the MMPcluster-B expression signature exhibited a high MMPscore, suggesting that the MMPscore is a reliable and powerful tool to comprehensively assess MMP expression patterns in individual tumors that can be used to further determine the TME infiltration pattern, also known as the tumor immunophenotype. More importantly, the MMPscore showed good assessment ability in terms of patient clinical characteristics, including tumor differentiation level, mutation burden, pathological stage, body weight, age, and clinical prognosis, and it could guide clinical treatment. Comprehensive analysis showed that MMPscore was an effective indicator of biological prognosis in CC. Our MMPscore has shown excellent predictive power in CC precision immunotherapy utilizing the immune escape feature.

The heterogeneity of solid tumors and the TME has been well mapped. In this study, we elucidated the heterogeneity and tumor infiltration pattern in CC through different MMP expression patterns. Some studies have shown that the CC TME is heterogeneous ⁵², and we further analyzed its biological and transcriptomic heterogeneity at the single-cell level. Different CC cell markers identified tumor cell heterogeneity, as well as apparent heterogeneity in MMP expression, compared with tumor stromal cells, which suggested that different transcriptomes of individual cells might reflect their tumor biological characteristics, further demonstrating the possibility that MMPs are involved in and influence tumor progression. By using scRNA-seq, we identified all cell types in CC, including those of neutrophils, T cells, smooth muscle cells, trophoblast progenitors, trophoblast stem cells, megakaryocytic progenitors, endothelial cells, myoepithelial cells, Purkinje cells, myeloid cells, and fibroblasts, as well as 12 cell types of unknown significance. Further dimension reduction analysis showed that the proportion of each cell subtype differed significantly between cancer tissues and normal tissue stroma. GSEA of the marker genes of each cancer cell cluster revealed some association with the MMP family, further validating the cell identity and biological pathways described above and elucidating specific gene expression signatures in CC cell types while comparing the different cell infiltration in tumor tissues and showing their MMP enrichment. Single-cell transcriptomics explain the potential link between tumor cells and MMPs by identifying rare cell subpopulations. In addition, the effectiveness of MMPscore was also well validated in CC tissue. Therefore, we anticipate that these findings will provide important clues for the developing applications of the MMP family in CC.

Unfortunately, our data did not show a significant correlation between MMPscore and tumor mutational burden. However, it is worth noting that patients in high MMPscore group showed better prognosis regardless of high or low tumor burden, and then it was speculated that survival analysis by MMPscore combined with tumor mutational burden could yield more nuanced prognostic results. Our results showed that MMP expression plays a non-negligible role in shaping different stromal and immune TMEs, implying that it may influence the therapeutic efficacy of immune checkpoint blockade. In addition, MMPscore integration with various biomarkers, including PD-1 expression and stromal and immune TME, may lead to a more effective CC immunotherapy strategy, which should be verified by further experimental results. Moreover, the correlation between MMPscore and tumor stage, degree of invasion, and prognosis analysis did not show significant effect in all grades. The correlation between MMPscore and patient weight and age factors and its survival predictive effect were also not obvious. With the power of single-cell transcriptomics to identify cellular subsets and interpret links between tumor cells and endothelial cells, scRNA-seq is inherently limited to transcript-level measurements, so the functional implications of each population require further investigation.

Scholars performing pan-cancer analysis showed that MMPs had prognostic value only in clear cell renal cancer ³⁶. This study elucidates the role of MMPs in cancer by developing an MMPscore scoring system, which may serve as an independent marker for predicting patient survival prognosis and provide new insights into CC immunotherapy. These new ideas may target MMP-related genes, reverse unfavorable TME cell infiltration characterization, and help develop novel drug combination strategies or novel immunotherapeutic agents in the future. We provide new ideas for improving patient clinical responses to immunotherapy, identifying distinct tumor immune phenotypes, and promoting personalized CC immunotherapy in the future. In conclusion, the MMPscore can be used in clinical practice to comprehensively evaluate the MMP expression pattern of individual patients and their corresponding TME cell infiltration characteristics, further determine the tumor immune phenotype, and guide more effective clinical practice.

In summary, this work demonstrates a broad regulatory mechanism of the CC TME by the MMP expression landscape. Differences in MMP expression patterns are a non-negligible factor contributing to the heterogeneity and complexity of individual TMEs. The MMPscore exhibited a strong predictive function in CC patient survival analysis that could provide guidance for clinical work-up. Single cell transcriptomics used to investigate the intratumoral heterogeneity of CC at the cellular level validated the strong correlation between MMPs and tumorigenesis and revealed the biological nature of tumorigenesis. Comprehensive assessment of MMP expression patterns in individual tumors will enhance our understanding of the characteristics of cellular infiltration in the TME. The correlation between MMPscore and immune checkpoints and immune cells may provide strategies and directions for subsequent immunotherapy research.

Ethics approval and consent to participate：For studies involving humans, the study was conducted in accordance with the Declaration of Helsinki, and approved by the Ethics Committee of the Second Hospital of Shanxi Medical University with Ethics Number: 2023 YX No.158 and received written informed consent from the participants.

Consent for publication: Informed consent was obtained from all subjects involved in the study. Written informed consent has been obtained from the patients to publish this paper.

Availability of data and materials: The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

Competing interests: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Funding: Please add: This study was funded by the project of Four "Batches" Innovation Project of Invigorating Medical through Science and Technology of Shanxi (2022XM23), the Shanxi Provincial Key Research and Development Project(201803D31111) and The Natural Science Foundation of Shanxi Province, China (202103021224394).

Authors' contributions: YP and ZW were responsible for project design. The first draft of the manuscript and the main experiments were performed by HW, XL, SZ and WG, and these authors contributed equally to this work and should be considered co‑first authors. ZW, HW, XL, SZ, WG, Zhao W, LS, ZZ and SZ: conceptualization, methodology, software, data curation, formal analysis, and validation. ZW, HW, XL: data analyzation, interpretation and visualization. YP and ZW: manuscript review and edit. YP and ZW contributed equally to the correspondence work. YP and ZW: Study coordination, guarantors, writing review. All authors contributed to the article and approved the submitted version.

Acknowledgments: We acknowledge the GEO and TCGA databases for providing their platforms and contributors to upload their meaningful datasets. We also thank the project of Four "Batches" Innovation Project of Invigorating Medical through Science and Technology of Shanxi (2022XM23), the Shanxi Provincial Key Research and Development Project(201803D31111) and The Natural Science Foundation of Shanxi Province, China (202103021224394).

Wang, Z. et al. Lower dietary mineral intake is significantly associated with cervical cancer risk in a population-based cross-sectional study. J. Cancer 12, 111–123, doi:10.7150/jca.39806 (2021).
Bosch, F. X., Lorincz, A., Muñoz, N., Meijer, C. J. & Shah, K. V. The causal relation between human papillomavirus and cervical cancer. J. Clin. Pathol. 55, 244–265, doi:10.1136/jcp.55.4.244 (2002).
Shields, T. S. et al. A case-control study of risk factors for invasive cervical cancer among U.S. women exposed to oncogenic types of human papillomavirus. Cancer Epidemiol. Biomarkers Prev. 13, 1574–1582 (2004).
Kjellberg, L. et al. Smoking, diet, pregnancy and oral contraceptive use as risk factors for cervical intra-epithelial neoplasia in relation to human papillomavirus infection. Br. J. Cancer 82, 1332–1338, doi:10.1054/bjoc.1999.1100 (2000).
Siegel, R. L., Miller, K. D., Fuchs, H. E. & Jemal, A. Cancer statistics, 2022. CA Cancer J. Clin. 72, 7–33, doi:10.3322/caac.21708 (2022).
Dyba, T. et al. The European cancer burden in 2020: Incidence and mortality estimates for 40 countries and 25 major cancers. Eur. J. Cancer 157, 308–347, doi:10.1016/j.ejca.2021.07.039 (2021).
Li, K., Xu, H., Wang, S., Qin, P. & Liang, B. Disparities in the increases of cervical cancer incidence rates: observations from a city-wide population-based study. BMC Cancer 22, 419, doi:10.1186/s12885-022-09531-2 (2022).
Noguchi, T. et al. Recent Increasing Incidence of Early-Stage Cervical Cancers of the Squamous Cell Carcinoma Subtype among Young Women. Int J Environ Res Public Health 17, doi:10.3390/ijerph17207401 (2020).
Chen, D., Zhang, X., Li, Z. & Zhu, B. Metabolic regulatory crosstalk between tumor microenvironment and tumor-associated macrophages. Theranostics 11, 1016–1030, doi:10.7150/thno.51777 (2021).
Mao, X. et al. Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: new findings and future perspectives. Mol. Cancer 20, 131, doi:10.1186/s12943-021-01428-1 (2021).
Zhou, J., Zhang, S. & Guo, C. Crosstalk between macrophages and natural killer cells in the tumor microenvironment. Int. Immunopharmacol. 101, 108374, doi:10.1016/j.intimp.2021.108374 (2021).
Li, C., Guo, L., Li, S. & Hua, K. Single-cell transcriptomics reveals the landscape of intra-tumoral heterogeneity and transcriptional activities of ECs in CC. Mol Ther Nucleic Acids 24, 682–694, doi:10.1016/j.omtn.2021.03.017 (2021).
Bu, L. et al. Biological heterogeneity and versatility of cancer-associated fibroblasts in the tumor microenvironment. Oncogene 38, 4887–4901, doi:10.1038/s41388-019-0765-y (2019).
Patel, A. K. & Singh, S. Cancer associated fibroblasts: phenotypic and functional heterogeneity. Front Biosci (Landmark Ed) 25, 961–978, doi:10.2741/4843 (2020).
Wu, F. et al. Single-cell profiling of tumor heterogeneity and the microenvironment in advanced non-small cell lung cancer. Nat Commun 12, 2540, doi:10.1038/s41467-021-22801-0 (2021).
Ahirwar, D. K. et al. Slit2 Inhibits Breast Cancer Metastasis by Activating M1-Like Phagocytic and Antifibrotic Macrophages. Cancer Res. 81, 5255–5267, doi:10.1158/0008-5472.Can-20-3909 (2021).
Goswami, K. K., Bose, A. & Baral, R. Macrophages in tumor: An inflammatory perspective. Clin. Immunol. 232, 108875, doi:10.1016/j.clim.2021.108875 (2021).
O'Grady, A. et al. Differential expression of matrix metalloproteinase (MMP)-2, MMP-9 and tissue inhibitor of metalloproteinase (TIMP)-1 and TIMP-2 in non-melanoma skin cancer: implications for tumour progression. Histopathology 51, 793–804, doi:10.1111/j.1365-2559.2007.02885.x (2007).
Kuivanen, T., Jeskanen, L., Kyllönen, L., Isaka, K. & Saarialho-Kere, U. Matrix metalloproteinase-26 is present more frequently in squamous cell carcinomas of immunosuppressed compared with immunocompetent patients. J. Cutan. Pathol. 36, 929–936, doi:10.1111/j.1600-0560.2009.01188.x (2009).
Hernández-Pérez, M., El-hajahmad, M., Massaro, J. & Mahalingam, M. Expression of gelatinases (MMP-2, MMP-9) and gelatinase activator (MMP-14) in actinic keratosis and in in situ and invasive squamous cell carcinoma. Am. J. Dermatopathol. 34, 723–728, doi:10.1097/DAD.0b013e31824b1ddf (2012).
Lederle, W. et al. MMP13 as a stromal mediator in controlling persistent angiogenesis in skin carcinoma. Carcinogenesis 31, 1175–1184, doi:10.1093/carcin/bgp248 (2010).
Kivisaari, A. K. et al. Matrix metalloproteinase-7 activates heparin-binding epidermal growth factor-like growth factor in cutaneous squamous cell carcinoma. Br. J. Dermatol. 163, 726–735, doi:10.1111/j.1365-2133.2010.09924.x (2010).
Chan, L. P. et al. Growth regulated oncogene-α contribute to EMT/MMPs pathway by binding its receptors in head and neck squamous cell carcinoma. Life Sci. 306, 120791, doi:10.1016/j.lfs.2022.120791 (2022).
Chuang, H. C. et al. Active matrix metalloproteinase-7 is associated with invasion in buccal squamous cell carcinoma. Mod. Pathol. 21, 1444–1450, doi:10.1038/modpathol.2008.99 (2008).
Roh, M. R. et al. Differential expression patterns of MMPs and their role in the invasion of epithelial premalignant tumors and invasive cutaneous squamous cell carcinoma. Exp. Mol. Pathol. 92, 236–242, doi:10.1016/j.yexmp.2012.01.003 (2012).
Pittayapruek, P., Meephansan, J., Prapapan, O., Komine, M. & Ohtsuki, M. Role of Matrix Metalloproteinases in Photoaging and Photocarcinogenesis. Int J Mol Sci 17, doi:10.3390/ijms17060868 (2016).
Vosseler, S. et al. Distinct progression-associated expression of tumor and stromal MMPs in HaCaT skin SCCs correlates with onset of invasion. Int. J. Cancer 125, 2296–2306, doi:10.1002/ijc.24589 (2009).
Sbardella, D. et al. Human matrix metalloproteinases: an ubiquitarian class of enzymes involved in several pathological processes. Mol. Aspects Med. 33, 119–208, doi:10.1016/j.mam.2011.10.015 (2012).
Wilkerson, M. D. & Hayes, D. N. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572–1573, doi:10.1093/bioinformatics/btq170 (2010).
Hänzelmann, S., Castelo, R. & Guinney, J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatics 14, 7, doi:10.1186/1471-2105-14-7 (2013).
Charoentong, P. et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep. 18, 248–262, doi:10.1016/j.celrep.2016.12.019 (2017).
Barbie, D. A. et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462, 108–112, doi:10.1038/nature08460 (2009).
Sotiriou, C. et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J. Natl. Cancer Inst. 98, 262–272, doi:10.1093/jnci/djj052 (2006).
Zeng, D. et al. Tumor Microenvironment Characterization in Gastric Cancer Identifies Prognostic and Immunotherapeutically Relevant Gene Signatures. Cancer Immunol Res 7, 737–750, doi:10.1158/2326-6066.Cir-18-0436 (2019).
Hazra, A. & Gogtay, N. Biostatistics Series Module 3: Comparing Groups: Numerical Variables. Indian J. Dermatol. 61, 251–260, doi:10.4103/0019-5154.182416 (2016).
Gobin, E. et al. A pan-cancer perspective of matrix metalloproteases (MMP) gene expression profile and their diagnostic/prognostic potential. BMC Cancer 19, 581, doi:10.1186/s12885-019-5768-0 (2019).
Abdel-Hamid, N. M. & Abass, S. A. Matrix metalloproteinase contribution in management of cancer proliferation, metastasis and drug targeting. Mol. Biol. Rep. 48, 6525–6538, doi:10.1007/s11033-021-06635-z (2021).
Mondal, S., Adhikari, N., Banerjee, S., Amin, S. A. & Jha, T. Matrix metalloproteinase-9 (MMP-9) and its inhibitors in cancer: A minireview. Eur. J. Med. Chem. 194, 112260, doi:10.1016/j.ejmech.2020.112260 (2020).
Zeng, Y., Gao, M., Lin, D., Du, G. & Cai, Y. Prognostic and Immunological Roles of MMP-9 in Pan-Cancer. Biomed Res Int 2022, 2592962, doi:10.1155/2022/2592962 (2022).
Chen, L., Zhang, J., He, Y. & Ding, X. Y. Matrix metalloproteinase-9 expression of GCTSC in peripheral tissue and central tissue of GCTB. J. Cell. Biochem. 119, 5805–5812, doi:10.1002/jcb.26766 (2018).
Huang, H. Matrix Metalloproteinase-9 (MMP-9) as a Cancer Biomarker and MMP-9 Biosensors: Recent Advances. Sensors (Basel) 18, doi:10.3390/s18103249 (2018).
Roy, R., Yang, J. & Moses, M. A. Matrix metalloproteinases as novel biomarkers and potential therapeutic targets in human cancer. J. Clin. Oncol. 27, 5287–5297, doi:10.1200/jco.2009.23.5556 (2009).
Tian, R. et al. Identification and validation of the role of matrix metalloproteinase-1 in cervical cancer. Int. J. Oncol. 52, 1198–1208, doi:10.3892/ijo.2018.4267 (2018).
Ma, B., Ran, R., Liao, H. Y. & Zhang, H. H. The paradoxical role of matrix metalloproteinase-11 in cancer. Biomed. Pharmacother. 141, 111899, doi:10.1016/j.biopha.2021.111899 (2021).
Lubowicka, E. et al. Diagnostic Power of Cytokine M-CSF, Metalloproteinase 2 (MMP-2) and Tissue Inhibitor-2 (TIMP-2) in Cervical Cancer Patients Based on ROC Analysis. Pathol. Oncol. Res. 26, 791–800, doi:10.1007/s12253-019-00626-z (2020).
Chen, D. S. & Mellman, I. Elements of cancer immunity and the cancer-immune set point. Nature 541, 321–330, doi:10.1038/nature21349 (2017).
Turley, S. J., Cremasco, V. & Astarita, J. L. Immunological hallmarks of stromal cells in the tumour microenvironment. Nat. Rev. Immunol. 15, 669–682, doi:10.1038/nri3902 (2015).
Gajewski, T. F. et al. Cancer immunotherapy strategies based on overcoming barriers within the tumor microenvironment. Curr. Opin. Immunol. 25, 268–276, doi:10.1016/j.coi.2013.02.009 (2013).
Gajewski, T. F. The Next Hurdle in Cancer Immunotherapy: Overcoming the Non-T-Cell-Inflamed Tumor Microenvironment. Semin. Oncol. 42, 663–671, doi:10.1053/j.seminoncol.2015.05.011 (2015).
Joyce, J. A. & Fearon, D. T. T cell exclusion, immune privilege, and the tumor microenvironment. Science 348, 74–80, doi:10.1126/science.aaa6204 (2015).
Salmon, H. et al. Matrix architecture defines the preferential localization and migration of T cells into the stroma of human lung tumors. J. Clin. Invest. 122, 899–910, doi:10.1172/jci45817 (2012).
De Gregorio, V., Urciuolo, F., Netti, P. A. & Imparato, G. In Vitro Organotypic Systems to Model Tumor Microenvironment in Human Papillomavirus (HPV)-Related Cancers. Cancers (Basel) 12, doi:10.3390/cancers12051150 (2020).

No competing interests reported.

Supplementarymaterial.docx

Download PDF

Version 1

posted

You are reading this latest preprint version

Integrated-analysis of single-cell RNA-seq and bulk RNA-seq reveals MMP-mediated expression patterns by distinct tumor microenvironment immune profiles in cervical cancer

Status:

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

1. Introduction

2. Materials and Methods

2.1 Single-cell sequencing

2.2 CC dataset origin and pretreatment

2.3 Unsupervised clustering of different MMP expression levels

2.4 Enrichment Analysis (GSVA) MMPs modifier mode and functional analysis of different phenotypes

2.5 Estimation of TME cell infiltration

2.6 Typing of MMP phenotype co-expression gene

2.7 Generation of MMP gene signature

2.8 Relevance of MMP gene traits to other relevant biological processes

2.9 Genomic and clinical information on immune checkpoint inhibitors

2.10 Patients and data collection

2.11 Hematoxylin-eosin staining(HE)

2.12 Quantitative real-time polymerase chain reaction (qRT-PCR)

2.13 Statistical analysis

3. Results

3.1 Single-cell transcription mapping and cell typing of the CC matrix

3.2 MMP genetic variation and expression landscape in CC

3.3 Immune cell infiltration characteristics in the TME under different MMP expression patterns

3.4 Generation of MMP gene signatures and functional annotations

3.5 Role of MMP expression patterns in immunotherapy

3.6 Characterization of MMP expression in TCGA molecular subtypes and tumor somatic cell mutations

3.7 Clinical features of MMP expression patterns

3.8 MMP expression is generally increased in CC tissue

4. Discussion

5. Conclusions

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1