Systematic identification of cancer-associated-fibroblast-derived genes in patients with colorectal cancer based on single-cell sequencing and transcriptomics

Colorectal cancer (CRC) has a high incidence rate and poor prognosis, and the available treatment approaches have limited therapeutic benefits. Therefore, understanding the underlying mechanisms of occurrence and development is particularly crucial. Increasing attention has been paid to the pathophysiological role of cancer-associated fibroblasts (CAFs) in the heterogeneous tumour microenvironment. CAFs play a crucial role in tumorigenesis, tumour progression and treatment response. However, routine tissue sequencing cannot adequately reflect the heterogeneity of tumours. In this study, single-cell sequencing was used to examine the fibroblast population in CRC. After cluster analysis, the fibroblast population was divided into four subgroups. The distribution and role of these four subgroups in CRC were found to be different. Based on differential gene expression and lasso regression analysis of the main marker genes in these subgroups, four representative genes were obtained, namely, TCF7L1, FLNA, GPX3 and MMP11. Patients with CRC were divided into the low- and high-risk groups using the prognostic risk model established based on the expression of these four genes. The prognosis of patients in different risk groups varied significantly; patients with low-risk scores had a greater response to PDL1 inhibitors, significant clinical benefits and significantly prolonged overall survival. These effects may be attributed to inhibition of the function of T cells in the immune microenvironment and promotion of the function of tumour-associated macrophages.

Introduction

As the third most common malignancy, colorectal cancer (CRC) causes more than 8% of all deaths worldwide each year (1, 2). Routine treatment of CRC includes surgery, radiotherapy and chemotherapy, which are invasive and may have a greater impact on the quality of life of patients (3). After comprehensive treatment, the 5-year survival rate of patients with early-stage CRC is 90%; however, treatment options for patients with advanced-stage CRC who are ineligible for surgery are limited (4). Immunotherapy may be beneficial for patients with advanced-stage CRC. Because of its strong anti-tumour activity, immunotherapy is used for treating several solid tumours, including melanoma, kidney cancer, non-small cell lung cancer and prostate cancer (5). In addition to targeted and anti-vascular therapies, immunotherapeutic strategies have gradually improved. PD-1/L1 and CTLA-4 are the main immunotherapeutic agents; however, their clinical efficacy remains unclear. Studies have shown that only patients with CRC with defective mismatch repair (dMMR) or high microsatellite instability (MSI-H) are eligible for checkpoint inhibition and may benefit from it (6, 7). Therefore, dMMR/MSI-H is considered a predictive biomarker for the application and efficacy of immunosuppressants (8). However, the efficacy of dMMR/MSI-H is only 30–40% (9), which considerably limits the application of immune checkpoint inhibitors for the treatment of colon cancer. Therefore, understanding the underlying mechanisms of the occurrence and development of CRC is necessary to screen for more effective predictors and improve the currently available treatment approaches.

In addition to the role of tumour cells, the tumour microenvironment (TME) is another major auxiliary factor in the onset and growth of tumours. Several studies have associated TME with the occurrence and growth of tumours, survival and clinical treatment sensitivity (10, 11). TME has an extremely complex system comprising stromal cells, tumour cells, various cytokines and an extracellular matrix (ECM) (12). Fibroblasts are the main cellular component of the matrix and are called cancer-associated fibroblasts (CAFs). They interact with cancer cells (13, 14) and are significantly associated with the prognosis of tumours (15). A recent study has demonstrated that CAFs play a significant role in various tumours. For example, matrix SOX2 upregulation promotes tumorigenesis by producing CAFs expressing SFRP1/2 (16), and Wnt-induced phenotypic transformation of CAFs inhibits EMT in CRC (17). However, most studies have focused only on the involvement of tumour cells in fibroblast remodelling or the effects of fibroblasts on tumour cells, and systematic analysis of tumours and TME including the whole fibroblast population is lacking.

In this study, we identified fibroblast subsets based on single-cell sequencing analysis and identified hub genes significantly related to fibroblasts by differential analysis, correlation analysis, univariate cox analysis and lasso cox analysis. Further, we analysed the roles of hub genes in tumors from various aspects by studying the mutations and immunity of these genes. Finally, we constructed a multi-gene signature and confirmed its role in predicting patient outcomes and immunotherapy predictions.

Materials and methods

Extraction and preprocessing of scRNA data

The read count expression profile data of 16 cancer tissues and 7 adjacent tissues were extracted from the single-cell sequencing dataset GSE200997 from the NCBI database Gene Expression Omnibus (GEO). First, the single-cell data were filtered by ensuring that each gene was expressed in at least three cells, and at least 250 genes were expressed per cell. The PercentageFeatureSet function was used to determine the proportion of mitochondria and rRNA and ensure that <3000 genes are expressed per cell and the Unique molecular identifier (UMI) of each cell is at least >100.

The data were standardised through log-normalisation, and highly variable genes were identified using the FindVariableFeatures function (variance-stabilising transformation was used to identify variable characteristics). Subsequently, the ScaleData function was used to scale all genes, and Principal components analysis (PCA) was used for dimensionality reduction to identify anchor points (dim = 40). The FindNeighbors and FindClusters functions (resolution = 0.2) were used to cluster the cells, and the RunTSNE function was used to reduce t-SNE dimensionality to screen for fibroblasts.

Extraction and preprocessing of the cancer genome atlas data

The clinical phenotype data of CRC were downloaded from TCGA database, and samples lacking data on survival time and survival status were removed. Samples were further filtered to ensure that the survival time in each sample was >0 days. In addition, the gene expression profile data were downloaded from TCGA database, and 431 tumour samples and 41 para-cancerous samples were selected for further analysis.

The copy number variation (CNV) of CRC samples were downloaded from TCGA database and integrated using the GISTIC2 software.

The single nucleotide variants (SNVs) data of TCGA-COAD cohort were downloaded from TCGA database and integrated using the Mutect2 software.

Extraction and preprocessing of GEO data

The GSE17536 and GSE17537 datasets were downloaded from GEO, and the probe IDs were converted to gene symbols according to the annotation files. A probe ID that corresponded to multiple genes was deleted, and the expression of several probes for a gene was averaged. Normal tissue samples were removed, and only tumour samples were retained. In addition, samples without clinical follow-up and OS data were removed to ensure that the survival time of all patients was >0 days. A total of 177 tumour samples and 21,655 genes were obtained from the GSE17536 dataset, and 55 tumour samples and 21,655 genes were obtained from the GSE17537 dataset.

Single-cell clustering dimensionality reduction

The R language Seurat package was first used to filter the single-cell data by setting each gene to be expressed in at least 3 cells, and each cell expresses at least 250 genes, calculating the proportion of mitochondria and rRNA through the PercentageFeatureSet function, and ensuring that each cell The expressed genes are less than 3000, and the UMI of each cell is at least greater than 100. Then, we normalized the data of 23 samples separately by log-normalization.The FindVariableFeatures function was used to find highly variable genes [identify variable features based on variance stabilizing transformation (“vst”)], then scaled all genes using the ScaleData function, and perform PCA dimensionality reduction to find anchors, we chose dim=40, pass The FindNeighbors and FindClusters functions cluster the cells (set Resolution=0.2), divided the subgroups, and used the RunTSNE function for TSNE dimensionality reduction,

Annotation and further segmentation of fibroblasts

The fibroblasts were screened with the four genes of ACTA2, FAP, PDGFRB and NOTCH3, and then the fibroblasts were extracted and clustered by the functions of FindNeighbors and FindClusters (setting Resolution=0.2), and the fibroblasts were further divided into 4 groups subpopulations and re-TSNE dimensionality reduction of fibroblasts using the RunTSNE function.

Identification of marker genes

The FindAllMarkers function of the Seurat package was used to identify marker genes of fibroblasts by LogFC=0.5, Minpct=0.35 (minimum expression ratio of differential genes) and identified marker genes with a corrected p<0.05.

Functional annotation of subgroups

KEGG enrichment analysis was performed on marker genes of fibroblast subpopulations using the compareCluster function of the clusterProfiler package in R language, and screening was performed with pvalue Cutoff=0.05.

Identification of malignant and non-malignant cells

Four fibroblast subpopulations were analyzed using the R language copykat package to differentiate between tumor cells/malignant cells and normal cells/non-malignant cells in each sample by changes in the cnv of the cells.

Copykat’s statistical workflow combines Bayesian methods with hierarchical clustering to calculate genomic copy number profiles of individual cells and to define clonal substructures from high-throughput 3’ scRNA-seq data. The workflow takes a gene expression matrix of Unique Molecular Identifier (UMI) counts as input to the calculation. Analysis begins with rows of gene annotations, ordered by their genomic coordinates. Freeman-Tukey transformation (FTT) was performed to stabilize variance, followed by polynomial dynamic linear modeling (DLM) to smooth out outliers in single-cell UMI counts. A subset of diploid cells with high confidence was then examined to infer baseline copy number values ​​for normal 2N cells. To do this, we pooled individual cells into several small hierarchical clusters and estimated the variance of each cluster using a Gaussian mixture model (GMM). By following strict classification criteria, the cluster with the smallest estimated variance was defined as “confident diploid cells”. Potential misclassification can occur when the data have only a few normal cells, or when tumor cells have near-diploid genomes and limited CNA events. In this context, Copykat provides a “GMM-defined” model to identify diploid normal cells one by one, where a mixture of three Gaussian models of gene expression in a single cell is assumed to represent genomic gain, loss, and neutral states. A single cell is defined as a “confident diploid cell” when the genes in the neutral state account for at least 99% of the expressed genes.

Tumour-related pathways.

As reported in a previous study, the 10 pathways related to tumours and genes associated with these pathways are shown in Supplementary Table 1. The scores of each cell for the 10 pathways were calculated via Single-sample GSEA (ssGSEA). The proportion of malignant and non-malignant cells and the MSI status in fibroblast subpopulations were compared via the chi-square test, and the scores of different fibroblast subpopulations associated with the 10 tumour-related pathways were compared via the Wilcoxon test.

Potential regulatory pathways of key genes

Using h.all.v7.5.1.symbols.gmt as a background, the enrichment scores of patients in TCGA cohort for each pathway were calculated using the GSVA package in R. Subsequently, the correlation between gene expression and pathway enrichment scores was analysed using the Hmisc package.

Construction of a risk model for predicting the response to PD-L1 inhibitor immunotherapy

The PD-L1 cohort (IMvigor210) was used to assess the relationship between risk scores and immunotherapy. The effects of PD-L1 inhibitors were different among 348 patients in the IMvigor210 cohort, which were characterised by stable disease (SD), progressive disease (PD), partial response (PR) and complete response (CR). In addition, differences between immunotherapy and chemotherapy were analysed in the IMvigor210 cohort. The risk model was used to evaluate the possible clinical outcomes of immunotherapy using the TIDE (http://tide.dfci.harvard.edu/) software. The likelihood of immune escape increased with increasing TIDE prediction scores, indicating that immunotherapy is less likely to benefit patients.

Statistical analyses

The Shapiro–Wilk test was used to compare the normality of variables between two groups. The unpaired Student’s t-test was used to determine the statistical significance of differences between normally distributed variables, and the Mann–Whitney U test was used to analyse non-normally distributed variables. The Kruskal–Wallis test and one-way ANOVA were employed as non-parametric and parametric methods, respectively, for comparing more than two groups. Spearman and distance correlation analyses were used to examine the correlation. The Kaplan–Meier method was used to compute survival rates, and the log-rank test was used to assess the significance of variations in survival curves.

Results

Identification of fibroblasts from scRNA-seq data

A total of 49,698 cells were obtained after filtering single-cell sequencing data. The PercentageFeatureSet function was used to calculate the proportion of mitochondria and rRNA, and 48,755 cells were obtained. As shown in Figure S1A, a significant correlation was observed between the number of UMI and mRNA but not between the number of UMI/mRNA and the content of mitochondrial genes. A violin diagram created before and after QC analysis is shown in Figures S1B, C.

Furthermore, the data of 23 samples were standardised via log-normalisation. A total of 16 subgroups were obtained, and the RunTSNE function was used to reduce t-SNE dimensionality. Fibroblasts were screened based on the expression of ACTA2, FAP, PDGFRB and NOTCH3. Because these four genes were mainly expressed by cells in subgroup 9, the cells were defined as fibroblasts (Figures S2A, B) and extracted for cluster analysis. These fibroblasts were further divided into four subgroups, and the RunTSNE function was used to reduce t-SNE dimensionality. The t-SNE map of the four fibroblast subpopulations and marker gene expression is shown in Figures S2C, D.

Figure 1A shows the t-SNE diagram of 23 samples, Figure 1B shows the t-SNE diagram of different tissues (cancer and adjacent tissues), Figure 1C shows the t-SNE diagram of the MSI status and Figure 1D shows the t-SNE diagram of fibroblast subsets after cluster analysis. The number of cells in each sample before and after data filtration is shown in Table 1.

FIGURE 1

Figure 1 (A) t-SNE diagram of 23 samples; (B) Distribution of t-SNE in cancer and adjacent tissues; (C) t-SNE distribution diagram of MSI status; (D) t-SNE map of four fibroblast subpopulations after cluster analysis; (E) Dot map of the expression of the top five marker genes in the subpopulations; (F) The proportion and cell number of subpopulations in cancer and adjacent tissues; (G) KEGG enrichment analysis of the four fibroblast subpopulations; (H) Distribution of t-SNE in malignant and non-malignant cells predicted using the copykat package.

TABLE 1

Table 1 Counting of cell counts before and after sample filtration.

The marker genes of the four subpopulations were identified using the FindAllMarkers function (logfc = 0.5 [difference multiple], Minpct = 0.35 [minimum expression ratio of different genes] and corrected p-value < 0.05). The expression of the top five marker genes with the most prominent contribution was analysed in each subgroup (Figure 1E).

Furthermore, the proportion of the four fibroblast subpopulations was analysed in each sample (Figure 1F), and the clusterprofiler package in R was used for KEGG enrichment analysis of marker genes in each subgroup (Figure 1G).

The copykat package in R was used to screen for tumour/malignant cells and normal/non-malignant cells in each sample based on CNVs (to ensure that normal cells were not included). A total of 297 cancer cells (malignant cells) and 491 normal cells (non-malignant cells) were eventually identified (Figure 1H).

Expression of fibroblasts in tumour-related pathways

Genes involved in 10 important pathways associated with tumorigenesis and development were extracted from previous studies. Figure 2A shows the enrichment of fibroblasts in the 10 tumour-related pathways. In addition, the proportion of malignant and non-malignant cells and the MSI status in the fibroblast subpopulations were compared (Figures 2B, C), and the scores of different fibroblasts in the 10 pathways were compared (Figures 2D–G).

FIGURE 2

Figure 2 (A) Heat map of the scores of 10 tumour-related pathways enriched in CAFs; (B) Comparison of CAF subpopulations in malignant and non-malignant cells; (C) Comparison of CAF subpopulations in terms of MSI status; (D) Comparison of the scores of 10 tumour-related pathways between malignant and non-malignant cells in the CAF_0 subgroup; (E) Comparison of the scores of 10 tumour-related pathways between malignant and non-malignant cells in the CAF_1 subgroup; (F) Comparison of the scores of 10 tumour-related pathways between malignant and non-malignant cells in the CAF_2 subgroup; (G) Comparison of the scores of 10 tumour-related pathways between malignant and non-malignant cells in the CAF_3 subgroup; (Wilcoxon test; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). ns, no significant.

Identification of key genes in fibroblasts

A total of 1424 upregulated and 1245 downregulated genes were identified in TCGA dataset using the limma package (FDR < 0.05 and |log2 (fold change)| > 1). Figure 3A shows a volcano map of differential analysis.

FIGURE 3

Figure 3 (A) Volcano map of differential gene expression between cancer and adjacent tissues in TCGA dataset; (B) The scores of the four fibroblast subgroups were compared between cancer and adjacent tissues (Wilcoxon test); (C) KM curve of the high- and low-score groups in the CAF_0 subgroup; (D) KM curve of the high- and low-score groups in the CAF_1 subgroup; (E) KM curve of the high- and low-score groups in the CAF_2 subgroup; (F) KM curve of the high- and low-score groups in the CAF_3 subgroup. **P < 0.01, ****P < 0.0001.

Based on the results of single-cell sequencing analysis, the scores of the CAF subgroups in TCGA dataset were calculated using ssGSEA to screen for marker genes in each subgroup. The results revealed that the scores of the CAF_0 subgroup were higher in cancer tissues, whereas those of CAF_1, CAF_2 and CAF_3 subgroups were higher in paracancerous tissues (Figure 3B). Subsequently, the survminer package was used to select optimal truncation based on the total survival time, and the scores of the four fibroblast subgroups were divided into the high- and low-score groups. The KM curve revealed that the high-score group of the four subgroups had a poor prognosis (Figures 3C–F).

Furthermore, the Hmisc package was used to examine the correlation between 2669 DEGs associated with tumorigenesis and development and the scores of the four CAF subgroups. A total of 248 key genes significantly associated with the four fibroblast subpopulations were identified (p < 0.001; cor > 0.7) and subjected to univariate cox analysis using the coxph function of the survival package. The results revealed 36 genes with a high prognostic impact, which were considered prognostic risk factors (p < 0.01) (Figure 4A).

FIGURE 4

Figure 4 (A) A total of 248 candidate DEGs were identified; (B) The locus of each independent variable changing with lambda; (C) Confidence interval under lambda; (D) Multivariate Cox regression analysis (coefficient of prognosis-related genes).

These 36 key genes were further filtered using lasso regression to decrease the number of genes used for constructing a risk model. Lasso regression is a compression estimation technique. By creating a penalty function, which causes certain coefficients to be compressed and some coefficients to be set to zero, lasso regression helps to create a more refined model. Therefore, lasso regression retains the benefit of subset contraction and is a biased estimation for analysing data with complex collinearity. It selects variables during parameter estimation and improves the method of dealing with multicollinearity in regression analysis. In this study, the R software package glmnet was used to perform lasso–Cox regression. The change in each independent variable was assessed (Figure 4B), and the number of independent variable coefficients tending to 0 was found to gradually increase with the increase in lambda. The risk model was constructed using 10-fold cross-validation, and the confidence interval of each lambda was evaluated (Figure 4C). The performance of the model was optimal at a lambda of 0.0251. The four genes obtained based on this value were selected as target genes for further analysis, and multivariate cox analysis revealed that the genes were prognostic risk factors (Figure 4D).

Mutation analysis of key genes

The SNVs of the four genes were examined in TCGA dataset, and FLNA was found to have the highest mutation frequency (Figure 5A). Subsequently, we examined the collinearity and mutual exclusiveness of these four and the top 10 genes with most mutations in CRC and found that the mutations of these four genes did not exhibit significant collinearity (Figure 5B). Furthermore, the CNVs of the four genes were analysed, and only a few samples were found to have copy number amplification/deletion (Figure 5C).

FIGURE 5

Figure 5 (A) Waterfall diagram of SNVs in the 4 key genes; (B) Collinearity and mutual exclusion analysis of the 4 key genes and 10 genes with the most mutations in CRC; (C) CNVs in the 4 key genes; (D) Heat map of the correlation between the 4 key genes and aneuploidy scores, homologous recombination defects, fraction altered, number of segments and non-silent mutation rates.

The molecular characteristics of TCGA-COAD cohort were obtained from previous pan-cancer studies. Correlation analysis revealed that MMP11 and TCF7L1 were significantly positively correlated with aneuploidy scores, homologous recombination defects and the fraction altered (Figure 5D).

Potential regulatory pathways of key genes

The enrichment scores of each pathway in TCGA cohort were calculated using the gsva package in R, and Pearson correlation analysis between the expression of the four genes and the pathway enrichment scores was performed using the Hmisc package in R. A total of 22 significantly related pathways were identified (|cor| > 0.4 and p < 0.001). Figure 6A shows a heat map of the relationship between the 4 genes and 22 pathways. Figure 6B shows a heat map of the enrichment scores of 22 pathways.

FIGURE 6

Figure 6 (A) Heat map of the correlation between genes and pathways; (B) Heat map of the enrichment scores of key pathways. *P< 0.05, **P < 0.01, ***P < 0.001.

Relationship between key genes and immunity

The immune scores of each sample in TCGA dataset were evaluated using the ESTIMATE algorithm and were found to have a significant positive correlation with the four genes (Figure 7A). The samples were divided into the high- and low-expression groups based on the median expression level of the four genes, and significant differences in immune scores were observed between the high- and low-expression groups (Figure 7B).

FIGURE 7

Figure 7 (A) Heat map of the correlation between key genes and immune scores predicted using the ESTIMATE algorithm; (B) Comparison of immune scores in the high- and low-expression groups (Wilcoxon test); (C) Heat map of the correlation between key genes and immune cell scores predicted using the CIBERSORT algorithm; (D) Comparison of the scores of 22 immune cells between the high- and low-expression groups (Wilcoxon test; *P < 0.05; **P < 0.01; ***P < 0.001; ****P < 0.0001). ns, no significant.

The CIBERSORT method was used to determine the immune cell scores of samples in TCGA dataset. Correlation analysis revealed that the expression of the four genes was significantly negatively correlated with T cell scores but was significantly positively correlated with macrophage-related scores (Figure 7C). The samples were divided into the high- and low-expression groups based on the median expression level of the four genes, and significant differences in some immune cell scores between the high- and low-expression groups (Figure 7D).

Construction of a risk model based on key genes

The results of multivariate Cox analysis are shown in Figure 4D. The risk scores of samples were calculated using the following formula: RiskScore = Σ βi × Expi, where i refers to the expression levels of the four key genes, and β is the multivariate Cox regression coefficient of the corresponding genes. The final formula for calculating risk scores based on the 4-gene signature is as follows:

$$\text{RiskScore}=0.173*\text{\hspace{0.17em}FLNA}+0.079*\text{\hspace{0.17em}MMP}11+0.146*\text{\hspace{0.17em}TCF}7\text{L}1+0.082*\text{\hspace{0.17em}GPX}3$$

TCGA cohort was used as the training dataset to determine the risk score of each sample. ROC analysis was performed to examine the efficiency of the risk model in predicting prognosis at 1–5 years using the R software package timeROC (Figure 8A). The AUC value for predicting prognosis at 4 and 5 years was 0.7. In addition, z-scores were evaluated for risk scores, and samples with risk scores of >0 were included in the high-risk group, whereas those with risk scores of <0 were included in the low-risk group. Subsequently, KM curves were plotted, and significant differences were observed between the two groups (p < 0.0001) (Figure 8B).

FIGURE 8

Figure 8 (A) ROC curve of the risk model constructed based on 4 genes in TCGA dataset; (B) KM curve of the risk model constructed based on 4 genes in TCGA dataset; (C) ROC curve of the risk model constructed based on 4 genes in the GSE17536 dataset; (D) KM curve of the risk model constructed based on 4 genes in the GSE17536 dataset; (E) ROC curve of the risk model constructed based on 4 genes in the GSE17537 dataset; (F) KM curve of the risk model constructed based on 4 genes in the GSE17537 dataset.

The GSE17536 dataset was used to verify the robustness of the model using the abovementioned method. A risk model was constructed, and its efficiency in predicting prognosis at 1–5 years was analysed using the R software package timeROC (Figure 8C). The AUC value for predicting prognosis at 1 year was 0.7. In addition, z-scores were evaluated for risk scores, and samples with risk scores of >0 were included in the high-risk group, whereas those with risk scores of <0 were included in the low-risk group. Subsequently, KM curves were plotted, and significant differences were observed between the two groups (p < 0.05) (Figure 8D).

The GSE17537 dataset was analysed using the same method. As shown in Figures 8E, F, the AUC value for predicting prognosis at 1–4 years was >0.7, and substantial differences were observed between the high- and low-risk groups.

Combination of risk scores and clinicopathological features to improve survival prediction

Multivariate and univariate Cox regression analyses of the risk score and clinicopathological features showed that the risk score was the most significant prognostic factor (Figures 9A, B). A nomogram integrating the risk scores and other clinicopathological parameters was constructed for quantifying the risk assessment and survival probability of patients with CRC (Figure 9C). The risk score had the most influence on survival rate prediction. The predictive accuracy of the risk model was further assessed using a calibration curve (Figure 9D). The calibration curve plotted for predicting prognosis at 1, 3 and 5 years and the standard curve yielded similar results, indicating the good predictive performance of the nomogram. Additionally, decision curve analysis was performed to assess the reliability of the model, and the benefits of the nomogram and risk score were found to be considerably greater than those of the extreme curve. The performance of the nomogram and risk score in predicting survival was superior to that of other clinicopathological features (Figures 9E, F).

FIGURE 9

Figure 9 (A, B) Univariate and multivariate Cox analyses of the risk score and clinicopathological features; (C) Nomogram model; (D) Calibration curve of the nomogram (1, 3 and 5 years); (E) Decision curve of the nomogram; (F) Compared with other clinicopathological features, the nomogram exhibited a superior capacity for survival prediction.

Prediction of the response to PD-L1 inhibitor immunotherapy via the risk model

The capability of the risk score to predict the response of patients to ICB therapy was assessed to study its association with immunotherapy. The results showed that patients with low risk scores had significant clinical benefits and prolonged OS in the anti-PD-L1 cohort (IMvigor210 cohort) (Figure 10C, p < 0.05). PD-L1 inhibitors had different effects among 348 patients in the IMvigor210 cohort, which were characterised by progressive disease (PD), stable disease (SD), partial response (PR) and complete response (CR). The risk scores of patients with SD/PD were higher than that of patients with other types of reactions (Figure 10A). Additionally, patients with low-risk scores experienced considerably superior treatment outcomes (Figure 10B). In addition, differences in survival among patients with different CRC stages in the IMvigor210 samples were analysed. The results revealed that stage I+II samples showed substantial survival differences (Figure 10D); however, stage III+IV samples did not show significant survival differences (Figure 10E).

FIGURE 10

Figure 10 (A) Differences in immunotherapy responses and risk scores in the IMvigor210 cohort; (B) Immunotherapy response among different risk groups in the IMvigor210 cohort; (C) Prognostic differences between different risk groups in the IMvigor210 cohort; (D) Prognostic differences between different risk groups of patients with early-stage CRC in the IMvigor210 cohort; (E) Prognostic differences between different risk groups of patients with middle- and late-stage CRC in the IMvigor210 cohort; (F) Differences in immunotherapy response and different risk scores in the IMvigor210 cohort were analysed using the TIDE software; (G) Differences in TIDE scores and immunotherapy responses in the IMvigor210 cohort; (H) Correlation analysis between the risk and TIDE scores in the IMvigor210 cohort. *P< 0.05, **P < 0.01, ****P < 0.0001.

Furthermore, differences in immunotherapy and chemotherapy responses were analysed among patients in the IMvigor210 cohort. The risk model was used to assess the potential clinical impacts of immunotherapy using the TIDE (http://tide.dfci.harvard.edu/) software. The likelihood of immune escape increased with increasing TIDE prediction scores, indicating that patients were less likely to benefit from immunotherapy. With regard to immunotherapy, the risk and TIDE scores of patients unresponsive to immunotherapy were found to be higher, which also showed that the high-risk group was less likely to benefit from immunotherapy (Figures 10F, G). In addition, Pearson correlation analysis revealed a strong positive correlation between the TIDE and risk scores (Figure 10H).

Discussion

The proliferation of connective tissue is one of the key hallmarks of tumours, and the components involved in proliferation include fibroblasts, macrophages, immune cells and dense ECM (18). Fibroblasts are the main cell type in ECM, which are called CAFs. Recently, a consensus statement was issued, which stated that cancer cells with slender morphology; a lack of mutations and negative markers of epithelial cells, endothelial cells and leukocytes may be considered CAFs (19). The characteristic markers of CAFs are α-SMA and fibroblast-activating protein (FAP), and the expression of fibroblast-specific protein 1 (FSP1), platelet-derived growth factor receptor (PDGFR)-α/β and vimentin is high in CAFs. These proteins are transcribed from ACTA2, FAP, PDGFRB and NOTCH3 genes, respectively. Because morphological features are subjective and not conducive to quantification, we used ACTA2, FAP, PDGFRB and NOTCH3 genes as markers to screen for CAFs in CRC samples via single-cell sequencing. Compared with single-cell sequencing technology, the traditional transcriptome sequencing technology (bulk RNA-seq) is based on tissue samples (cell population), which reflects the average expression level of genes in the cell population. However, several studies have indicated that CAF is heterogeneous, and certain CAF subtypes stimulate tumour growth, whereas some inhibit it. For instance, in a study by Costa et al., CAF subgroup 1 created an immunosuppressive microenvironment by suppressing CD4+CD25+ T cells in breast cancer (20). Su et al. (21) reported that the new subset, CD10+GPR77+ CAFs, can facilitate the formation of tumours in patients with breast and lung cancers. Therefore, conventional sequencing technology cannot reflect the role of CAFs in tumours. In this study, cells in subgroup 9 mainly expressed ACTA2, FAP, PDGFRB and NOTCH3 and were, therefore, defined as fibroblasts. The fibroblasts of subgroup 9 were extracted, subjected to cluster analysis and further divided into four subgroups. KEGG enrichment analysis of marker genes in each subgroup revealed that the genes were mainly enriched in pathways associated with ‘ECM’ and ‘focal adhesion’, which play an important role in tumours. However, this finding does not indicate that the four CAF subgroups play the same role in tumours.

Consistent with previous studies, this study revealed that the four CAF subpopulations may play different or contradictory roles in tumours. The distribution of malignant and non-malignant cells among the CAF subpopulations was significantly different. In the CAF_0 subpopulation, the proportion of malignant cells was higher, and that of cells with MSI-H was lower. However, in the other three subpopulations, the proportion of malignant cells was lower, and that of cells with MSI-H was higher. Furthermore, single-cell sequencing was used to screen for marker genes in the CAF subgroups, and the scores of the subgroups in TCGA dataset were calculated via ssGSEA. The results showed that the scores of the CAF_0 subgroup were higher in cancer tissues, and those of CAF_1, CAF_2 and CAF_3 subgroups were higher in adjacent tissues, which was consistent with the previous results, that is, the proportion of malignant tumour cells was higher in the in CAF_0 subgroup and lower in the other three subgroups. However, no significant differences in prognosis were observed among the four subgroups, and subgroups with high gene expression had a better prognosis. This finding indicates that CAFs in the same subgroup have some heterogeneity and hence cannot adequately predict the survival of patients in different subgroups.

Several studies have shown that CAFs promote tumour progression in various ways, such as by remodelling ECM (22, 23), interfering with drug delivery (24), producing collagen in ECM and regulating the hardness of the tumour matrix (25). CAFs can secrete chemokines (26, 27) and cytokines (28), leading to lymphatic angiogenesis (29), so as to promote the endocrine function of cancer cells. In addition, they change the immune cell environment by recruiting immunosuppressive cells and inhibiting the activity of immune effector cells (30). In this study, the role of different CAF subtypes in tumorigenesis and development of CRC was examined, and the scores of 10 tumour-related pathways in 4 CAF subtypes were compared between malignant and non-malignant cells. The PI3K pathway was found to be highly expressed in malignant cells. Studies have shown that the PI3K pathway promotes tumour progression. The EphA2–PI3K signal can simulate angiogenesis induced by CAFs in gastric cancer cells (31). CAF-derived HGF promotes cell proliferation and drug resistance by upregulating the c-Met/PI3K/Akt and GRP78 signalling pathways in ovarian cancer cells (32). The results of this study are consistent with those of previous studies, suggesting that CAFs promote tumour progression through the PI3K pathway.

To decrease the heterogeneity among subgroups, the marker genes of different CAF subgroups were used to classify CAFs. After differential expression analysis, four genes were selected via lasso regression analysis, namely, TCF7L1, FLNA, GPX3 and MMP11. TCF7L1 is a member of the TCF/lymphoid enhancer (LEF) family of transcription factors, which is involved in maintaining stem cell pluripotency (33) and skin epithelial tissue homeostasis (34). Studies have shown that ectopic TCF7L1 expression impairs the growth and invasion of highly metastatic breast cancer cells (35). In addition, overexpression of TCF7L1 can induce the growth of colorectal tumour cells (36). FLNA, the most abundant and widely distributed member of the filamin family, is a non-muscle actin filament cross-linked protein (37). Some studies have shown that FLNA is associated with multiple functional non-cytoskeletal proteins and participates in several related pathways regulating cell migration and adhesion (38). FLNA acts as a pro-oncoprotein in various human malignancies, including metastatic melanoma and hepatocellular carcinoma (39, 40). However, the expression of FLNA is decreased in breast cancer, which is negatively correlated with lymph node metastasis. FLNA knockout can promote cell migration and invasion (41). In CRC, FLNA promotes chemotherapy resistance by inducing epithelial–mesenchymal transformation and the Smad2 signalling pathway (42). Therefore, the controversial role of FLNA in human malignant tumours has been reported in several studies. GPX3 is a tumour suppressor gene and the main antioxidant enzyme in plasma. It plays an important role in scavenging hydrogen peroxide and other oxygen free radicals and protecting cells from oxidative stress-induced damage (43–45). As an important member of the MMP family, MMP11 regulates a series of physiological processes and signalling events, manipulates some bioactive molecules on the cell surface, changes the biological behaviour of cells and plays an important role in TME (46, 47). In addition, studies have shown that MMP is closely related to tumorigenesis. The most important MMP is MMP11, which is overexpressed in tumours and participates in the proliferation and malignant development of tumour cells (48, 49). However, according to previous studies, CAFs can also degrade ECM by releasing MMPs and synthesising new matrix proteins to provide structural support for tumour invasion and angiogenesis (50, 51). Therefore, MMP11 can be used for the evaluation of prognosis.

The four genes identified via lasso regression were subjected to enrichment analysis, and 22 significantly related pathways were identified including those associated with ‘angiogenesis’, ‘apical junction’, ‘apoptosis’ and ‘IL2–STAT5’. The four key genes were used to establish a prognostic risk model, which had good stability and accuracy in predicting prognosis in both training and validation sets. The prognosis of patients in the high-risk group was worse. To quantify the risk assessment and survival probability of patients, the risk score was combined with other clinicopathological features, and it was found that the risk score adequately predicted clinicopathological features, especially the M stage, indicating that patients with high risk scores may be more predisposed to distant metastasis. In addition, to examine the relationship between the risk score and immunotherapy, the ability of risk score to predict the response of patients to ICB therapy was examined. Patients with low risk scores had significant clinical benefits and significantly prolonged OS in the anti-PD-L1 cohort. Furthermore, mutation analysis of the four genes in TCGA cohort revealed that FLNA had the highest mutation frequency, and there was no significant collinearity among the mutations of the four genes. Moreover, only a few samples had copy number amplification/deletion. Because the mutation frequency of the four genes is not significant, their role may be directly realised through their expression levels.

Furthermore, the correlation between the prognostic risk model and infiltrating immune cells was analysed, and a significant positive correlation was observed between the four genes and immune scores, indicating that high gene expression increased the abundance of infiltrating immune cells in ECM. Moreover, these four genes had a significant negative correlation with T cell-related scores. Therefore, CAFs labelled by these genes can promote tumour progression by inhibiting T-cell function. This result is consistent with that of previous studies. CAFs can induce immune evasion of cancer cells (52, 53) and restrict the recruitment of immune effector cells (such as CD8+ T cells) to tumour tissues by secreting different chemokines (54). In this study, a significant positive correlation was observed between the four genes and the score of macrophages, which is consistent with the finding of a previously reported study, indicating that CAF can induce M2 polarisation (55). These results suggest that the interaction between stromal cells and immune-related cells in TME promotes tumour progression.

However, this study has certain limitations. First, the results of single-cell sequencing were not verified in actual clinical samples. The screened key genes lack basic in vivo and in vitro experimental verification, and the prognostic model should be verified in actual clinical samples, which is our next research direction. In addition, there are some contradictory and unexplained results. For example, the distribution of different CAF subpopulations among malignant and normal cells is different; however, the prognosis among these populations was not different. Whether their distribution in malignant cells also plays an important role warrants further investigation and verification.

In conclusion, the fibroblast population screened via single-cell sequencing in CRC was divided into four subpopulations through cluster analysis. The distribution and role of these four subpopulations are different in CRC. In addition, by analysing the differential expression of the main marker genes in these subpopulations, four representative genes were identified via lasso regression, namely, TCF7L1, FLNA, GPX3 and MMP11. Using the prognostic risk model constructed based on the expression of these four genes, patients with CRC were divided into the high- and low-risk groups. Patients with low risk scores had significant clinical benefits from immunotherapy and had significantly prolonged OS, which may be attributed to inhibition of T-cell function in the immune microenvironment and promotion of the function of tumour-associated macrophages.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

JZ performed the study and wrote the paper. YC edited and proofread the paper. All authors contributed to the article and approved the submitted version.

Conflict of interest

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.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2022.988246/full#supplementary-material

Supplementary Figure 1 | (A) Correlation between mitochondrial genes and UMI/mRNA quantity as well as between UMI and mRNA quantity; (B) Correlation among the mRNA/UMI/mitochondrial content/rRNA content of samples before filtration; (C) Correlation among the mRNA/UMI/mitochondrial content/rRNA content of samples after filtration; (D) Dimensionality reduction and identification of anchor points via PCA.

Supplementary Figure 2 | (A) Distribution of subpopulations of all cells after cluster analysis; (B) t-SNE map of marker gene expression in fibroblasts; (C) Distribution of fibroblast subgroups after re-clustering; (D) t-SNE map of marker gene expression in four small fibroblast subpopulations.

References

1. Castells A. Hereditary forms of colorectal cancer. Gastroenterol y Hepatol (2016) 39 (Suppl 1):62–7. doi: 10.1016/S0210-5705(16)30176-5

CrossRef Full Text | Google Scholar

2. Ma H, Brosens LAA, Offerhaus GJA, Giardiello FM, de Leng WWJ, Montgomery EA. Pathology and genetics of hereditary colorectal cancer. Pathology (2018) 50(1):49–59. doi: 10.1016/j.pathol.2017.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Venook A. Critical evaluation of current treatments in metastatic colorectal cancer. Oncol (2005) 10(4):250–61. doi: 10.1634/theoncologist.10-4-250

CrossRef Full Text | Google Scholar

4. Chibaudel B, Tournigand C, Bonnetain F, Richa H, Benetkiewicz M, André T, et al. Therapeutic strategy in unresectable metastatic colorectal cancer: An updated review. Ther Adv Med Oncol (2015) 7(3):153–69. doi: 10.1177/1758834015572343

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Topalian SL, Hodi FS, Brahmer JR, Gettinger SN, Smith DC, McDermott DF, et al. Safety, activity, and immune correlates of anti-PD-1 antibody in cancer. N Engl J Med (2012) 366(26):2443–54. doi: 10.1056/NEJMoa1200690

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Lichtenstern CR, Ngu RK, Shalapour S, Karin M. Immunotherapy, inflammation and colorectal cancer. Cells (2020) 9(3):618. doi: 10.3390/cells9030618

CrossRef Full Text | Google Scholar

7. Asaoka Y, Ijichi H, Koike K. PD-1 blockade in tumors with mismatch-repair deficiency. N Engl J Med (2015) 373(20):1979. doi: 10.1056/NEJMoa1500596

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Ganesh K, Stadler ZK, Cercek A, Mendelsohn RB, Shia J, Segal NH, et al. Immunotherapy in colorectal cancer: rationale, challenges and potential. Nat Rev Gastroenterol Hepatol (2019) 16(6):361–75. doi: 10.1038/s41575-019-0126-x

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Le DT, Durham JN, Smith KN, Wang H, Bartlett BR, Aulakh LK, et al. Mismatch repair deficiency predicts response of solid tumors to PD-1 blockade. Sci (New York NY) (2017) 357(6349):409–13. doi: 10.1126/science.aan6733

CrossRef Full Text | Google Scholar

10. Demircan NC, Boussios S, Tasci T, Öztürk MA. Current and future immunotherapy approaches in ovarian cancer. Ann Trans Med (2020) 8(24):1714. doi: 10.21037/atm-20-4499

CrossRef Full Text | Google Scholar

11. Angell HK, Bruni D, Barrett JC, Herbst R, Galon J. The immunoscore: Colon cancer and beyond. Clin Cancer Res: An Off J Am Assoc Cancer Res (2020) 26(2):332–9. doi: 10.1158/1078-0432.CCR-18-1851

CrossRef Full Text | Google Scholar

12. Catalano V, Turdo A, Di Franco S, Dieli F, Todaro M, Stassi G. Tumor and its microenvironment: A synergistic interplay. Semin Cancer Biol (2013) 23(6 Pt B):522–32. doi: 10.1016/j.semcancer.2013.08.007

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Xing F, Saidou J, Watabe K. Cancer associated fibroblasts (CAFs) in tumor microenvironment. Front biosci (Landmark edition) (2010) 15(1):166–79. doi: 10.2741/3613

CrossRef Full Text | Google Scholar

14. Cirri P, Chiarugi P. Cancer associated fibroblasts: The dark side of the coin. Am J Cancer Res (2011) 1(4):482–97.

PubMed Abstract | Google Scholar

15. Paulsson J, Micke P. Prognostic relevance of cancer-associated fibroblasts in human cancer. Semin Cancer Biol (2014) 25:61–8. doi: 10.1016/j.semcancer.2014.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Kasashima H, Duran A, Martinez-Ordoñez A, Nakanishi Y, Kinoshita H, Linares JF, et al. Stromal SOX2 upregulation promotes tumorigenesis through the generation of a SFRP1/2-expressing cancer-associated fibroblast population. Dev Cell (2021) 56(1):95–110.e110. doi: 10.1016/j.devcel.2020.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Mosa MH, Michels BE, Menche C, Nicolas AM, Darvishi T, Greten FR, et al. A wnt-induced phenotypic switch in cancer-associated fibroblasts inhibits EMT in colorectal cancer. Cancer Res (2020) 80(24):5569–82. doi: 10.1158/0008-5472.CAN-20-0263

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Feig C, Gopinathan A, Neesse A, Chan DS, Cook N, Tuveson DA. The pancreas cancer microenvironment. Clin Cancer Res: An Off J Am Assoc Cancer Res (2012) 18(16):4266–76. doi: 10.1158/1078-0432.CCR-11-3114

CrossRef Full Text | Google Scholar

19. Giaquinto AN, Miller KD, Tossas KY, Winn RA, Jemal A, Siegel RL. Cancer statistics for African American/Black people 2022. CA: Cancer J Clin (2022) 72(3):202–29. doi: 10.3322/caac.21718

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Costa A, Kieffer Y, Scholer-Dahirel A, Pelon F, Bourachot B, Cardon M, et al. Fibroblast heterogeneity and immunosuppressive environment in human breast cancer. Cancer Cell (2018) 33(3):463–479.e410. doi: 10.1016/j.ccell.2018.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Su S, Chen J, Yao H, Liu J, Yu S, Lao L, et al. CD10(+)GPR77(+) cancer-associated fibroblasts promote cancer formation and chemoresistance by sustaining cancer stemness. Cell (2018) 172(4):841–856.e816. doi: 10.1016/j.cell.2018.01.009

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Gaggioli C, Hooper S, Hidalgo-Carcedo C, Grosse R, Marshall JF, Harrington K, et al. Fibroblast-led collective invasion of carcinoma cells with differing roles for RhoGTPases in leading and following cells. Nat Cell Biol (2007) 9(12):1392–400. doi: 10.1038/ncb1658

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Goetz JG, Minguet S, Navarro-Lérida I, Lazcano JJ, Samaniego R, Calvo E, et al. Biomechanical remodeling of the microenvironment by stromal caveolin-1 favors tumor invasion and metastasis. Cell (2011) 146(1):148–63. doi: 10.1016/j.cell.2011.05.040

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Jacobetz MA, Chan DS, Neesse A, Bapiro TE, Cook N, Frese KK, et al. Hyaluronan impairs vascular function and drug delivery in a mouse model of pancreatic cancer. Gut (2013) 62(1):112–20. doi: 10.1136/gutjnl-2012-302529

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Olivares O, Mayers JR, Gouirand V, Torrence ME, Gicquel T, Borge L, et al. Collagen-derived proline promotes pancreatic ductal adenocarcinoma cell survival under nutrient limited conditions. Nat Commun (2017) 8:16031. doi: 10.1038/ncomms16031

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Scherz-Shouval R, Santagata S, Mendillo ML, Sholl LM, Ben-Aharon I, Beck AH, et al. The reprogramming of tumor stroma by HSF1 is a potent enabler of malignancy. Cell (2014) 158(3):564–78. doi: 10.1016/j.cell.2014.05.045

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Kumar V, Donthireddy L, Marvel D, Condamine T, Wang F, Lavilla-Alonso S, et al. Cancer-associated fibroblasts neutralize the anti-tumor effect of CSF1 receptor blockade by inducing PMN-MDSC infiltration of tumors. Cancer Cell (2017) 32(5):654–668.e655. doi: 10.1016/j.ccell.2017.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Biffi G, Oni TE, Spielman B, Hao Y, Elyada E, Park Y, et al. IL1-induced JAK/STAT signaling is antagonized by TGFβ to shape CAF heterogeneity in pancreatic ductal adenocarcinoma. Cancer Discov (2019) 9(2):282–301. doi: 10.1158/2159-8290.CD-18-0710

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Cadamuro M, Brivio S, Mertens J, Vismara M, Moncsek A, Milani C, et al. Platelet-derived growth factor-d enables liver myofibroblasts to promote tumor lymphangiogenesis in cholangiocarcinoma. J Hepatol (2019) 70(4):700–9. doi: 10.1016/j.jhep.2018.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Öhlund D, Handly-Santana A, Biffi G, Elyada E, Almeida AS, Ponz-Sarvise M, et al. Distinct populations of inflammatory fibroblasts and myofibroblasts in pancreatic cancer. J Exp Med (2017) 214(3):579–96. doi: 10.1084/jem.20162024

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Kim HS, Won YJ, Shim JH, Kim HJ, Kim BS, Hong HN. Role of EphA2-PI3K signaling in vasculogenic mimicry induced by cancer-associated fibroblasts in gastric cancer cells. Oncol Lett (2019) 18(3):3031–8. doi: 10.3892/ol.2019.10677

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Deying W, Feng G, Shumei L, Hui Z, Ming L, Hongqing W. CAF-derived HGF promotes cell proliferation and drug resistance by up-regulating the c-Met/PI3K/Akt and GRP78 signalling in ovarian cancer cells. Biosci Rep (2017) 37(2):BSR20160470. doi: 10.1042/BSR20160470

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Yi F, Pereira L, Hoffman JA, Shy BR, Yuen CM, Liu DR, et al. Opposing effects of Tcf3 and Tcf1 control wnt stimulation of embryonic stem cell self-renewal. Nat Cell Biol (2011) 13(7):762–70. doi: 10.1038/ncb2283

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Nguyen H, Merrill BJ, Polak L, Nikolova M, Rendl M, Shaver TM, et al. Tcf3 and Tcf4 are essential for long-term homeostasis of skin epithelia. Nat Genet (2009) 41(10):1068–75. doi: 10.1038/ng.431

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Slyper M, Shahar A, Bar-Ziv A, Granit RZ, Hamburger T, Maly B, et al. Control of breast cancer growth and initiation by the stem cell-associated transcription factor TCF3. Cancer Res (2012) 72(21):5613–24. doi: 10.1158/0008-5472.CAN-12-0119

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Murphy M, Chatterjee SS, Jain S, Katari M, DasGupta R. TCF7L1 modulates colorectal cancer growth by inhibiting expression of the tumor-suppressor gene EPHB3. Sci Rep (2016) 6:28299. doi: 10.1038/srep28299

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Hartwig JH, Stossel TP. Isolation and properties of actin, myosin, and a new actinbinding protein in rabbit alveolar macrophages. J Biol Chem (1975) 250(14):5696–705. doi: 10.1016/S0021-9258(19)41235-0

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Zhou X, Borén J, Akyürek LM. Filamins in cardiovascular development. Trends Cardiovasc Med (2007) 17(7):222–9. doi: 10.1016/j.tcm.2007.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Zhang K, Zhu T, Gao D, Zhang Y, Zhao Q, Liu S, et al. Filamin a expression correlates with proliferation and invasive properties of human metastatic melanoma tumors: Implications for survival in patients. J Cancer Res Clin Oncol (2014) 140(11):1913–26. doi: 10.1007/s00432-014-1722-3

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Kircher P, Hermanns C, Nossek M, Drexler MK, Grosse R, Fischer M, et al. Filamin a interacts with the coactivator MKL1 to promote the activity of the transcription factor SRF and cell migration. Sci Signaling (2015) 8(402):ra112. doi: 10.1126/scisignal.aad2959

CrossRef Full Text | Google Scholar

41. Xu Y, Bismar TA, Su J, Xu B, Kristiansen G, Varga Z, et al. Filamin a regulates focal adhesion disassembly and suppresses breast cancer cell migration and invasion. J Exp Med (2010) 207(11):2421–37. doi: 10.1084/jem.20100433

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Cheng M, Jiang Y, Yang H, Zhao D, Li L, Liu X. FLNA promotes chemoresistance of colorectal cancer through inducing epithelial-mesenchymal transition and smad2 signaling pathway. Am J Cancer Res (2020) 10(2):403–23.

PubMed Abstract | Google Scholar

43. Baez-Duarte BG, Mendoza-Carrera F, García-Zapién A, Flores-Martínez SE, Sánchez-Corona J, Zamora-Ginez I, et al. Glutathione peroxidase 3 serum levels and GPX3 gene polymorphisms in subjects with metabolic syndrome. Arch Med Res (2014) 45(5):375–82. doi: 10.1016/j.arcmed.2014.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Liu K, Jin M, Xiao L, Liu H, Wei S. Distinct prognostic values of mRNA expression of glutathione peroxidases in non-small cell lung cancer. Cancer Manage Res (2018) 10:2997–3005. doi: 10.2147/CMAR.S163432

CrossRef Full Text | Google Scholar

45. Min SY, Kim HS, Jung EJ, Jung EJ, Jee CD, Kim WH. Prognostic significance of glutathione peroxidase 1 (GPX1) down-regulation and correlation with aberrant promoter methylation in human gastric cancer. Anticancer Res (2012) 32(8):3169–75.

PubMed Abstract | Google Scholar

46. Pittayapruek P, Meephansan J, Prapapan O, Komine M, Ohtsuki M. Role of matrix metalloproteinases in photoaging and photocarcinogenesis. Int J Mol Sci (2016) 17(6):868. doi: 10.3390/ijms17060868

CrossRef Full Text | Google Scholar

47. Motrescu ER, Rio MC. Cancer cells, adipocytes and matrix metalloproteinase 11: A vicious tumor progression cycle. Biol Chem (2008) 389(8):1037–41. doi: 10.1515/BC.2008.110

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Scheau C, Badarau IA, Costache R, Caruntu C, Mihai GL, Didilescu AC, et al. The role of matrix metalloproteinases in the epithelial-mesenchymal transition of hepatocellular carcinoma. Analytical Cell Pathol (Amsterdam) (2019) 2019:9423907. doi: 10.1155/2019/9423907

CrossRef Full Text | Google Scholar

49. Zhang X, Huang S, Guo J, Zhou L, You L, Zhang T, et al. Insights into the distinct roles of MMP-11 in tumor biology and future therapeutics (Review). Int J Oncol (2016) 48(5):1783–93. doi: 10.3892/ijo.2016.3400

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Fullár A, Dudás J, Oláh L, Hollósi P, Papp Z, Sobel G, et al. Remodeling of extracellular matrix by normal and tumor-associated fibroblasts promotes cervical cancer progression. BMC Cancer (2015) 15:256. doi: 10.1186/s12885-015-1272-3

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Eble JA, Niland S. The extracellular matrix in tumor progression and metastasis. Clin Exp Metastasis (2019) 36(3):171–98. doi: 10.1007/s10585-019-09966-1

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Martinez-Outschoorn UE, Lisanti MP, Sotgia F. Catabolic cancer-associated fibroblasts transfer energy and biomass to anabolic cancer cells, fueling tumor growth. Semin Cancer Biol (2014) 25:47–60. doi: 10.1016/j.semcancer.2014.01.005

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Farhood B, Najafi M, Mortezaee K. Cancer-associated fibroblasts: Secretions, interactions, and therapy. J Cell Biochem (2019) 120(3):2791–800. doi: 10.1002/jcb.27703

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Ene-Obong A, Clear AJ, Watt J, Wang J, Fatah R, Riches JC, et al. Activated pancreatic stellate cells sequester CD8+ T cells to reduce their infiltration of the juxtatumoral compartment of pancreatic ductal adenocarcinoma. Gastroenterology (2013) 145(5):1121–32. doi: 10.1053/j.gastro.2013.07.025

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Zhang A, Qian Y, Ye Z, Chen H, Xie H, Zhou L, et al. Cancer-associated fibroblasts promote M2 polarization of macrophages in pancreatic ductal adenocarcinoma. Cancer Med (2017) 6(2):463–70. doi: 10.1002/cam4.993

PubMed Abstract | CrossRef Full Text | Google Scholar