Integrative Analysis From Multicenter Studies Identifies a WGCNA-Derived Cancer-Associated Fibroblast Signature for Ovarian Cancer.

Cancer-associated fibroblasts (CAFs) are a major contributor to tumor stromal crosstalk in the tumor microenvironment (TME) and boost tumor progression by promoting angiogenesis and lymphangiogenesis. This study aimed to identify prognostic genes associated with CAFs that lead to high morbidity and mortality in ovarian cancer (OC) patients. We performed bioinformatics analysis in 16 multicenter studies (2,742 patients) and identified CAF-associated hub genes using the weighted gene co-expression network analysis (WGCNA). A machine learning methodology was used to identify COL16A1, COL5A2, GREM1, LUM, SRPX, and TIMP3 and construct a prognostic signature. Subsequently, a series of bioinformatics algorithms indicated risk stratification based on the above signature, suggesting that high-risk patients have a worse prognosis, weaker immune response, and lower tumor mutational burden (TMB) status but may be more sensitive to routine chemotherapeutic agents. Finally, we characterized prognostic markers using cell lines, immunohistochemistry, and single-cell sequencing. In conclusion, these results suggest that the CAF-related signature may be a novel pretreatment guide for anti-CAFs, and prognostic markers in CAFs may be potential therapeutic targets to inhibit OC progression.

Introduction

Cancer-associated fibroblasts (CAFs) play a key role in the tumor microenvironment (TME) and influence tumor progression and metastasis through multiple pathways, including remodeling of the extracellular matrix (ECM), producing growth factors, and promoting angiogenesis (1). Meanwhile, ovarian cancer (OC) is a heterogeneous disease characterized by a propensity for peritoneal spread. Due to the complex interconnected signaling network and the unique peritoneal TME, cancer cells can interact with CAFs, adipocytes, immune cells, and chemokines (2). As a result, tumor migration and immune evasion frequently occur in OC patients, and immunotherapy has little effect (3).

The ECM is composed and reconstituted by CAFs, a barrier that supports tumor cell invasion and inhibits infiltration of antitumor immune cells, thus leading to immune evasion and chemoresistance (4, 5). Several researchers have explored different CAF subgroups with varying CAF marker expressions, such as alpha-smooth muscle actin (α-SMA), fibronectin attachment protein (FAP), and pelleted growth factor receptor (PDGFR) (6, 7). For example, in oral cancer, WNT2+ CAFs were negatively correlated with CD8+ T-cell activity (8). In pancreatic cancer, knocking down α-SMA+ CAFs unexpectedly enhanced tumor infiltration and increased Regulatory T cells (Tregs) abundance, leading to enhanced disease progression and reduced survival rates in mice (9). In breast and colon cancer, DNA-based vaccines targeting FAP induced the killing of CAFs by CD8+ T cells (10). Therefore, targeting the CAF-mediated immunosuppressive stromal microenvironment in combination with immunotherapy is expected to improve immune checkpoint inhibitor (ICI) response (11).

Weighted gene co-expression network analysis (WGCNA) is a systematic bioinformatics algorithm that enables the integration of highly correlated genes into several modules (12). This is a novel method to explore the relationship between numerous genes and clinical phenotypes. WGCNA has been applied to identify CAF markers, such as in gastric cancer (13), bladder cancer (14), and renal cell carcinoma (15). However, to date, CAFs have not been analyzed by WGCNA in large-sample multicenter OC cohorts. In this study, we integrated 16 multicenter studies that included 2,742 patients with complete follow-up information for bioinformatics analysis. We explored the hub modules most relevant for CAF infiltration and identified COL16A1, COL5A2, GREM1, LUM, SRPX, and TIMP3 as prognostic CAF markers. Subsequently, CAF signatures capable of predicting prognosis and treatment response were constructed, and the predictive ability was validated in multiple cohorts. In addition, we characterized markers using cell lines, immunohistochemistry, and single-cell sequencing. In conclusion, our results imply that the CAF signature may be a novel anti-CAF therapeutic approach in OC.

Methods

Datasets and Data Preprocessing

The fragments per kilobase of transcript per million mapped reads (FPKM) format RNA sequencing (RNA-seq) data with complete follow-up information of 372 samples were downloaded from The Cancer Genome Atlas (TCGA) database (16). Except for the samples without survival follow-up information, we still retained the samples with other clinical information missing. The somatic mutation data were also acquired from TCGA database. The tumor mutational burden (TMB) value of each sample was calculated via the tmb algorithm in the “maftools” package (17). We performed log2 [transcripts per million (TPM) + 1] transformation on the above raw data (18). In the Gene Expression Omnibus (GEO) database (19), we integrated the multiple datasets (RNA-seq or microarray) based on the GPL570 platform (GSE19829, GSE18520, GSE9891, GSE26193, GSE30161, and GSE63885; n = 597), GPL96 platform (GSE3149, GSE23554, GSE26712, and GSE14764; n = 409), GPL7759 platform (GSE13876, n = 415), GPL2986 platform (GSE49997, n = 194), and GPL14951 platform (GSE140082, n = 380). In addition, we downloaded anti-programmed death-1 (PD-1) dataset (IMvigor, n = 348) and anti-PD-L1 dataset (GSE78220, n = 27) based on immunotherapy. Cell line RNA-seq data from 47 fibroblast origins and 37 OC origins were acquired from the Cancer Cell Line Encyclopedia (CCLE) database (20). Immunohistochemical (IHC) staining images in OC tissues were downloaded from the Human Protein Atlas (HPA) database (21). Batch effects from meta-cohorts (GPL570 or GPL96) were corrected using the ComBat algorithm in the “sva” package (22). CAF markers were collected from previous references (23).

In conclusion, we integrated 16 multicenter studies and included 2,742 patients with complete follow-up information for our bioinformatics analysis.

Cancer-Associated Fibroblasts and Stromal Analysis

CAF abundances and stromal scores were computed using four methods: Estimate the Proportion of Immune and Cancer cells (EPIC) algorithm (24), xCell algorithm (25), microenvironment cell populations-counter (MCP-counter) algorithm (26), and Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression data (ESTIMATE) algorithm (27). We used “IOBR” package to invoke the above algorithm (28). The CAF abundances calculated by EPIC and MCP-counter were defined as phenotypic data for subsequent WGCNA. The data calculated by other algorithms were used for validation.

Weighted Gene Co-Expression Network Analysis

The “WGCNA” package (12) screened hub genes that were significantly associated with CAF scores. The expression profiles of the top 25% of the variance in the GPL570 meta-cohort and TCGA-OV cohort first were as the input. Then, according to our previous study (29), a soft threshold was determined, an adjacency matrix was clustered, and a hub module was determined. The strongest positive correlation was selected for further analysis by calculating the Pearson correlation coefficient between the modules and CAF scores. Then, we measured gene significance (GS) for each gene’s traits and module membership (MM) in the hub module. Finally, genes in the module were screened as potential CAF-related genes using MM >0.6 and GS >0.6 as thresholds.

Enrichment Analysis

The h.all.v7.4.symbols gene set was downloaded from the MSigDB database (30) for enrichment analysis in “GSVA” package (31). The adj.P value <0.05 was considered statistically significant. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were conducted using “clusterProfiler” package (32). The adj.P-value <0.05 and adj.q-value <0.05 were considered statistically significant.

Construction and Validation of the Cancer-Associated Fibroblast Signature

The GPL570 meta-cohort with a larger sample size was used as the training cohort, and other cohorts were used as the validation cohort. Univariate Cox regression analysis was performed on common hub genes in 16 multicenter studies (P-value <0.05). In the least absolute shrinkage and selection operator (LASSO) regression analysis (33), 1,000 iterations were performed to reduce the genes, and subsequently, the above genes were subjected to multivariate Cox regression analysis to obtain the coefficients. CAF risk score was derived using the same formula as in our previous study (34, 35). The OC patients in each cohort were divided into high-risk and low-risk groups, and the cutoff value for each cohort was used as the threshold.

Chemotherapy Response Predictions

The “pRRophetic” package (36) was used to calculate IC50 value (bleomycin, cisplatin, docetaxel, gemcitabine, doxorubicin, and etoposide) of different samples based on gene expression.

Single-Cell Sequencing Analysis

We analyzed single-cell RNA-sequencing (scRNA-seq) data (GSE118828) from OC tissues based on the Tumor Immune Single Cell Hub (TISCH) database (37), and the whole cells were annotated into six clusters: fibroblasts, myofibroblasts, endothelial, malignant, monocyte or macrophage (Mono/Macro), and conventional CD4 T cell (CD4Tconv).

Immunofluorescence Staining

In total, two formalin-fixed paraffin-embedded (FFPE) tissue (primary tumors and recurrent tumors) blocks were selected from the Zhongda Hospital Southeast University. Immunofluorescence staining was also done by a commercial entity (Servicebio, Wuhan, Hubei, China). According to the company, detailed methods are available in a previous publication (38). Antibody staining order always remains the same, all slices with 4,6-diamidino-2-phenylindole (DAPI) (Servicebio, G1012) finally after dyeing. Monoclonal antibodies in immunofluorescence panels were CD8 (Servicebio, GB13068-2, 1:500), Foxp3 (Servicebio, GB112325, 1:3,000), and α-SMA (Servicebio, GB13044, 1:1,000). Slices were placed under a scanner to collect images, and image data were obtained using CaseViewer software.

Statistical Analysis

All statistical analyses were performed using the R software (v.4.0.1). The Wilcoxon test was applied for pairwise comparisons. The Kaplan–Meier analysis with the log-rank test was adopted for overall survival comparisons. More detailed statistical methods for transcriptome data processing are covered in the above section. P < 0.05 was considered statistically significant.

Results

Cancer-Associated Fibroblasts and Stromal Score Could Be Considered Prognostic Markers for Ovarian Cancer

We integrated the multidatasets based on the GPL570 platform (GSE19829, GSE18520, GSE9891, GSE26193, GSE30161, and GSE63885), and the Uniform Manifold Approximation and Projection (UMAP) analysis showed the distribution of each dataset before and after removal of batch effect ( ). The expression density plot also revealed that the batch effect of the GPL570 meta-cohort was well removed ( ). Finally, we normalized the expression profiles of 597 samples with complete follow-up information ( ). Previous references have reported the ability of CAFs to recruit immunosuppressive cells, so we performed immunofluorescence staining using immunofluorescence in patients with primary tumors and in patients with recurrent tumors (39). Interestingly, there was a recruitment of Treg cells (green) around the CAF cells (pink) in patients with primary tumors ( ), especially in the recurrent samples, where a larger number of Treg cells clustered to the prominent part of the CAF cells ( ). CD8+ cells (red) were rarely seen around CAF cells in both samples. Subsequently, the CAF infiltration score was predicted by EPIC, xCell, and MCP-counter algorithms based on the GPL570 meta-cohort (n = 597) and TCGA-OV cohort (n = 372), and the stromal score was calculated by ESTIMATE algorithm. We divided all samples into a high CAF/stromal score group and a low CAF/stromal group according to the cutoff values of the scores calculated by the four bioinformatics algorithms. In the GPL570 meta-cohort, the results showed that higher CAF infiltration and stromal score were significantly associated with poorer overall survival (OS) in OC patients ( ). Similarly, it could also be used as a predictive biomarker in TCGA-OV cohort ( ). Our study defined the CAF abundances calculated by EPIC and MCP-counter as phenotypic data for subsequent WGCNA. The data calculated by other algorithms were used for validation.

Figure 1

Normalization process based on the GPL570 platform dataset. (A) UMAP plot of the six datasets before normalization. (B) UMAP plot of the six datasets after normalization. (C) Expression density plot of the six datasets before normalization. (D) Expression density plot of the six datasets after normalization. (E) Expression distribution plots for the six datasets after normalization. UMAP, Uniform Manifold Approximation and Projection.

Figure 2

CAFs spatially associate with Treg cells and survival analysis-based CAFs and stromal score. (A) Immunofluorescence staining of the original ovarian tissue samples. (B) Immunofluorescence staining results of recurrent ovarian tissue samples. (C) The Kaplan–Meier analysis of GPL570 meta-cohort, including CAF_EPIC (i), Fibroblasts_MCPcounter (ii), Fibroblasts_xCell (iii), and StromalScore_estimate (iv). (D) The Kaplan–Meier analysis of TCGA-OV cohort, including CAF_EPIC (i), Fibroblasts_MCPcounter (ii), Fibroblasts_xCell (iii), and StromalScore_estimate (iv). CAFs, Cancer-Associated Fibroblasts; Tregs, Regulatory T cells.

Co-Expression Network of Cancer-Associated Fibroblast Scores

WGCNA was performed using the expression profiles of the top 25% of variance in the GPL570 meta-cohort and TCGA-OV cohort. The soft threshold power in the GPL570 meta-cohort was 3 ( ); similarly, the threshold for TCGA-OV cohort was also 3 ( ). Subsequently, dynamic module identification was performed in the different cohorts, with the number of genes per module not less than 50 ( ). For the GPL570 meta-cohort, 9 co-expression modules were clustered, with the brown module having the strongest positive correlation with CAFs_EPIC score (Cor = 0.88, P = 3e-208) and Fibroblasts_MCPcounter score (Cor = 0.9, P = 5e-234) ( ). For TCGA-OV cohort, the 9 co-expression modules were clustered, with the blue module having the strongest positive correlation with CAFs_EPIC score (Cor = 0.76, P = 2e-71) and Fibroblasts_MCPcounter score (Cor = 0.92, P = 3e-157) ( ). In the brown module, positive correlations between CAFs_EPIC score (Cor = 0.96) and Fibroblasts_MCPcounter score (Cor = 0.97) were observed between MM and GS ( ); in the black module, positive correlations between CAFs_EPIC score (Cor = 0.87) and Fibroblasts_MCPcounter score (Cor = 0.97) were also observed between MM and GS ( ). Finally, 120 genes in the brown module and 160 genes in the blue module were screened as potential CAF-related genes using MM > 0.6 and GS > 0.6 as thresholds.

Figure 3

WGCNA in the GPL570 meta-cohort and TCGA-OV cohort. (A) Scale independence and mean connectivity in the GPL570 meta-cohort. (B) Scale independence and mean connectivity in TCGA-OV cohort. (C) Gene dendrogram and modules before merging in the GPL570 meta-cohort. (D) Gene dendrogram and modules before merging in TCGA-OV cohort. (E) Pearson correlation analysis of merged modules and CAF score in the GPL570 meta-cohort. (F) Pearson correlation analysis of merged modules and CAF score in TCGA-OV cohort. (G) Scatterplot of MM and GS from the brown module in the GPL570 meta-cohort, including CAFs_EPIC (i) and Fibroblasts_MCPcounter (ii). (H) Scatterplot of MM and GS from the blue module in TCGA-OV cohort, including CAFs_EPIC (i) and Fibroblasts_MCPcounter (ii). WGCNA, Weighted Gene Co-expression Network Analysis; CAFs, Cancer-Associated Fibroblasts; GS, Gene Significance; MM, Module Membership.

Functional Analyses of Cancer-Associated Fibroblast-Related Genes

The above CAF-related genes were overlapped and screened to 95 hub genes ( ). Regulation of small GTPase-mediated signal transduction, extracellular matrix, collagen-containing extracellular matrix, and metallopeptidase activity were the main enriched GO terms ( ). Fatty acid degradation, glycolysis/gluconeogenesis, regulation of lipolysis in adipocytes, peroxisome proliferator activated receptor (PPAR) signaling pathway, vascular smooth muscle contraction, and cyclic guanosine monophosphate (cGMP)/protein kinase G (PKG) signaling pathway were the mainly enriched KEGG pathways ( ).

Figure 4

Functional analyses and construction of the CAF-based signature. (A) The hub CAF-related genes were overlapped in the brown module and the blue module. (B) GO enrichment analysis. (C) KEGG enrichment analysis. (D) Univariate Cox regression analysis of common hub genes (P < 0.001). (E) LASSO regression analysis. (F) Kaplan–Meier analysis of different cohorts. On the left is GPL570 meta-cohort; on the right is TCGA-OV cohort. CAFs, Cancer-Associated Fibroblasts; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, Least Absolute Shrinkage and Selection Operator.

Construction of the Cancer-Associated Fibroblast-Based Signature

The GPL570 meta-cohort with a larger sample size was used as the training cohort, and TCGA-OV cohort was used as the validation group. Univariate Cox regression analysis was performed on common hub genes in the training cohort ( ), with OS and survival time as dependent variables, and 63 prognostic genes (P < 0.05) were screened, and only some with P < 0.001 were shown in . The 63 prognostic genes were subjected to LASSO regression analysis to determine the minimum λ value ( ). Finally, 6 genes were identified for the CAF-based signature: CAF risk score = COL16A1 expression * 0.0924 + COL5A2 expression * -0.0031 + GREM1 expression * 0.0847 + LUM expression * 0.0069 + SRPX expression * 0.0649 + TIMP3 expression * 0.0425. The OC patients in each cohort were divided into high-risk and low-risk groups, and the cutoff for each cohort was used as the threshold value (GPL570 meta-cohort = 1.257016302; TCGA-OV cohort = 0.415034301). Kaplan–Meier curves showed that patients in the high-risk group had worse OS than that of those in the low-risk group ( ). These results suggested that the CAF signature was the hub prognostic marker for OC patients.

Cancer-Associated Fibroblast-Based Signature Genes Were Correlated With Cancer-Associated Fibroblast Markers

Spearman correlation analyses were performed between the CAF risk score and the CAF score predicted by the other methods (xCell, EPIC, ESTIMATE, and MCP-counter). Subsequently, we observed a strong and positive correlation between risk scores and CAF infiltration and stromal score in both GPL570 meta-cohort ( ) and TCGA-OV cohort ( ). Moreover, CAF marker genes from previous references had a higher expression in the high-risk group ( ). In addition, the expression levels of 6 genes in the signature also were highly and positively correlated with CAF marker expression ( ).

Figure 5

Genes involved in the signature were correlated with CAF markers. (A) Correlation analysis of CAF score, stromal score, and risk score in the GPL570 meta-cohort. (B) Correlation analysis of CAF score, stromal score, and risk score in TCGA-OV cohort. (C) Heatmaps of expression of CAF markers in different risk groups (GPL570 meta-cohort). (D) Heatmaps of expression of CAF markers in different risk groups (TCGA-OV cohort). (E) Correlation analysis of genes involved in signature and CAF markers (GPL570 meta-cohort). (F) Correlation analysis of genes involved in signature and CAF markers (TCGA-OV cohort). CAFs, Cancer-Associated Fibroblasts.

Multidimensional Validation in Multicenter Studies

To further validate the prognostic value of the CAF-based signature, we integrated the GPL96 meta-cohort (GSE3149, GSE23554, GSE26712, and GSE14764) according to the method described above, which included a total of 409 patients ( ). Meanwhile, the datasets based on GPL7759 (GSE13876, n = 415), GPL2986 (GSE49997, n = 194), and GPL14951 (GSE140082, n = 380) platforms were downloaded for external validation. The risk scores of each cohort were calculated with the same formula and stratified by their respective cutoff values (GPL96 meta-cohort = 0.976888643; GPL7759 cohort = 3.088372669; GPL2986 cohort = 0.147731773; GPL14951 cohort = 2.479072527). Unsurprisingly, risk score stratified patients by survival risk in multicenter studies, and OS was shorter in the high-risk group, such as in the GPL96 meta-cohort ( , P = 0.004), GPL7759 cohort ( , P = 0.006), GPL2986 cohort ( , P < 0.001), and GPL14951 cohort ( , P = 0.002).

Figure 6

Multidimensional validation for risk score. (A) Kaplan–Meier analysis in the GPL96 meta-cohort. (B) Kaplan–Meier analysis in the GPL7759 cohort. (C) Kaplan–Meier analysis in the GPL2986 cohort. (D) Kaplan–Meier analysis in the GPL14951 cohort. (E) Histogram of anti-PD-L1 response distribution in different risk groups. (F) Box plot of risk score in different anti-PD-L1 response groups. (G) Kaplan–Meier analysis in the anti-PD-L1 cohort. (H) Histogram of anti-PD-1 response distribution in different risk groups. (I) Box plot of risk score in different anti-PD-1 response groups. (J) Kaplan–Meier analysis in the anti-PD-1 cohort. PD-1, Programmed Death-1; PD-L1, Programmed Cell Death 1 Ligand.

Cancer-Associated Fibroblast-Based Signature in the Role of Immunotherapy

Immunotherapy represented by PD-L1 and PD-1 blockade has undoubtedly become a breakthrough in cancer treatment, so we investigated whether the CAF-based signature could predict response to anti-PD-1 and anti-PD-L1 based on two immunotherapy cohorts. In the anti-PD-L1 cohort (IMvigor210), the high-risk group had a higher percentage of stable disease (SD)/progressive disease (PD). In contrast, more patients in the low-risk group were in complete response (CR)/partial response (PR) ( ). Moreover, patients with a low risk score exhibited a markedly prolonged survival ( ). In the anti-PD-1 cohort (GSE78220), the significant therapeutic advantages and clinical response in patients with a low score also were confirmed ( ). However, due to the small sample size of the anti-PD-1 cohort, there was no significant difference in survival time between different groups ( ).

Correlation Between the Cancer-Associated Fibroblast-Based Signature and Somatic Variation

Preclinical research has shown that patients with higher TMB are associated with enhanced immunotherapy response and lasting clinical benefits when treated with immune checkpoint blockade. Therefore, we investigated the discriminatory ability of the CAF-based signature in the somatic mutation data of TCGA-OV cohort. Firstly, we screened the most differentially mutated genes in different risk groups, including KMT2C, WDFY3, CACNA1S, etc. ( ). We found no significant differences between the two groups in CAF marker mutations, but TNC (15.0%) and COL3A1 (11.7%) exhibited a higher frequency of mutations in the whole TCGA-OC cohort ( ). Subsequently, we observed that TMB values were higher in the low-risk group than those in the high-risk group ( ). However, Spearman analysis showed no statistically significant correlation between CAF risk score and TMB values ( ). However, TMB values were negatively correlated with stromal score and CAF-activating factors transforming growth factor beta (TGF-β), suggesting that higher TMB might have intense tumor-killing effects via modulating a fibroblast-weak TME (40) ( ).

Figure 7

Functional analyses and construction of the CAF-based signature. (A) (i) Differentially mutated genes in different risk groups. (ii) CAF marker mutations in different risk groups. (B) TMB values in different risk groups. (C) Spearman analysis between CAF risk score and TMB values. (D) Correlation analysis between TMB values, stromal/CAF scores, and CAF-activating factors. CAFs, Cancer-Associated Fibroblasts; TMB, Tumor Mutational Burden.

GSEA of the Cancer-Associated Fibroblast-Based Signature

Gene Set Enrichment Analysis (GSEA) was performed in two datasets (GPL570 meta-cohort and TCGA-OV cohort) to explore the pathways involved in different risk groups. Allograft rejection, apical junction, and epithelial–mesenchymal transition were significantly enriched ( ). The ssGSEA score also showed that the CAF risk score was positively correlated with TNFA signaling via nuclear factor-kappaB (NF-kappaB), hypoxia, and Wnt beta catenin signaling pathway ( ).

Figure 8

GSEA in cancer hallmark gene set. (A) GSEA plot in the GPL570 meta-cohort. (B) GSEA plot in TCGA-OV cohort. (C) Correlation analysis between risk score and TNF signaling (i), hypoxia (ii), EMT (iii) in the GPL570 meta-cohort. (D) Correlation analysis between risk score and TNF signaling (i), hypoxia (ii), and EMT (iii) in TCGA-OV cohort. GSEA, Gene Set Enrichment Analysis; TNF, Tumor Necrosis Factor; EMT, Epithelial to Mesenchymal Transition.

Sensitivity of Chemotherapy Between Different Risk Groups

Maintenance therapy and chemotherapy after debulking surgery for OC patients are crucial, and the mutation of the Breast Cancer Susceptibility Genes (BRCA) is relevant to the efficacy of olaparib. Therefore, we explored the distribution of mutations in the BRCA under different risk groups. BRCA1 may be more distributed in the high-risk group, but there was no significant difference in BRCA2. Interestingly, the combined BRCA mutation status and risk score allowed for better survival prediction ( ). In addition, Wilcoxon analysis revealed significant differences in IC50 values between different risk groups. Among them, high-risk patients may be more sensitive to bleomycin ( ), cisplatin ( ), docetaxel ( ), and gemcitabine ( ). Still, the IC50 values of doxorubicin ( ) and etoposide ( ) were not significantly different between groups.

Figure 9

Sensitivity of chemotherapy between different risk groups. (A) Histogram of BRCA1 state distribution and Kaplan–Meier analysis of integrated groupings. (B) Histogram of BRCA2 state distribution and Kaplan–Meier analysis of integrated groupings. IC50 values between different risk groups, including bleomycin (C), cisplatin (D), docetaxel (E), doxorubicin (F), etoposide (G), and gemcitabine (H).

Validation in Cell Lines, scRNA-Seq, and Immunohistochemistry

To validate that the CAF-related genes involved in the signature were the primary origins in CAFs, we performed a multidimensional validation, including cell lines, single-cell sequencing, and immunohistochemistry. We collected cell line RNA-seq data from 47 fibroblast origins and 37 OC origins. We found that all six genes (COL16A1, COL5A2, GREM1, LUM, SRPX, and TIMP3) were overexpressed in fibroblasts by the “limma” package ( ) and Wilcoxon test ( ). Meanwhile, we annotated the scRNA-seq into 6 clusters: fibroblasts, myofibroblasts, endothelial, malignant, Mono/Macro, and CD4Tconv ( ). The differential analysis results showed that most CAF-related genes were highly expressed in fibroblasts or myofibroblasts, while lower expression was observed in malignant ( ). Moreover, the single-cell GSEA was consistent with the bulk-RNA GSEA, showing significant enrichment of upregulated genes of fibroblasts in the EMT pathway ( ). We analyzed IHC images from the HPA database, and the section showed that GREM1 and LUM proteins were deeply stained in the stroma ( ). Unfortunately, the other four genes did not have corresponding IHC images in the HPA database. These verifications implied that these six genes might be CAF-specific markers.

Figure 10

Multidimensional expression validation. (A) Heatmap of gene expression in different cell lines based on “limma” package. (B) Wilcoxon test of gene expression in different cell lines. (C) Major cell type in single-cell seq. (D) Differential distribution of gene expression at the single-cell level. (E) GSEA of upregulated genes in different cell types. (F) IHC images of OC tissues from the HPA database. GSEA, Gene Set Enrichment Analysis; OC, Ovarian cancer; IHC, Immunohistochemistry; HPA, Human Protein Atlas.

Discussion

The CAF is regarded as an essential factor in promoting tumor progression by interacting with cancer cells in the TME (1). Meanwhile, for a specific mesenchymal subtype of OC, it is characterized by frequent generation of desmoplastic stroma (41). The generation of desmoplastic stroma is associated with a lower OS and resistance to platinum (42). Consistently, we observed that higher CAF and stromal scores were associated with poorer OS in OC patients and represented a poorer immunotherapy response. This is the first study with a large sample and using WGCNA as a starting point for exploring markers associated with CAFs. A 6-gene prognostic (COL16A1, COL5A2, GREM1, LUM, SRPX, and TIMP3) signature was constructed and validated using Cox and LASSO regression algorithms. With the cutoff value as a threshold, we observed that patients with a high CAF risk score were more sensitive to numerous chemotherapeutic agents. Furthermore, we revealed that lower risk scores were associated with improved immunotherapy outcomes and higher TMB value. Based on our results, we propose an alternative mechanism by which higher TMB may also enhance tumor killing by modulating the microenvironment of stromal fibroblasts, similar to previous findings. It reported that cancer cells with high levels of somatic mutation are more easily recognized by the immune system (43). However, we need more in vitro and in vivo experiments to elucidate the above crosstalk in the future.

Compared to the traditional differential gene expression (DEG) approach for screening hub CAF markers (44), we used different bioinformatics algorithms to assess the abundance of CAFs and biomarkers in each OC sample to ensure the robustness (EPIC and MCP-counter for WGCNA network construction; xCell and ESTIMATE for correlation validation). Similarly, to ensure the robustness of the prognostic signature, different cohorts were used for construction and validation (GPL570 meta-cohort for construction; TCGA-OV cohort, GPL96 meta-cohort, GPL7759 cohort, GPL2986 cohort, GPL14951 cohort, IMvigor210 cohort, and GSE78220 cohort for validation). With the above approach, we confirmed that our model closely correlated with CAF infiltration and CAF markers from references. Meanwhile, to differentiate identified genes from tumor cells to highlight gene heterogeneity in CAFs, we confirmed significantly higher expression in fibroblast cell lines, higher staining of proteins in the stroma, and higher mRNA expression of the CAFs at the single-cell level.

For the six genes involved in the risk signature, the relevant references have reported on the role in tumor cells and TME. COL16A1 was indicated in the study of Pan and Ma (45) to be involved in a risk model and could be considered a prognostic marker in OC patients. Renner et al. (46) sought to determine the ECM composition of benign fallopian tubes and the changes associated with tubal intraepithelial carcinomas and identified seven proteins that had not been identified in previous studies (COL2A1, COL4A5, COL16A1, elastin, LAMA5, annexin A2, and PAI1). Interestingly, they suggested that the seven proteins mentioned above accompany tubal intraepithelial carcinoma formation and cause ECM changes. Head and neck squamous cell carcinoma (HNSCC) cell lines were cocultured with their patient-matched CAFs in 2D and 3D in vitro models, and GREM1 was upregulated (47). In addition, related studies have also found that GREM1 binds to miR-205-5p (48) or miR-206 (49) to regulate metastasis of cervical cancer and non-small cell carcinoma. As part of the ECM, collagen family proteins, together with elastins, fibronectins, and laminins, play a key role in tissue organization, tissue resistance, and its primary shape. The collagen family, including COL5A2, is overexpressed in various types of epithelial cancers and is associated with poorer OS. Furthermore, inhibition of gene expression decreases cell proliferation and invasion (50, 51). SRPX, also known as SRPX1 (52), ETX1 (53), and DRS (54), is a suppressor that has been found to be downregulated in a range of human tumor cells and tissues. Unlike other soluble members of the Tissue Inhibition of Matrix Metalloproteinase (TIMP) family, TIMP3 is tightly sequestered in the ECM. TIMP-3 is also the only TIMP capable of inhibiting tumor necrosis factor alpha (TNF-α), ADAMTS4, and ADAMTS5, as well as syndecan sheddase (55). Nevertheless, functional validation about the six genes involved in the risk signature in the CAFs of OC is not much, which requires us to conduct further experiments on the six CAF markers in the future.

In conclusion, the CAF risk score can be used in clinical practice to comprehensively evaluate the corresponding cellular infiltration of CAFs in patients to further define the immunophenotype. We have also demonstrated that risk score can be used to assess the clinicopathological characteristics of patients. Similarly, risk score also can be used as a biomarker to predict survival and the efficacy of adjuvant chemotherapy and the response to anti-PD-1/PD-L1 immunotherapy. More importantly, this study may help to leverage the future development of new drug combination strategies or new immunotherapeutic agents. Our findings provide new ideas to facilitate future individualized cancer immunotherapy.

Data Availability Statement

The original contributions presented in the study are included in the article/ . Further inquiries can be directed to the corresponding author.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics Committees and Institutional Review Boards of Zhongda Hospital Southeast University (ZDSYLL187-P04). The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

Author Contributions

SF, YX, and ZD conceived and designed the study. SY was responsible for materials. SF drafted the article. SY, HY, and KZ revised the article critically. All authors had final approval of the submitted versions.

Funding

This study was supported by National Natural Science Foundation of China (No. 82072078), Jiangsu Province Key Research and Development Project (SBE2020741118), and Postgraduate Research & Practice Innovation Program of Jiangsu Province (SJCX22_0070).

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.