A Novel Molecular Signature of Cancer-Associated Fibroblasts Predicts Prognosis and Immunotherapy Response in Pancreatic Cancer

Cancer-associated fibroblasts (CAFs), a prominent population of stromal cells, play a crucial role in tumor progression, prognosis, and treatment response. However, the relationship among CAF-based molecular signatures, clinical outcomes, and tumor microenvironment infiltration remains largely elusive in pancreatic cancer (PC). Here, we collected multicenter PC data and performed integrated analysis to investigate the role of CAF-related genes (CRGs) in PC. Firstly, we demonstrated that α-SMA+ CAFs were the most prominent stromal components and correlated with the poor survival rates of PC patients in our tissue microarrays. Then, we discriminated two diverse molecular subtypes (CAF clusters A and B) and revealed the significant differences in the tumor immune microenvironment (TME), four reported CAF subpopulations, clinical characteristics, and prognosis in PC samples. Furthermore, we analyzed their association with the immunotherapy response of PC patients. Lastly, a CRG score was constructed to predict prognosis, immunotherapy responses, and chemosensitivity in pancreatic cancer patients. In summary, these findings provide insights into further research targeting CAFs and their TME, and they pave a new road for the prognosis evaluation and individualized treatment of PC patients.


Introduction
The poor prognosis of pancreatic cancer (PC) urges us to more deeply understand its potential molecular mechanism and seek better therapies. Cancer-associated fibroblasts (CAFs), a significant fraction of the pancreatic cancer stroma, contribute to a dense stromal accumulation in PC [1,2]. Previous studies have demonstrated that CAFs can facilitate the malignant phenotypes of tumors, particularly tumorigenesis and invasion, inflammation, and extracellular matrix (ECM) remodeling [3,4]. As the leading participant of the desmoplastic stroma in PC, CAFs play a crucial part in diverse clinical responses, drug tolerance, and the tumor immunosuppressive environment by producing ECM proteins and cytokines and interacting with cancer cells [5][6][7][8]. The intratumoral heterogeneity of CAFs in the stroma of PC has been extensively studied. Four major distinct subpopulations of CAFs have been demonstrated in PC: (1) myCAF, the myofibroblastic subset (myCAF) characterized by smooth muscle actin expression, high transforming growth factor (TGF) signaling, and ECM components [1,9]; (2) iCAF, the inflammatory subset characterized by high expressions of inflammatory mediators [10]; (3) apcCAF, the antigen-presenting subset characterized by the expression of CD74 and MHC class II [11]; and (4) meCAF, the highly activated metabolic subset characterized by high expression of PLA2G2A and CRABP2 [12]. Some CAF markers have been studied separately in the past few years, such as α-smooth muscle actin (α-SMA), fibroblast activation protein, CD29, fibroblast-specific protein 1, platelet-derived growth factor receptor B, and podoplanin [13]. For instance, tumors accumulated with the α-SMA + fibroblasts have worse prognoses and higher invasiveness, and they can affect therapeutic reactions [13,14]. Fibroblast activation protein-positive CAFs can lead to immunosuppression and resistance to immunotherapy [15]. However, whether the CAF-mediated tumor microenvironment (TME) is associated with tumor characteristics and the underlying molecular mechanism remains unclear [16,17]. In addition, the guiding significance of current pathological and molecular classification for PC treatment is limited [18]. Although existing strategies targeting the stroma have suppressed tumor growth and enhanced treatment responses in the mouse model, clinical trials have not yet produced promising results [6,19]. Some of these strategies even lead to tumor recurrence and metastasis [20], which suggests that accurate identification of the CAF molecular subtypes in PC is necessary to apply stromal-targeting therapies efficiently.
This study determined that PC patients with accumulated α-SMA + CAFs had a poor prognosis regarding tissue microarrays, and that myCAF, apcCAF, and meCAF subsets were highly enriched in PC. Then, we thoroughly estimated the expression profiles of CAF-related genes (CRGs) and their influences on prognosis, clinical features, and immune cell infiltration in PC patients. Furthermore, we constructed a CRG score to predict PC patients' prognoses, clinical outcomes, immunotherapy responses, and chemosensitivity. Our findings could deepen our understanding of CRGs and smooth the way for prognosis evaluation and personalized therapy strategies in PC patients.

α-SMA + CAFs Accumulate in PC Tissues with Worse Prognoses
Previous studies have shown that α-SMA is a marker of activated CAFs and an efficacy evaluation indicator of targeted CAF therapy [21,22]. The Kaplan-Meier curves showed worse overall survival (OS) in patients with high α-SMA + CAFs accumulated in PC ( Figure 1F). Immunofluorescence staining, immunohistochemistry, and Masson staining were used to further confirm the CAF population in PC tissue microarrays ( Figure 1A-E). The results explain that α-SMA + CAFs, as a prominent desmoplastic stroma, were remarkably enriched in PC tissues. To define the substantial proportion of CAFs in human PC tissues and mousederived allografts, we excluded hematopoietic and epithelial cells using flow cytometry analysis with CD45 and EpCAM markers. We also identified CAFs using human fibroblast markers (integrin b1/CD29) and mouse fibroblast marker (PDPN). The fresh human samples included PC patients at the time of surgery before any treatment. Figure 1G-J illustrate that the CAFs accounted for about 30% of all cellular populations in human and mouse tumor tissues. We further assessed the enrichment scores of four reported CAF subtypes in pancreatic tumors and normal pancreas samples from TCGA-PAAD and GTEx cohorts, and we found that myCAF, apcCAF, and meCAF were abundant in the tumor samples ( Figure 2A). These results suggest that CAFs are essential components of TME in PC, which may modulate tumorigenesis and progression.

Genetic Mutation Landscape of CRGs in PC
We first determined the expression levels of the 25 CRGs in tumor specimens and normal specimens, and we observed that almost all CRGs were abundant in the tumor specimens ( Figure 2B). To reveal the interaction of CRGs, we performed a PPI analysis. Figure 2C displayed that COL1A1, COL11A1, COL3A1, COL5A2, COL1A2, FN1, FAP, CDH1, POSTN, COMP, COL5A1, COL10A1, and THBS2 were hub genes. Furthermore, we identified the total frequency of somatic mutations and copy number variations (CNVs) of the 25 CRGs in PC. As depicted in Figure 2D, 16 of 158 (10.13%) PC samples emerged with genetic mutations. Figure 2D also indicated that, among the 25 CRGs, VCAN, FN1, and COL11A1 were the genes with the highest mutation rate, followed by COL5A1 and CDH1. In addition, we demonstrated evident CNV alterations of the 25 CRGs ( Figure 2E). We also analyzed the CNV alteration location of the 25 CRGs on chromosomes using the "circlize" R package ( Figure 2F). We concluded that CNV might act in regulating the expression of 25 CRGs. These findings reveal significant differences in the genomic background and expression levels of the 25 CRGs between PC and normal samples, implying the latent roles of the 25 CRGs in PC tumor progression.

Genetic Mutation Landscape of CRGs in PC
We first determined the expression levels of the 25 CRGs in tumor specimens and normal specimens, and we observed that almost all CRGs were abundant in the tumor specimens ( Figure 2B). To reveal the interaction of CRGs, we performed a PPI analysis. Figure 2C displayed that COL1A1, COL11A1, COL3A1, COL5A2, COL1A2, FN1, FAP, CDH1, POSTN, COMP, COL5A1, COL10A1, and THBS2 were hub genes. Furthermore, we identified the total frequency of somatic mutations and copy number variations (CNVs) of the 25 CRGs in PC. As depicted in Figure 2D, 16 of 158 (10.13%) PC samples emerged with genetic mutations. Figure 2D also indicated that, among the 25 CRGs, VCAN, FN1, and COL11A1 were the genes with the highest mutation rate, followed by COL5A1 and CDH1. In addition, we demonstrated evident CNV alterations of the 25

Identification of CAF Subtypes and Characteristics of the TME in PC
To better understand the expression pattern of CRGs in tumorigenesis, we performed a subsequent analysis of 160 PC patients from TCGA-PAAD. Table S2 lists detailed information about these patients. We further performed a consensus clustering analysis to investigate the relationships between the expression pattern of CRGs and PC subtypes, and we classified PC patients according to the expression levels of these CRGs. Our findings indicate that k = 2 is an optimal choice to divide the entire cohort into CAF cluster A (n = 130) and CAF cluster B (n = 30) ( Figures 3A and S1). Moreover, we used the ICGC cohort to verify the repeatability of the clustering. We also conducted a consensus clustering analysis on this cohort and classified the cohort into two distinct subtypes ( Figure S2A,B). Patients with CAF cluster A had worse OS than patients with CAF clus-ter B in both TCGA and ICGC cohorts ( Figures 3B and S2C). We further dissected the CAF signature of the patients in two CAF subtypes. The expression of the CAF signature in CAF cluster A was substantially higher than in cluster B ( Figure 3C). Figure 3D presents the relevant networks of CRG interactions and regulator connections. It also illustrates the prognostic value of CRGs (Table S3) and the enrichment of the CRG-related KEGG pathways (Table S4) in PC patients. Additionally, significant differences in the genomic expression of CRGs and clinical variables were observed between the two CAF clusters ( Figure 3E). Apart from the differences in prognosis and genome between CAF cluster A and CAF cluster B, there were also distinct discrepancies in immune cell infiltration and TME score between them. Firstly, we observed higher enrichment scores of myCAF and apcCAF in the CAF cluster A group ( Figure 3F). To investigate the roles of CRGs in the TME of PC, Correlations between the two CAF clusters and TME score. (H) The infiltration abundance of 33 TME cells of two CAF subtypes in PC. The Wilcoxon test analyzed the statistical differences between the two clusters (*** p < 0.001, ** p < 0.01, * p < 0.05, and not significant (p > 0.05)).
Apart from the differences in prognosis and genome between CAF cluster A and CAF cluster B, there were also distinct discrepancies in immune cell infiltration and TME score between them. Firstly, we observed higher enrichment scores of myCAF and apcCAF in the CAF cluster A group ( Figure 3F). To investigate the roles of CRGs in the TME of PC, we then evaluated the association among the two CAF clusters, 33 immune cell subtypes, and the TME score (Table S5). Compared with CAF cluster B, CAF cluster A had higher immune and stromal scores ( Figure 3G) and higher infiltration levels of immunosuppressive cells, such as regulatory T cells (Tregs), MDSC cells, and DC cells, and other immunosuppressive factors, such as TGF-β-associated ECM ( Figure 3H). More importantly, we detected a higher enrichment score of anti-PD-1-resistant signatures and a lower enrichment score of nivolumab-responsive signatures in CAF cluster A ( Figure 3H), indicating that patients in the CAF cluster A group may be less sensitive to immunotherapy. These results imply that the CAF cluster A group may be closely associated with stromal activation and immunosuppression features.

Establishment and Verification of the Prognostic CRG Score
The CRG score was created according to the LASSO and multivariate Cox (multiCox) analysis for 25 CRGs. Eventually, we obtained five hub genes (VCAN, COL1A2, ZNF469, SPARC, and FNDC1). The CRG score was calculated as follows: Figure 4A displays the distribution of patients in the two CAF clusters and two CRGscore groups. Compared with alive patients, the CRG score was significantly elevated in patients who died during follow-up ( Figure 4B), and CAF subtype A had higher CRG scores ( Figure 4C). The risk plot of the CRG score indicated that, with an increasing CRG score, OS time decreased while mortality rose ( Figure 4D,G). Patients with higher CRG scores in both categories were associated with worse survival rates ( Figure 4E). Additionally, the AUC values of 1-, 2-, and 3-year OS were 0.63, 0.659, and 0.638, respectively ( Figure 4F). Moreover, the CRG score retained excellent predictability in assessing the prognosis of PC patients ( Figure 4H). Among multiple clinical features, multivariate Cox regression modeling proved that the CRG score was the only independent risk factor for the OS of PC patients in the TCGA cohort ( Figure 4I).

Characteristics of the TME and Function Enrichment in Distinct Subgroups
To examine the association between CRG score and the TME of PC, we analyzed their immune microenvironment in detail. As confirmed by different methods, the CRG score was positively associated with M1 macrophages and neutrophils, whereas it was negatively related to B cells, NK cells, CD8 T cells, and CD4 T cells ( Figure 5A). Moreover, we sought to explore the potential pathways related to the CRG score using GSVA. Several cancer-associated pathways (P53, Notch, and ERBB pathways) were most closely correlated with the CRG score ( Figure 5B). Consistently, we found that the enrichment levels of B cells, plasma cells, CD8 T cells, and CD4 T cells were markedly higher in the low-CRGscore group than in the high-CRG-score group ( Figure 5C). Figure 5D reveals a higher enrichment score of meCAF in the low-CRG-score group. Furthermore, time-dependent receiver operating characteristic (tROC) analysis showed that the CRG score was the most accurate predictor for overall survival compared with other single-CAF subsets in PC ( Figure 5E). These findings indicate that patients with lower CRG scores had higher meCAF accumulation and more immune cell infiltration.

Characteristics of the TME and Function Enrichment in Distinct Subgroups
To examine the association between CRG score and the TME of PC, we analyzed their immune microenvironment in detail. As confirmed by different methods, the CRG score was positively associated with M1 macrophages and neutrophils, whereas it was negatively related to B cells, NK cells, CD8 T cells, and CD4 T cells ( Figure 5A). Moreover, we sought to explore the potential pathways related to the CRG score using GSVA. Several cancer-associated pathways (P53, Notch, and ERBB pathways) were most closely correlated with the CRG score ( Figure 5B). Consistently, we found that the enrichment levels of B cells, plasma cells, CD8 T cells, and CD4 T cells were markedly higher in the low-CRG-score group than in the high-CRG-score group ( Figure 5C). Figure 5D reveals a

Association of the CRG Score with Tumor Mutation Burden (TMB) and Mutation
Previous studies have indicated that TMB is a valuable predictor of survival outcomes and immunotherapy response in tumor patients [23]. We explored the distribution alternations of somatic mutations between two CRG-score groups in TCGA cohort ( Figure 6A,B). Patients with high CRG scores had substantially higher frequencies of TP53, KRAS, CDKN2, SMAD4, and TTN mutations than patients with low CRG scores, implying that these gene mutations were in charge of the poor prognosis of PC patients with high CRG scores. However, we observed opposite results regarding the mutation levels of RNF43, MUC16, and RYR1 ( Figure 6A,B). In addition, our analysis of the mutation data demonstrated a higher TMB score in the high-CRG-score group compared with the low-CRG-score group ( Figure 6C).
higher enrichment score of meCAF in the low-CRG-score group. Furthermore, time-dependent receiver operating characteristic (tROC) analysis showed that the CRG score was the most accurate predictor for overall survival compared with other single-CAF subsets in PC ( Figure 5E). These findings indicate that patients with lower CRG scores had higher meCAF accumulation and more immune cell infiltration. The enrichment score of four CAF subsets between the low-CRG-score group and the high-CRG-score group. (E) tROC analysis showed that the GRC score was an accurate variable for survival prediction. The Wilcoxon test analyzed the statistical differences between the two clusters (** p < 0.01 and not significant (* p > 0.05)).

Association of the CRG Score with Tumor Mutation Burden (TMB) and Mutation
Previous studies have indicated that TMB is a valuable predictor of survival outcomes and immunotherapy response in tumor patients [23]. We explored the distribution alternations of somatic mutations between two CRG-score groups in TCGA cohort ( Figure  6A,B). Patients with high CRG scores had substantially higher frequencies of TP53, KRAS, CDKN2, SMAD4, and TTN mutations than patients with low CRG scores, implying that these gene mutations were in charge of the poor prognosis of PC patients with high CRG scores. However, we observed opposite results regarding the mutation levels of RNF43, MUC16, and RYR1 ( Figure 6A,B). In addition, our analysis of the mutation data demonstrated a higher TMB score in the high-CRG-score group compared with the low-CRGscore group ( Figure 6C). The enrichment score of four CAF subsets between the low-CRG-score group and the high-CRG-score group. (E) tROC analysis showed that the GRC score was an accurate variable for survival prediction. The Wilcoxon test analyzed the statistical differences between the two clusters (** p < 0.01 and not significant (* p > 0.05)).

Clinical Outcomes and Drug Susceptibility Analysis
We investigated the CRG score's ability to predict the impact of initial surgical treatment in PC patients. As displayed in Figure 6D,E, among the patients receiving the initial therapy of surgery, those with lower CRG scores showed significant treatment advantages.
Subsequently, to explore the efficacy of the CRG score as a biomarker for predicting chemotherapeutic susceptibility in PC patients, we assessed the semi-inhibitory concentration of 138 chemotherapeutic drugs commonly used to treat tumors. We identified 27 drugs more sensitive to patients with low CRG scores (Table S6), including EHT.1864 and PD.173074 (p < 0.01; Figure 6F,G). Nevertheless, 15 drugs responded better to patients with high CRG scores (Table S7), including paclitaxel and lapatinib (p < 0.01; Figure 6H,I). In brief, these findings suggest that the CRG score is associated with drug sensitivity.

Protein Expression Level of CAF-Related Risk Genes and Survival Analysis
To validate the tissue expression of risk CRGs in pancreatic normal and tumor tissues, we obtained immunohistochemical results from the Human Protein Atlas (HPA). Except for ZNF469, which is not available in the HPA database, consistent with the mRNA level in Figure 2B, protein expressions of VCAN ( Figure 7A

Clinical Outcomes and Drug Susceptibility Analysis
We investigated the CRG score's ability to predict the impact of initial surgical treatment in PC patients. As displayed in Figure 6D,E, among the patients receiving the initial

Discussion
Immune and stromal cells, the essential TME components, are associated with the clinical features and prognosis of PC [24][25][26]. Extensive stromal involvement is a crucial hallmark of PC, which makes it challenging to obtain accurate tumor-specific molecular information [24]. Early studies identified that CAFs, a substantial portion of the tumor microenvironment, drove tumorigenesis and treatment resistance [6,27,28]. Previous studies revealed how CAF patterns affect the characteristics of TME and the efficacy of immunotherapy in triple-negative breast cancer (TNBC) [29]. With the development of tumor immunology and molecular biology research, immunotherapies, such as immune checkpoint inhibitors, have become new treatments for various tumors [30,31]. Recently, anti-PD-1/PD-L1 therapy has led to outstanding achievements in many malignancies [32,33]. However, due to the dense extracellular matrix acting as a physical barrier, PC patients remain poorly responsive to PD-1 antibodies [32,34]. Moreover, single-cell analysis also revealed that TGF-β-myCAF subtypes are related to the resistance to immunotherapy in breast cancer [35]. Whether analyzing CAF molecular subtypes improves the clinical response of PC remains to be determined [25]. Despite several studies having identified various biomarkers and clinical factors to predict PC prognosis [24,36], the relationship among CAF-based molecular signature, clinical outcome, and tumor microenvironment infiltration remains largely elusive in PC. Here, our study found an abundance of myCAF, apcCAF, and meCAF in the tumor tissues of PC. We also identified the alterations in genomic backgrounds and expression levels of CAF-related genes based on TCGA, GTEx, and ICGC cohorts. Most of the expressions of CRGs were increased in PC tumor tissues and correlated with prognosis. The aggregation of gene mutations leads to carcinogenesis, and gene mutations in PC may significantly impact immunotherapy response [37]. Among 25 CRGs, VCAN, FN1, and COL11A had the highest mutational intensity. However, there are currently no reports that these mutations are associated with carcinogenesis or fibrosis.
Additionally, we divided PC patients into two CAF clusters and observed discrepant prognoses, clinical characteristics, and immune infiltrations between them. The interaction of CAFs and immunity is a critical feature of tumorigenesis, which can serve as a therapeutic target for PC. Diverse CAF subsets play distinct roles in tumor immunosuppression of breast cancer. Their effects are achieved by Tregs regulating the proliferation of effector T cells [8]. Our findings showed that CAF cluster A, with a high enrichment of myCAF and apcCAF, had significantly higher stromal and immune scores than CAF cluster B. Cluster A also had higher infiltration levels of immunosuppressive cells, such as Tregs, MDSC cells, and DC cells, and other immunosuppressive factors, such as TGF-βassociated-ECM. Previous studies have shown that myCAF is the main component of the ECM [10]. Furthermore, apcCAF potentially modulates the immune response in pancreatic tumors [11]. Our results imply that CAF cluster A group may be closely associated with stromal activation and immunosuppression features, and myCAF and apcCAF abundance may be the main factors underlying such an immunosuppressive microenvironment. More interestingly, the higher enrichment score of anti-PD-1-resistant signatures and the lower enrichment score of nivolumab-responsive signatures were also observed in the CAF cluster A group, indicating that patients in the CAF cluster A group may be less sensitive to immunotherapy. Furthermore, the CRG score was constructed to quantify CAF subtypes. The CAF subtype A with worse survival had higher CRG scores. Patients with higher CRG scores also had worse OS, implying that high CRG scores could predict an unfavorable prognosis. By integrating the CRG scores and clinical characteristics, we demonstrated that the CRG score was a unique, independent risk factor of OS. Moreover, we found a higher enrichment of meCAF in PC patients with low CRG scores. Our previous research showed that PDAC patients with abundant meCAF had a dramatically better response to immunotherapy [12]. Consistent with the previous conclusion, the CAF cluster B with better survival had a lower CRG score, and it also had a lower enrichment score of anti-PD-1-resistant signatures and a higher enrichment score of nivolumab-responsive signatures. These findings indicate that CRGs may participate in tumor immunosuppression; therefore, patients with low CRG scores can benefit from immunotherapy.
Due to PC patients' distinctive molecular and clinical features, it is necessary to classify them precisely. We further identified potentially sensitive drugs in patients in different CRG-score groups. We expected that targeting CAFs combined with these drugs could reduce drug resistance and improve clinical outcomes.

Human Tissue Specimens
The Human Ethics Committee of Shanghai Renji Hospital, Shanghai Jiao Tong University School of Medicine (Shanghai, China), reviewed and approved research on human pancreatic cancer under informed consent from all patients.

Cell Lines and Mouse Pancreatic Cancer Model
Mouse pancreatic cancer cell lines KPC1199 were obtained from Jing Xue lab (Shanghai, China) and cultured in DMEM with 10% FBS. A total of 1 × 10 6 KPC1199 cells were resuspended in 100 µL of PBS and injected subcutaneously into 6-week-old female C57BL/6 mice from Shanghai Laboratory Animal Center. The tumor tissues were ultimately weighed in 15 days and collected for flow cytometry analysis.

Masson's Trichrome Staining
Formalin-fixed tissues were immersed in paraffin, and 5 µM sections were stained with Masson trichrome reagent to show collagen. First, the samples were partially dewaxed and rehydrated, fixed in Bouin's liquor, and then washed and rinsed in distilled water overnight. Next, the slides were stained in Meyer hematoxylin solution for 5 min, and then placed in 0.5% hydrochloric acid and 70% ethanol for 5 s. After the specimens were washed three times and dissolved in 1% phosphomolybdic acid aqueous solution, the slides were stained with aniline blue or bright green for 5 min. Subsequently, we dehydrated the samples in 95% ethanol 10 times and added xylene to make them transparent. Finally, all slides were scanned and digitized using the digital pathological slice scanner system (Leica Biosystems Wetzlar, Germany). The collagen fibers were dyed blue, the nuclei were black, and the background was red.

Immunohistochemistry and Tissue Microarrays
Tissue microarrays included 91 PC samples. The histopathology of all cancer specimens was reassessed, and representative regions were labeled. Table S1 lists the prognosis information of individual patients in tissue microarrays. The immunohistochemistry staining was carried out to identify α-SMA and EpCAM marker expressions in tissue microarrays. As previously described, we conducted a semi-quantitative scoring system according to the distinct percentages of positively stained cells and staining intensity [38]. The frequency of positively stained cells was defined as 1+ (less than 25%), 2+ (25% to 50%), 3+ (50% to 75%), or 4+ (greater than 75%). Additionally, the intensity was scored as 0 (negative), 1+ (weak), 2+ (moderate), or 3+ (strong).

Flow Cytometry Analysis
Briefly, the fresh tissue samples from two human PC tissues and seven mouse-derived allografts of KPC1199 cells were mechanically chopped and digested by collagenase IV at 37 • C for 30 min. The digested suspension was combined with DNase at room temperature for 5 min, washed twice with phosphate-buffered saline buffer containing 2% serum, and then filtered through the 100 µm filter. We used markers of CD29, PDPN, EpCAM, and CD45 to separate tumor epithelial cells (CD45 − CD29 − or PDPN − EpCAM + ), CAFs (EpCAM − CD45 − CD29 − or PDPN + ), and tissue leukocytes (CD29 − or PDPN − EpCAM − CD45 + ) in human and mouse tumor specimens, respectively. The digested single cells were washed twice and centrifuged for 5 min at 500× g, and then 1 µg/mL of antibody was added. Then, the samples were kept at 4 • C for 30 min in a dark place. Flow cytometry was employed using a BD Flow Cytometry Analysis Celesta cell sorter (Becton Dickinson, New York, NY, USA). The sidescatter width versus side-scatter region and the forward-scatter width versus forward-scatter height were applied to remove dead cells and cell clumps. Antibodies including anti-EpCAM-PerCR/Cy5.5 (BioLegend, San Diego, CA, USA, #324214), anti-CD45-APC-Cy7 (BioLegend, #368515), and anti-CD29-Alexa Fluor ® 488 (BioLegend, # 303015) were verified according to the manufacturer's website.

Consensus Clustering Analysis of CRGs
Initially, 25 CRGs were identified from previous studies [24]. Then, we used the "ConsensusClusterPlus" package [39,40] to perform a consensus clustering analysis by the k-means algorithm to identify different CAF-associated subtypes. Furthermore, a proteinprotein interaction (PPI) analysis through the string website (https://cn.string-db.org/ (accessed on 1 January 2022) was constructed to determine the interplay of CRGs.

Correlations of Molecular Patterns with the Clinical Features and Prognosis of PC
We explored the correlation of molecular subtypes, clinical variables, and survival outcomes to assess the clinical value of the two CAF subtypes. The clinical features included age (≥65 and <65 years), gender (male and female), tumor location (left and right side), TNM stage (stage I-IV), KRAS mutation status (abnormal and normal), and TP53 mutation status (abnormal and normal). In addition, the differences in OS between the two subtypes were estimated by Kaplan-Meier analysis using the "survival" and "survminer" packages [41].

Association of Molecular Subtypes with Tumor Immune Microenvironment of PC
We assessed PC patients' immune and stromal scores using the ESTIMATE algorithm [42]. Then, the infiltrating fractions of 33 immune cell subtypes and four CAF subsets of each patient were computed with a single-sample gene set enrichment (ssGSEA) analysis algorithm [43].

Relationship of CAF Subtypes with Immunotherapy Responses in PC
As CAFs play a vital role in regulating tumor immune evasion [44], to further determine the association of CAF subtypes with immunotherapy responses of PC, we performed ssGSEA analysis to dissect the gene expression profiles of immunotherapy responses in the PC patients according to the nivolumab-responsive and anti-PD-1-resistant signatures [45].

Development of the CAF-Related Gene Risk Signature
We performed univariate and multivariate Cox regression analysis in TCGA-PAAD cohort to construct a novel CRG-based signature for predicting prognosis. Initially, the univariate Cox regression analysis was employed to assess the prognostic values of 25 CRGs in PC patients. Then, p < 0.05 was selected as a screening threshold, and 22 prognostic CRGs related to the survival of patients with pancreatic cancer were screened out. We further performed a multivariate Cox regression analysis and the penalized Cox regression model with the least absolute shrinkage and selection operator (LASSO) according to 22 prognostic CRGs. We obtained five hub prognostic CRGs (VCAN, COL1A2, ZNF469, SPARC, and FNDC1) and their corresponding coefficients. Finally, a scoring algorithm named CRG score was established to quantify the CAF state at the transcriptomic level. The CRG score was calculated as follows: CRG score = Σ (expression × correlative coefficient).
Additionally, we executed tROC curve analysis and independent prognostic analysis to validate the predictive capability of this novel CRGs-based signature and other clinical variables at 1-, 2-, and 3-year OS using the R package "survivalROC".

Expression Level Validation of CAF-Related Risk Gene Expression
Immunohistochemical results of CRGs involved in the risk signature were obtained from the HPA database (https://www.proteinatlas.org/ (accessed on 1 January2022) to validate CRG expression in normal and tumor tissue.

Identification of Immune Microenvironment Affected by CRGs
We investigated the association between CRG score and immune microenvironment with several common methods, including XCELL, TIMER, QUANTISEQ, MCPOUNTER, EPIC, CIBERSORT-ABS, and CIBERSORT.

Correlation of CRG Score Signature with Signal Pathways, Tumor Mutation, and Chemosensitivity
To identify the differences in somatic mutations of PC patients between high-and low-CRG-score groups, the mutation annotation format was created with the "maftools" R package [46]. We further examined the dependence of the CRG score, clinical outcome, and TMB. Moreover, to explore diversities in chemotherapy drug efficacy between the two subgroups, we estimated the half maximal inhibitory concentration values (IC 50 ) of chemotherapy drugs for each patient using the "pRRophetic" package [47], which is based on drug sensitivity data from the Genomics of Drug Sensitivity in Cancer dataset (https://www.cancerrxgene.org/ (accessed on 1 January2022)). A ridge regression model fitted the standardized expression data using predictor genes and the drug sensitivity (IC 50 ) values as the outcome variables [48].

Additional Bioinformatics and Statistical Analyses
We applied the Wilcoxon test to analyze the inter-group differences, conducted Spearman analysis for correlation tests, and performed the log-rank and Kaplan-Meier tests to draw survival curves. The R 3.6.3 software (Bell Laboratories, New York, NY, USA) and its corresponding packages were used to process, analyze, and present the data. By comparing different groups, p < 0.05 was considered to indicate statistical significance.

Conclusions
In this study, we systematically analyzed the genomic backgrounds and expression levels of CRGs and inferred their latent role in PC patients' prognosis and tumor microenvironment. We also constructed a novel CAF-associated gene signature as a robust biomarker to predict the prognosis, chemotherapeutic drug sensitivity, and immunotherapy impacts in PC. These results reveal the vital clinical significance of CRGs and put forward new ideas about the molecular classification of PC, which may be applied to precision medicine.