The process of this study is shown in the flow chart ( Figure 1 ). The gene expression data and clinicopathological data were obtained from TCGA (https://www.cancer.gov/tcga) and GEO (https://www.ncbi.nlm.nih.gov/geo/) databases. Multiple analyses were performed using R (http://www.r-project.org, version 4.1) in this study. Sequencing data from different GEO datasets were batch-normalized using the “sva” package (20) and subsequently analyzed. The GEO accession codes of microarray data were all base on an identical platform (Affymetrix Human Genome U133 Plus 2.0 Array) and utilized in this study (GSE76039, GSE66030, GSE29265, GSE53157, and GSE65144). PDTC and ATC samples from GEO were pooled together in cross-comparison with PTC samples. Clinical information of 37 PDTC/ATC patients from GSE76039 was downloaded from cBioPortal (http://www.cbioportal.org/) (21). The median RNA sequencing value was chosen as the cut-off value of the cohorts included in this study.

The R package “limma” (22) was used to screen DEGs with the adjusted P-value <0.05. The “umap” package (23) was used to perform the sample heterogeneity analysis between PTC and PDTC/ATC. The “ComplexHeatmap” package (24) was used for customizing the heatmap and the “ggplot2” package (25) was used for visualization.

We downloaded the original single-cell RNA sequencing data from five ATC patients and six PTC patients on the GEO database (GSE148673 and GSE191288). After standard data quality control, batch effect adjustment, and normalization using the “Seurat” package (26), we clustered all cells using the Uniform Manifold Approximation and Projection (UMAP) method (27) as four basic types: tumor/epithelial cells, immune cells, endothelial cells, and fibroblasts via cell markers ( Supplementary Figure 1A ). MMP1 expression was analyzed in PTC/ATC and all kinds of cells.

The “DESeq2” and “corrplot” packages (28, 29) were used to perform a spearman correlation analysis between the expression levels of MMP1 and other different indicators, and a P-value <0.05 was selected as a cutoff criterion. The receiver operator characteristic (ROC) curve analysis was performed using the “pROC” package (30) to evaluate the diagnostic value of MMP1. Logistics regression and Cox regression with 95% confidence intervals (CIs) were performed in R to calculate the odds ratios (ORs) and hazard ratios (HRs) to assess the value of MMP1 in predicting some pathological characteristics and outcomes of PTC patients from TCGA, respectively,

The survival data were obtained from the TCGA-THCA dataset. Considering the favorable prognosis of most PTC patients, the correlation between the number of death events and overall survival (OS) was pretty low, and we focused on the progression-free survival (PFS) and disease-free interval (DFI) of the patients. The median of MMP1 expression was defined as the cutoff point for dividing the samples into high and low-expression groups. The survival probability was estimated via Statistical Product and Service Solutions (SPSS, version 26.0) using the Kaplan-Meier method, and a log-rank P-value <0.05 was considered statistically significant. Data from Gene Set Cancer Analysis (GSCA) database (31) (http://bioinfo.life.hust.edu.cn/GSCA/#/) were also adopted into survival analysis. Although there is a difference in the prognosis between ATC and PDTC that we cannot ignore, their outcomes are still very poor compared with PTC. Referring to the study (32) of Wen et al., we pooled ATC and PDTC from GSE76039 together to analyze the relationship between their prognosis and MMP1 level. Cox proportional hazard regression model was used to assess the survival difference and hazard ratio (HR) of MMP1 in PTC patients.

Enrichment analysis has been widely used in recent years to identify gene properties, based on the hypothesis that genes with similar expression profiles may be regulated by common pathways and involved in related functions (33). In this study, Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, and Gene set enrichment analysis (GSEA) were performed in different sample groupings to find potential molecular mechanisms. The “clusterProfiler” package (34) was used, and nominal P <0.05 and false discovery rate (FDR) <0.25 were selected as cutoff criteria.

Multiple gene set signature-based methods for comprehensively estimating the abundance of different immune cell types were utilized in this study. To estimate the variation of involved immune cells over PTC samples, a gene set enrichment analysis called Gene Set Variation Analysis (GSVA) was performed (35). The “ESTIMATE” package (36) was applied to calculate the immune scores of each PTC sample, and the correlation between immune cell infiltrate and the GSVA enrichment score of PTC was further verified. The immune infiltration landscape of PTC was also conducted via the “ssGSEA” algorithm in the “GSVA” package (37), “TIMER” and “xCELL” algorithms in the “IOBR” package (38, 39) to identify immune cells that might be significantly associated with MMP1.

The stemness level of a tumor could be quantified based on the RNA expression and DNA methylation signature (40), and the Epigenetically regulated RNA expression-based Stemness Score (EREG.EXPss, 103 probes) and DNA methylation-based Stemness Score (DNAss, 219 probes) were utilized in this study to assess the stemness of PTC from TCGA. BRAF ^V600E-RAS score (BRS) was developed by the TCGA group to quantify the extent to which the gene expression profile of PTC resembles either the BRAF ^V600E- or RAS-mutant profiles (41). Thyroid differentiation Score (TDS) was also developed by the TCGA group to quantify the PTC differentiation level (41) and it was calculated based on the expression levels of 16 thyroid metabolism and function genes (DIO1, DIO2, DUOX1, DUOX2, FOXE1, GLIS3, NKX2-1, PAX8, SLC26A4, SLC5A5, SLC5A8, TG, THRA, THRB, TPO, TSHR). Data of stemness scores, BRS, and TDS were obtained from TCGA and their relationships with MMP1 in PTC were also further explored.

The PTC cell lines TPC-1 and K1, PDTC cell line KTC-1, and ATC cell line CAL-62 used in this study were purchased from the American Type Culture Collection (ATCC). All cells were cultured in the 37°C and 5% CO₂ culture environment, and in specific mediums (Gibco, USA) suggested by ATCC with 10% fetal bovine serum (FBS).

The total protein of cells was extracted with protein extraction reagent Radio-immunoprecipitation Assay buffer (Beyotime, China) containing 1mM Phenylmethylsulfonyl fluoride (Beyotime, China), and quantified by the BCA Protein Assay Kit (Beyotime, China). Equal amounts of protein were subjected to 10% SDS-PAGE and then transferred to a PVDF membrane (Millipore, USA). The immunoblots were incubated with primary antibodies against MMP1 (1: 800 dilution, Proteintech, China), and GAPDH (1: 3000 dilution, Cell Signaling Technology, USA) as the internal control. The protein signals were visualized with the ChemiDoc XRS+ System (Bio-Rad, USA) using the ECL detection kit (Beyotime, China).

Short interference RNA (siRNA) for MMP1 and corresponding siRNA negative control were purchased from RiboBo (China). Transient transfection was performed using Lipofectamine 3000 (Thermo Fisher, USA) according to the manufacturer′s protocol. The transfection efficiency was evaluated using western blot and the siRNA target sequence for MMP1 was as follows: 5′- ACACAAGAGCAAGATGTGG-3′; The cells were harvested at 48 hours after transfection and then used for further experiments.

The cell viability assays were analyzed by the Cell Counting Kit-8 (CCK8) and colony formation assay to evaluate the cell proliferation of different cells according to the manufacturer’s protocol.

CCK8 assay: A total of 1000 cells were seeded into 96-well plates with different concentrations of MMP1 inhibitor Triolein (42). After 48 hours, the medium was removed and a 100 μL serum-free medium with 10% CCK8 solution (Dojindo, Japan) inside was added to each well of the plate. After incubation for 2 hours at 37°C, the spectrometric absorbance of each well at 450 nm was measured on a microplate reader (Thermo Fisher, USA).

Colony formation assay: A total of 1000 cells were seeded into 6-well plates with 50 μM Triolein. All cells were cultured in their corresponding medium, and the medium was renewed every three days over the next 15 days. Cell clones were stained with 0.2% crystal violet (Beyotime, China) and then photographed.

SPSS 26.0 (IBM, USA), R 4.1 (Lucent Technologies, USA), and Graphpad Prism 8.0 (GraphPad Software, USA) were used for statistical analysis. ImageJ 1.53k (NIH, USA) was used for colony formation counting. Unless stated otherwise, the Student’s t-test, two-way analysis of variance, or Chi-square test was performed to compare the differences between different groups, respectively. Results of P <0.05 were considered significant: NS means not significant, *P <0.05, **P <0.01, ***P <0.001, and ****P <0.0001.

Bulk gene expression data were extracted from 5 GEO datasets and samples were listed in Figure 2A . 73 advanced thyroid tumors (22 PDTCs and 51 ATCs) and 84 PTCs met the sequencing quality standards and are included in this study. The sample-to-sample heterogeneity between PDTC/ATC and PTC was detected by principal component analysis (PCA) (43)( Figure 2B ). 351 DEGs between these two groups were screened (Details in Supplementary Table 1 ) while the cutoff criteria were |log₂ (fold change)| >1. As shown in Figure 2C , MMP1 was the only highly expressed gene in PDTC/ATC relative to PTC when we chose |log₂ (fold change)| >2 as the cutoff criteria. Another 16 down-regulated genes were listed in Figure 2D , among which TG and TSHR are characteristic genes for PTC. In addition, MMP1 and other top 40 DEGs were exhibited in a heat map ( Figure 2E ).

Even if the above bulk RNA data from GEO showed MMP1 as a significantly upregulated gene in PDTC/ATC, MMP1 can be secreted by many cells including immune cells and stromal cells. To illustrate the expression of MMP1 in different cells, we used single-cell sequencing data in five ATC and six PTC patients ( Figure 2F ). After analyzing a total of 32168 cells, we found 29 clusters of different cells ( Supplementary Figure 1B ) and defined all cells as four basic types ( Figure 2G ). MMP1 was expressed mostly in tumor cells but rarely in other cells as shown in Figure 2H . Comparing its expression in PTC and ATC tumor cells, we found that MMP1 showed major expression in ATC tumor cells but lower expression in PTC tumor cells ( Figure 2I ). This reminds us that MMP1 expression in tumor cells might be related to the heterogeneity of different TC cells, and MMP1 was upregulated in ATC cells.

Data of MMP1 mRNA expression levels in PTC tissues and normal thyroid tissues were obtained from TCGA to demonstrate the role of MMP1 expression in PTC tumorigenesis. As shown in Figures 3A, B , MMP1 expression in PTC tumor tissues was significantly higher than that in normal thyroid tissues, whether they were paired samples or not. ROC curves were plotted to investigate the diagnostic value of MMP1, and the area under the curve (AUC) was 0.840 ( Figure 3C ) and 0.705 ( Figure 3D ) respectively, which can be used to distinguish the normal tissues from PTC, as well as PTC from PDTC/ATC, respectively. The clinicopathological characteristics of PTC patients from TCGA were summarized in Table 1 . MMP1 expression was significantly associated with gender (P=0.029), T stage (P =0.004), N stage (P =0.005), extrathyroidal extension (P =0.025), histological type (P <0.001) and the thyroid gland disorder history (P <0.001).

Kaplan-Meier survival analysis with a log-rank test was applied to determine the association between patients’ survival and MMP1 expression. As shown in Figures 3E, F , the MMP1 low-expressing group had significantly longer PFS and DFI in PTC patients from TCGA (P =0.043 and 0.031, respectively). Additionally, we collected OS data of 17 PDTC patients and 20 ATC patients from the GSE76039 dataset, and the prognosis of the MMP1 low-expressing group was still significantly better ( Figure 3G , P =0.001), while the level of MMP1 in ATC is significantly higher than that in PDTC (8.24 versus 3.90, P <0.01). Next, univariate logistic regression analyses were performed to explore the relationship between MMP1 and some clinicopathological characteristics of the PTC patients from TCGA ( Figure 3H ), and MMP1 was positively related with N stage (OR =1.720, 95% CI = 1.191-2.492, P =0.004), and extrathyroidal extension (OR =1.583, 95% CI =1.079-2.033, P =0.023), respectively. The HRs for DFI of PTC patients from TCGA were explored to investigate the survival significance of MMP1. Only features with P <0.1 in univariate Cox regression were included in multivariate regression analysis. M1 stage, Pathologic stage III/IV, and high MMP1 were all confirmed to be independent risk factors for DFI in PTC patients ( Table 2 ).

Multiple functional enrichment analyses were performed using the GEO cohort and TCGA cohort to explore the different pathogenesis involved in this study. GO enrichment analysis ( Figure 4A ) indicated that the DEGs between PDTC/ATC and PTC samples from the GEO cohort were mostly enriched in extracellular matrix-related pathways (Details in Supplementary Table 2 ), which were reported to be significantly associated with the MMPs family (44), including MMP1. GSEA of the DEGs was also performed ( Figure 4B ), and the top 5 pathways between PDTC/ATC and PTC were cell cycle checkpoints (normalized enrichment score (NES) =2.73, P <0.01), GTPases active formins (NES= 2.73, P <0.01), resolution of sister chromatid cohesion (NES= 2.92, P <0.01), separation of sister chromatids (NES= 2.92, P <0.01), and mitotic metaphase and anaphase (NES= 2.84, P <0.01) (45) (Details in Supplementary Table 3 ). To more precisely explore the pathway associated with MMP1, we analyzed genes associated with MMP1 expression from the TCGA PTC samples, and the top 10 positively associated genes and top 10 negatively associated genes with MMP1 were shown in the co-expression heatmap ( Figure 4C ). Genes with |spearman correlation coefficient| >2 and p <0.05 were included in the subsequent GO and KEGG analysis, and the representative pathways are shown in Figure 4D (Details in Supplementary Table 4 ).

Immune infiltration in the tumor is a complex microenvironment that interacts with tumorigenesis, progression, and metastasis (46). We further explored the infiltration of immune cells in PTC, and the correlation between the GSVA scores and immune cell infiltrates over the PTC samples from TCGA is shown in Figure 5A . In addition, we explored the relationship between infiltrating immune cells and MMP1 expression levels. ImmuneScore is an indicator to characterize the immune landscape, and it was reported to be significantly correlated with PTC progression and dedifferentiation (47). MMP1 was identified to be positively correlated with ImmuneScore in this study ( Figure 5B ), which suggested that MMP1 might be involved in immune infiltration. ssGSEA ( Figure 5C ), TIMER ( Figure 5D ), and xCELL ( Figure 5E ) were then used to explore specific immune cells associated with MMP1, and dendritic cells (DCs), macrophages, and neutrophils were immune cells confirmed to be significantly associated with MMP1 in PTC via the 3 algorithms in this study. As shown in Figures 5A, C, D , Treg cells were also found to be significantly positively correlated with MMP1. Taken together, the above results demonstrated that MMP1 might influence the infiltration of immune cells in the PTC microenvironment.

The mutation indicators BRS of PTC samples from TCGA was displayed in Figures 6A, B , and the higher MMP1, the lower BRS, which means the high MMP1 group had a higher propensity for BRAF-like mutations. In other words, MMP1 may be involved in the BRAF ^V600E- or RAS-mutant profiles in PTC. The stemness (self-renewal, differentiation, and fate determination) of a tumor could participate in multi-step tumorigenesis, recurrence, and metastasis (48). We quantified the stemness level of PTC samples from TCGA and explored its relationship with MMP1. The expression level of MMP1 was significantly positively correlated with the PTC stemness via both EREG.EXPss method ( Figure 6D ) and DNAss method ( Figure 6E ). As for the differentiation indicator, we used TDS to assess the differentiation level of PTC samples from TCGA. The results revealed that the higher MMP1 expression group had the lower TDS ( Figure 6C ). Specific associations between MMP1 and the 16 TDS-related genes were also explored, and the expression levels of MMP1 and most TDS genes were negatively correlated in both the GEO cohort of PDTC/ATC and PTC ( Figure 6F ) and TCGA cohort of PTC ( Figure 6G ). Given the above, MMP1 was considered to be involved in the dedifferentiation of PTC and could serve as a potential indicator for the differentiation level.

Since our findings above were based on RNA sequencing data obtained from public databases, we further explored and verified the expression and function of MMP1 in cell lines. PTC cell lines TPC-1 and K1, PDTC cell line KTC-1, and ATC cell line CAL-62 were used to perform western blotting assays. As shown in Figure 7A , high levels of MMP1 protein were detected in KTC-1 and CAL-62, especially in the ATC cell line CAL-62. Exogenous introduction of MMP1 siRNA was used to knock down the expressing levels of MMP1 in TPC-1, K1, KTC-1, and CAL-62 ( Figure 7B ). To further explore the potential of MMP1 as a therapeutic target in PDTC/ATC or PTC, a specific inhibitor triolein for MMP1 was used in this study. CCK8 assays were performed to determine cell viability in cell lines treated with different concentrations of triolein, and the half-maximal inhibitory concentrations (IC₅₀) of triolein in different cell lines were also detected. The IC₅₀ of triolein was 72.52 μM for TPC-1, 79.97 μM for K1, 51.80 μM for KTC-1, and 64.52 μM for CAL-62 ( Figure 7C ), respectively. Colony formation assays were also performed to analyze the cell viability of tumor cells under triolein treatment and MMP1 knockdown, where we observed that triolein and MMP1 knockdown both inhibited cell proliferation of all the TC cell lines ( Figures 7D, E ). Collectively, these results provide us with robust evidence indicating that MMP1 may serve as a novel biomarker and therapeutic target for PDTC/ATC and probably as a potential therapeutic target for PTC with the risk of progression or dedifferentiation.