A novel signature predicts prognosis and immunotherapy in lung adenocarcinoma based on cancer-associated fibroblasts

We obtained scRNA-seq data from the Gene Expression Omnibus (GEO) database (accession number GSE149655↗), which comprised four samples: two primary lung adenocarcinoma (LUAD) samples and two normal tissue samples. We filtered out single cells expressing fewer than 250 genes or those with any gene expressed in fewer than three cells. We also evaluated the percentage of mitochondria and rRNA using the PercentageFeatureSet function in the Seurat R package (25, 26). This resulted in a total of 12,554 cells for further analysis.

We collected transcriptome data, copy number variants (CNV), single-nucleotide variants (SNV), and corresponding clinical data of LUAD from The Cancer Genome Atlas (TCGA) database. We excluded samples lacking survival data or outcome status and included 500 tumor samples and 59 normal samples in the analysis. We also utilized two external validation cohorts: GSE72094↗ cohort with 398 samples and GSE26939↗ cohort with 115 samples after removing samples without follow-up.

From the literature, we identified ten cancer-associated pathways (HIPPO, Cell Cycle, MYC, NRF1, NOTCH, PI3K, RAS, TP53, TGF-Beta, and WNT) and analyzed their gene expression profiles in our dataset (27).

The Seurat package was used to re-analyze the scRNA-seq data of LUAD (28), with the aim of systematically characterizing the CAF signature. Firstly, expressed genes were log-normalized after removing cells with below 250 or over 6000 expressed genes. Then, the FindIntegrationAnchors function was employed to remove batch effects for the four samples. Non-linear dimensional reduction was performed using the uniform manifold approximation and projection method, with a resolution of 0.2 and 30 principal components selected. Subsequently, single cells were clustered into different subgroups using the FindNeighbors and FindClusters (dim = 30 and resolution =0.2) functions. UMAP dimensional reduction was performed using the RunUMAP function. Fibroblasts were annotated based on four marker genes, including FAP, PDGFRB, ACTA2, and NOTCH3. The fibroblasts were re-clustered using the same algorithm of FindClusters and FindNeighbors functions. Marker genes for each CAF cluster were defined with the FindAllMarkers function by comparing different clusters with minpct = 0.35, logFC =0.5, and adjust p-value<0.05. We also used the CopyKAT R package (29) to analyze the CNV characteristics of CAFs clusters and distinguish between tumor cells and normal ones. Finally, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis on the marker genes using the clusterProfiler package (30).

Firstly, the limma package (31, 32) was used to identify differentially expressed genes (DEGs) between normal and tumor tissue based on |log2(FoldChange)|>1 and a false discovery rate (FDR)<0.05. Then, correlations between CAF clusters and DEGs were evaluated, followed by the identification of key CAF-related genes with p<0.01 and cor>0.4. To identify prognosis-related genes, univariate Cox regression analysis was conducted using the survival package (33). To reduce the number of genes, the least absolute shrinkage and selection operator (lasso) was performed (34–36). Multivariate Cox regression analysis was conducted using the stepwise regression method to establish a CAF-based risk signature, which was calculated using the formula: 0.123CLEC3B+0.114EXO1 + 0.103CCNB1+-0.177CD302. Patients were classified into low- and high-risk groups using zero-mean normalization. The predictive value of the risk signature was evaluated using receiver operating characteristic curve (ROC) analysis with the timeROC package (37, 38).

Following the univariate and multivariate Cox regression analysis based on the risk signature and clinicopathological features (39), we constructed a novel nomogram to predict the prognosis of LUAD using variables with p<0.05 in the multivariate Cox model. We evaluated the predictive accuracy of the model by generating a calibration curve.

We comprehensively assessed the correlation between the risk signature and the tumor immune microenvironment (TIME) using several algorithms, including CIBERSORT, EPIC, MCPCOUNTER, and TIMER. The “estimate” R package was used to calculate stromal scores, immune scores, and estimate scores (stromal scores + immune scores) to evaluate differences in the tumor microenvironment of patients. Additionally, we estimated the proportions of 22 immune cell subtypes using the CIBERSORT algorithm based on the TCGA cohort.

Anti-PD-1 or anti-PD-L1 checkpoint inhibition therapy has gained increasing attention as a crucial component of immunotherapy. To evaluate the performance of the risk signature in predicting responsiveness to immunotherapy (immune checkpoint blocks), we collected transcriptomic data as well as corresponding clinical data from patients who received anti-PD-L1 therapy from the IMvigor210 cohort. We also downloaded transcriptomic data from the GSE78220↗ cohort, which included melanoma patients who received anti-PD-1 checkpoint inhibition therapy before treatment.

All patients conferred their informed consent before being enrolled in the study. Sample collections were conducted following procedures approved by the Ethics Committee of Jiangsu Province People’s Hospital (2019-SR-156). Lung adenocarcinoma cell lines including A549 and H1299 cells was purchased from ATCC. All cells were cultured using Ham’s F-12K and RPMI 1640 medium (Gibco, USA), supplemented with 10% FBS (HyClone Sera, USA) and 1% penicillin‐streptomycin (Sangon Biotech, China), and maintained in an atmosphere containing 5% CO2 at 37°C. The EXO1 siRNA expression vector and scrambled siRNA nontarget control were obtained from Genewiz (China). Plasmids were transfected using Lipofectamine 3000 (Thermo Scientific, USA), as per the manufacturer’s protocols (40).

Total RNA was extracted from the cell lines using TRIzol in accordance with the manufacturer’s instructions (15596018, Thermo). Subsequently, cDNA was synthesized utilizing the PrimeScript TMRT kit (R232-01, Vazyme). Real-time polymerase chain reaction (RT-PCR) was performed using SYBR Green Master Mix (Q111-02, Vazyme). The expression levels of each mRNA were standardized to the level of mRNA GAPDH, and the quantification of expression levels was executed using the 2^–ΔΔCT method (41).

The suspension of cells was seeded in 96-well plates at a density of 5×10³ cells per well. After adding 10 μl of CCK-8 labeling agent (A311-01, Vazyme) to each well, the plate was incubated for 2 hours in the dark at 37°C. Cell viability was evaluated by measuring absorbance at 450 nm at 0, 24, 48, 72, and 96 hours using an enzyme-labeled meter (A33978, Thermo). The experiment was performed using a 96-well plate with 2×10⁴ treated cells in each well, after the cells had adhered to the wall. The 5-Ethynyl-2’-deoxyuridine (EdU) assay was performed according to the manufacturer’s instructions (Ribobio, China), and cell proliferation was quantified using an inverted microscope.

The transfected cells were seeded in 6-well plates and cultured in a cell incubator until they reached 95% confluence. Each well was gently scraped using a sterile 200 μl plastic pipette tip, and any unattached cells and debris were rinsed twice with PBS. The breadth of the scratch wounds was measured using the Image J program, and photographs were taken at 0 h and 48 h. For the cell invasion and migration experiments, treated A549 and H1299 cells (2×10⁴) were incubated in the upper chamber of 24-well plates for 48 hours. The top surface of the plate was either precoated with matrigel solution (BD Biosciences, USA) or left untreated to assess the cells’ ability to invade and migrate. After removing the cells from the top surface, the remaining cells on the bottom layer were fixed with 4% paraformaldehyde and stained with 0.1% crystal violet (Solarbio, China).

All statistical analyses were performed using R software (version 4.1.0). Wilcoxon test was used for comparing two groups, while Spearman or Pearson correlation was used for correlation matrices. The Log-rank test was used to assess survival differences through K-M curves, where statistical significance was defined as p-value < 0.05.

The flow chart of our study is depicted in Figure 1. Following initial screening, a total of 12,554 cells were obtained from scRNA-seq data. Detailed results of data preprocessing are presented in Figure S1. Firstly, after performing log-normalization and dimensionality reduction, we identified 31 subpopulations (Figure 2A). Using four marker genes (FAP, FDGFRB, ACTA2, and NOTCH3), we further identified five CAF populations, as shown in Figure 2B. Cells collected from the four CAF populations were then separated for clustering and dimensionality reduction in subsequent research. Utilizing the same clustering algorithm, we discovered five clusters, as depicted in Figure 2B. Moreover, after performing the R package ‘FindVariableFeatures’, we obtained 756 DEGs from the five CAF clusters. The top 5 DEGs, which were characterized as CAF cluster marker genes, are exhibited in Figures 2C, D. Histograms illustrating the proportion of the five clusters in each cohort are shown in Figure 2E. Furthermore, KEGG analysis revealed that these DEGs were enriched in divergent pathways such as ECM-receptor interaction, focal adhesion, proteoglycans in cancer, PI3K-Akt signaling pathway, and others, as presented in Figure 2F. Finally, the distribution of tumor and normal cells based on the five CAF clusters according to the CNV characteristics is shown in Figure 2G.

We explored the characteristics of ten tumor-related pathways in the five CAF clusters to elucidate the correlations between tumor progression and the CAF clusters. GSVA was employed to investigate the underlying mechanisms involved in the progression and prognostication of LUAD. GSVA scores of those pathways were calculated based on different CAF clusters, and the results are presented in Figure 3A. As shown in Figure 3B, the percentage of normal cells in CAF_0 cluster was the highest, while the ratio of malignant cells from CAF_1 was significantly higher than that in the others. Significant differences were only identified among the CAF_0, CAF_2, and CAF_4. Furthermore, GSVA scores were analyzed based on the ten tumor-related pathways between normal and malignant cells in each CAF cluster, with slight differences observed in CAF_2 and CAF_4 (Figures 3C–F). (The results of GSVA score analysis in CAF_0 is shown in Figure S2).

Moreover, to explore the correlations between the CAF clusters and crucial clinicopathologic features, we analyzed the ssGSEA score of the marker genes (the top 5 DEGs referred to in Figures 2C, D) of each CAF cluster according to the TCGA cohort. The results showed that tumor samples had a significantly higher score compared with normal ones in each cluster (Figure S2B). Using the survminer R package, LUAD samples of TCGA cohort were divided into low-and-high CAF score groups based on the optimal cut-off value. The samples in the low-CAF score group had a significantly worse prognosis in CAF_0, CAF_1, and CAF_3. Although no significant differences were observed in the other two clusters, there was a trend that patients in high-CAF score groups shared a better prognosis (Figure S2C). Furthermore, other clinicopathologic features, including T. stage, N. stage, M. stage, and stage, were included in the analysis. Only slight differences were observed between low-and-high CAF group patients (Figure S3). However, patients in high-CAF groups tended to share favorable clinicopathologic features.

To construct a risk signature, we first screened for DEGs between normal and tumor tissues. A total of 1731 DEGs were obtained, with 725 up-regulated and 1006 down-regulated (Figure 4A). Out of these, 492 genes were significantly correlated with the prognosis-related CAF clusters. Using univariate Cox regression analysis, we further evaluated the prognosis of each gene, identifying 49 genes related to risk factors and 62 genes exhibiting protective values (Figure 4B). To reduce the number of genes, we conducted Lasso Cox regression analysis, resulting in 4 genes with lambda=0.074 (Figure 4C). Finally, we used the stepwise regression method to construct the risk signature after multivariate Cox regression analysis. The signature consists of 4 genes (Figures 4D, E), namely Exonuclease 1 (EXO1), Cyclin B1 (CCNB1), C-Type Lectin Domain Family 3 Member B (CLEC3B), and Type I C-type lectin receptor CD302. The final signature formula is as follows: RiskScore = -0.123CLEC3B+0.114EXO1 + 0.103CCNB1+-0.177CD302. Using z-mean normalization, we calculated the risk score for each sample, dividing patients into high and low-risk groups. Kaplan-Meier survival analysis showed that low-risk patients had significantly better survival outcomes compared to high-risk patients, not only in the TCGA cohort (Figure 4F) but also in the GSE72094↗ (Figure 4G) and GSE26939↗ (Figure 4H) cohorts. Additionally, based on the TCGA and GEO cohorts, the AUC values of the signature for 1-3-5-year survival were satisfying, indicating the model’s excellent predictive power (Figures 4F-H). We also presented the distribution of risk score, patient survival status, and expression of hub genes in the TCGA cohort in Figure S3A. Similarly, the results of GSE72094↗ and GSE26939↗ were shown in Figures S3B-C.

To construct a more accurate predictive model, we integrated the risk score with clinicopathological characteristics using both univariate and multivariate Cox regression analyses. The multivariate analysis demonstrated that the risk signature was the most significant independent prognostic factor for lung adenocarcinoma (p-value < 0.001), followed by N-stage (Figures 5A, B). Accordingly, we developed a novel nomogram incorporating T-stage, N-stage, and risk score (Figure 5C), which demonstrated strong predictive power for actual survival outcomes according to calibration plot analysis (Figure 5D). TimeROC analysis in the TCGA cohort confirmed that the AUC of the nomogram and risk score exceeded that of other indicators (Figure 5E).

After examining the clinicopathologic features (Age, Gender, Stage, T-stage, N-stage, and M-stage) between high- and low-risk groups, we found that gender, T-stage, N-stage, and stage were significantly associated with the risk signature (). These findings were consistent with previous studies that identified gender as a risk factor for LUAD, with males being more likely to be in the high-risk group (). Moreover, patients in the high-risk group tended to have more advanced clinical stages (). 4 4 4

Exploring SNV mutations in lung adenocarcinoma based on the TCGA cohort to investigate the SNV mutations in lung adenocarcinoma, the top 20 genes with the highest mutation frequency were analyzed based on the TCGA cohort, as shown in Figure 6A. Subsequently, the SNV mutations of the four genes in the risk signature were examined. As displayed in Figure 6B, EXO1 had the highest mutation frequency, with Missense-Mutation being the most common type of mutation. Conversely, no mutations were observed in CD302. Additionally, the probability of co-occurrence of the 10 most mutated genes and the risk genes (except for CD302) was assessed, and the results indicated a low likelihood of co-occurrence of mutations in these three genes. However, EXO1 was found to significantly co-occur with MUC16, CSMD3, RYR2, ZFHX4, and USH2A (Figure 6C). Further analysis revealed that only a few samples had loss/gain of CNV based on the four genes (Figure 6E). The fraction of the pathway affected by these risk genes was also calculated in the TCGA cohort (Figure 6D). Moreover, the relationships between the risk genes and several molecular signatures of LUAD were explored to demonstrate the links between the risk genes and LUAD. The results indicated that EXO1 and CCNB1 were positively correlated with molecular signatures such as Aneuploidy Score, Homologous Recombination Defects, and Fraction Altered, while CLEC3B and CD302 were negatively correlated with these signatures (Figure 6F).

Based on these four genes from the risk signature, Gene Set Enrichment Analysis was performed. The results showed that 16 pathways were significantly correlated with these four genes in total (Figures 7A, B), such as the p53 signaling pathway, cell cycle, and DNA replication. Similar to the results obtained previously, EXO1 and CCNB1 were positively correlated with these pathways, while CD302 and CLEC3B were negatively related to them. The GSEA score was estimated based on the high-and-low-risk subgroups (Figures 7C, D). Centromere Complex Assembly, Cell Cycle Checkpoint Signaling, and Cell Cycle G2 Phase Transition were significantly enriched in the high-risk group. In contrast, Positive Regulation of Lipase Activity, Axoneme Assembly, and Cilium Movement were significantly enriched in the low-risk group. Finally, the results of KEGG and GO analysis are shown in Figure 7E.

After investigating the landscape of immune and stromal cell infiltrations in both low- and high-risk groups, we found that Figure 8A illustrates how patients in the low-risk group have higher proportions of immune and stromal cell infiltrations compared to those in the high-risk group. Moreover, using the CIBERSORT algorithm (42, 43), we calculated the immune cell proportions between the high- and low-risk groups (Figure 8B) and found that patients in the high-risk group significantly shared higher proportions of CD8 T cells, activated memory CD4 T cells, activated NK cells, Macrophages (M0), and Macrophages (M1). On the other hand, B cells, resting memory CD4 T cells, Monocytes, resting dendritic cells, and Activated mast cells were significantly enriched in the low-risk group.

We then explored the relationship between risk genes and immunity. Our results showed that EXO1 and CCNB1 had a significantly negative relationship with the majority of T cells, while CLEC3B and CD302 were remarkably correlated with Macrophages (Figure 8C). Additionally, correlation analysis presented that EXO1 and CCNB1 were negatively linked with Stromal Score, Immune Score, and ESTIMATE Score. In contrast, the other risk genes exhibited the opposite trend (Figures 8D, E, G). Finally, Figure 8F revealed the correlation between the four risk genes and the 75 immune-related genes.

We analyzed the response to PD-L1 blockade immunotherapy in the IMvigor210 and GSE78220↗ cohorts after assessing immune infiltrations. The 348 patients from the IMvigor210 cohort presented different responses to anti-PD-L1 receptor blockers, including stable disease (SD), partial response (PR), complete response (CR), and progressive disease (PD). We found that CR/PR patients had lower risk scores than SD/PD patients (Figure 9B). Additionally, in the low-risk group, the proportion of SD/PD patients was lower than that in the high-risk group (Figure 9C). Our analysis of the IMvigor210 cohort revealed that patients in the low-risk group had significantly better clinical outcomes than those in the high-risk group (Figure 9A). Furthermore, we identified significant survival differences between different risk groups not only in Stage I+II but also in Stage III+IV patients (Figures 9D, E).

To confirm our findings, we enrolled the GSE78220↗ cohort into further analysis. Corresponding with the results obtained in IMvigor210, PR/CR patients had lower risk scores and shared a lower percentage in the high-risk group (Figures 9F–H).

In order to further elucidate the function of EXO1 in LUAD, we conducted in vitro research to scrutinize the function of EXO1 in LUAD cells. We gauged the degree of EXO1 expression after 24 hours of transfection via qRT-PCR to assess the efficacy of siRNA-mediated EXO1 knockdown in A549 and H1299 cell lines. As compared to the NC group, we observed a marked reduction in the expression of EXO1 in A549 and H1299 cells upon treatment with siRNA (Si-1 and Si-2) sequences (P < 0.001) (Figures 10). Correspondingly, the CCK-8 assay revealed that suppression of EXO1 significantly curbed the viability of A549 and H1299 cells as compared to control cells (Figure 10A). The findings of the EdU staining assay provided further evidence that inhibition of EXO1 expression impeded the proliferation of A549 and H1299 cells relative to the NC group (Figure 10B). This implies that EXO1 might play an indispensable role in the proliferation of LUAD cells. The transwell experiments confirmed that EXO1 knockdown considerably reduced the migration and invasion of A549 and H1299 cells (Figure 10C). The scratch-wound healing experiment also produced congruent results, wherein decreased EXO1 expression led to a noteworthy deceleration in the rate of wound healing in cells (Figure 10D). To ensure the accuracy and consistency of the results, all tests were performed in two LUAD cell lines (A549 and H1299), and all data were presented as means with standard deviations from three independent experiments. *P < 0.05, **P < 0.01, ***P < 0.001.

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 within the article/. 1