Anoikis-related CTNND1 is associated with immuno-suppressive tumor microenvironment and predicts unfavorable immunotherapeutic outcome in non-small cell lung cancer

Dataset acquisition

The normalized RNA-sequencing profile and clinical annotations of patients in TCGA-NSCLC cohort were downloaded from the UCSC Xena website (https://xenabrowser.net/datapages/). GSE42127 ²², the validation cohort, were downloaded from the Gene Expression Omnibus (GEO) portal (https://www.ncbi.nlm.nih.gov/geo/). Furthermore, the normalized gene expression profile of NSCLC clinical cohorts with anti-PD-1 therapy (GSE126044 ²³ and GSE135222 ²⁴) were also obtained from the GEO database. Also, immunotherapy cohorts of other epithelial tumors, including breast cancer (breast cancer: GSE173839 ²⁵), melanoma (PRJEB23709 ²⁶), and gastric carcinoma (PRJEB25780 ²⁷), were also obtained from public database. Samples with overall survival (OS) above zero-day were included in this research. For the immunotherapy cohort, diagnostic patients who received immunotherapy were selected for further analysis. In addition, the the anoikis-related genes were obtained from the Harmonizome portals ²⁸ (https://maayanlab.cloud/Harmonizome/, accessed on 12 October 2022). The information of datasets used in this study was listed in Table S1.

Single-cell RNA sequencing datasets analysis

The single-cell RNA sequencing datasets of 12 patients with NSCLC from the GSE150660 ²⁹, GSE127465 ³⁰ and GSE117570 ³¹ were downloaded. All additional analyses were performed using the Seurat (4.0.4, http://satijalab.org/seurat/) R toolkit ³², including quality control and all subsequent analyses. To eliminate the influence of abnormal cells and technical background noise on downstream analysis, cells were reserved if the expression of mitochondrial genes was greater than 10% or with detected genes less than 200 or greater than 5,000. Finally, a total of 47,359 cells were used for further analysis.

In order to minimize the technical batch effects among individuals and experiments, we used the “RunHarmony” function in R package harmony ³³ to integrate 47,359 cells from 12 NSCLC patients. The top 4,000 variable genes were used for principal component analysis (PCA) to reduce dimensionality. The dimensionality of the scaled integrated data matrix was further reduced to two-dimensional space based on the first 30 principal components (PCs) and visualized by t-Distributed Stochastic Neighbor Embedding (t-SNE). The cell clusters were identified based on a shared nearest neighbor (SNN) modularity optimization-based clustering algorithm with a resolution of 1, and all cells were divided into 26 clusters (Figure 1B). In order to recognize the types of these cells, some known markers, such as VWF for endothelial cells, EPCAM for epithelial cells, DCN for fibroblasts, CD3D for T cells, were used to verify the annotation of cell types (Figure 1C).

To validate the results found from the integrated scRNA-seq datasets of the GSE150660 ²⁹, GSE127465 ³⁰ and GSE117570 ³¹, another scRNA-seq datasets (GSE131907 ³⁴) including 11 NSCLC patients were downloaded from the GEO website. Cell filtration, integration, and annotation were followed the criteria used above.

Identification of differential expressed genes (DEG)

To identify the tumor cell-specific genes, “FindAllMarkers” function was performed. Genes with the |fold-change (FC)| ≥ 1.2, pct.1 ≥ 0.4, pct.2 ≤ 0.1, and adjusted P-values < 0.05 were identified as the tumor cell-specific genes.

In order to recognized the immunotherapeutic-related genes, the the R package “limma” ³⁵ was used to perform the differential expression analysis between the NSCLC patients who received (R) and not received remission (NR) after immunotherapy in the GSE126044 cohort. Genes with the | FC| ≥ 1.5 and adjusted P-values < 0.05 were identified as the immunotherapeutic-related genes.

In order to identify the DEGs for CTNND1-high and low groups respectively, the R package “limma” ³⁵ was used to perform the differential expression analysis. Genes with the FC ≥ 1.5 and adjusted P-values < 0.05 were defined as up-regulated genes for CTNND1-high group, while genes with the FC ≤ -1.5 and adjusted P-values < 0.05 were recognized as up-regulated for CTNND1-low group.

Assessment of immunological characteristics of the TME

The associations between CTNND1 and the immunological features of the TME was evaluated ³⁶. In order to assess the immunological characteristics of the TME, the ESTIMATE algorithm ³⁷, a method inferring tumor purity and stromal and immune cell from tumor samples based on bulk transcriptomic profile, was performed to assess tumor purity, ESTIMATE score, immune score, and stromal score. Besides, the information of immunomodulators including MHC signatures, receptors, chemokines, and immune-stimulators was collected from the previous studies ³⁸. To further deconstruct the immunological status of each patient, a set of signature genes of 29 immune cell types and immune-related pathways ³⁹ was used to estimate the infiltration levels of different immune cell populations and the activities of immune-related pathways and functions of each patient were calculated by utilizing the single-sample gene sets enrichment analysis (ssGSEA) in the R package “GSVA” ⁴⁰.

Cell-cell communication analysis

Cell-cell communications mediated by ligand-receptor complexes were critical to diverse biological processes, such as inflammation and tumorigenesis. To investigate the molecular interaction networks between different cell types, we used “CellPhoneDB” ⁴¹, a software to infer cell-cell communication from the combined expression of multi-subunit ligand-receptor complexes, to analyze the interactions between tumor cells and microenvironment cell subpopulations. The ligand-receptor pairs with a P value < 0.05 were remained for the assessment of relationship among different cell clusters.

Immunohistochemistry and semi-quantitative analysis

Lung cancer tumor microarray (TMA) HLugC120PT01, which was purchased from Outdo BioTech, contained 60 paired tumor and para-tumor samples. A total of 58 tumor samples were included in our research after removing the samples separated from the TMA, and the detailed clinic parameters of enrolled in-house patients were exhibited in Table S2. The use of the TMA was approved by the Clinical Research Ethics Committee in Outdo Biotech (Shanghai, China). The TMA was submitted for immunohistochemistry (IHC) assay to define the protein expression of CTNND1 in tumor and para-tumor tissues. The sections were then washed with xylene for three 5-min. The sections were rehydrated by successive washes in 100, 90 and 70% graded ethanol. Hydrogen peroxidase (0.3%) was used to block endogenous peroxidase activity for 20 min. The EDTA antigen repair solution was used for antigen repair. The primary antibody utilized in the study was anti-CTNND1 (1:200 dilution, Cat. sc-23873, Santa Cruz) and anti-CD8A (ready-to-use, Cat. PA577, Abcarta). Antibody staining was visualized with DAB and hematoxylin counterstain. Stained TMA was evaluated to define CTNND1 expression by two independent senior pathologists according to the immunoreactivity score standard ⁴². For the assessment of tumor-infiltrating CD8+ T cells, two senior pathologists estimated the CD8 score according to the criterion established by The Cancer Genome Atlas Network ⁴³. A CD8 score defined as the sum of the distribution and density scores (0-6) was calculated for each case. Samples with the CD8 score ≥ 3 (3, 4, 5, 6) are considered to be immune-hot, and samples with the CD8 score ≤ 2 (0, 2) are considered to be immune-cold.

Statistical analysis

All statistical analyses were handled using R software (version 4.0.4). The significant difference in continuous variables between the two groups was assessed using the Wilcoxon rank-sum test, while fisher exact test was used to measure the difference among categorical variables. For all analyses, a two-paired p-value < 0.05 was deemed to be statistically significant, and labeled with *p-value < 0.05, **p-value < 0.01, ***p-value < 0.001, and ****p-value < 0.0001.

Anoikis-related genes expressed specifically in tumor cells of NSCLC at single-cell atlas

Using the “Seurat” package to preprocess the single cell transcriptome data, a total of 47,359 cells were obtained (Figure 1A-1B). Markers of various cell types of NSCLC were collected, and the cells were annotated. The cells were divided into epithelial cells (EPCAM), T cells (CD3E), B cells (CD19), fibroblasts (COL1A1), endothelial cells (VWF), myeloid cells (CD14), and plasma cell (SLAMF7) (Figure 1C). Perform feature recognition on each annotated cluster of cell subpopulations to verify the accuracy of our cell type annotation (Figure 1D-1E). Furthermore, we used the ssGSEA algorithm to calculate the ANOIKIS score for the genes related to anoikis, and compared their differences between tumor cells and non-tumor cells identified previously. The results showed that the ANOIKIS score was significantly higher in tumor cells than in non-tumor cells (Figure 2A), thus confirming that anoikis is an important feature of tumor cells of NSCLC.

CTNND1 was identified as a tumor specific biomarker of immunotherapy for NSCLC

In order to explore the characteristic markers, and remodeling effect of tumor anoikis on the immune microenvironment of NSCLC, we analyzed the immune cohort data GSE126044. We obtained 3979 differential genes significantly related to the efficacy of immunotherapy through differential analysis. Next, we intersected immunotherapy-associated differential genes with anoikis-related genes, and the genes significantly expressed in tumor cells, obtaining a special biomarker: CTNND1, which can be inferred as an immunotherapy and tumor anoikis related biomarker (Figure 2B). The CTNND1 gene was validated in single cell atlas, and the results showed that its expression in tumor cells was significantly higher than that in non-tumor cells (Figure 2C-D). We binarized tumor cells based on whether the expression of CTNND1 was greater than 0, and compared the proportion of positive and negative cell rates in tumor and non-tumor cells. The results showed that the rate of CTNND1+ cells was significantly higher in the tumor cells than in non-tumor cells (Figure 2E). Finally, in the GSE126044 data, the expression of CTNND1 in the immunotherapy sensitive group was significantly higher than that in the ineffective group (Figure 2F).

Furthermore, another scRNA-seq dataset (GSE131907) including 11 NSCLC patients was used to validate the results found at the single-cell levels. Unsupervised clustering and cell annotation, we classified the cells into eight cell types (Figure 3A-3B, Supplementary Figure 1A-1B). Consistently, CTNND1 almost expressed on tumor cells (Figure 3C-3D). Meanwhile, after binarizing the cells into CTNND1+ and CTNND1- groups, the rate of CTNND1+ cells were significantly higher in the tumor cells than in non-tumor cells (Figure 3E), further supporting the viewpoint that CTNND1 was the biomarker of tumor cells.

CTNND1+ tumor cells communicated with microenvironment subpopulations frequently

Having observed that the expression of CTNND1 in the immunotherapy sensitive group was significantly lower than that in the ineffective group, we next sought to investigate the correlations between CTNND1 expression and the TME status of the NSCLC patients. Benefiting from the advantages of scRNA-seq technology, we could perform a high-resolution dissection of interactions among various subgroups of CTNND1+/- tumor cells and microenvironment subpopulation in the NSCLC patients based on the combining expression of multi-subunit ligand-receptor complexes. The number of interactions among different subpopulations was compared between CTNND1+ and CTNND1- tumor cells. Results showed that the CTNND1+ tumor cells presented significantly more interactions than CTNND1- tumor cells (Figure 4A-4B). Notably, we calculated the difference in interaction numbers among various CTNND1+/- tumor cells and microenvironment subpopulations. The CTNND1+ tumor cells showed more frequently crosstalk with other subpopulations, especially the immune subsets (Figure 4C). Results from the GSE131907 scRNA-seq dataset was in-kept with these findings (Supplementary Figure 2). These results collectively suggested that compared with CTNND1- tumor cells, CTNND1+ tumor cells had activated cell-cell communications, especially interactions with immune cells, which potentially take part in the formation of an immunosuppressive TME ⁴⁴.

We further identified the significant ligand-receptor interactions between CTNND1+ tumor cells and immune subsets using the “CellphoneDB” tool ⁴¹ (Figure 4D-4E). Results showed that CTNND1+ tumor cells communicated with myeloid cells via ANXA1-FPR1/3, SIRPA-CD47, and LGALS9-CD44, which have been reported to involve in the inhibition of anti-tumor response ^45-47. In addition, the CTNND+ tumor cells mediated the dysfunction of T cells via inhibitory ligand-receptor pairs, such as ANXA1-FPR3 ^{48, 49}. In addition, compared with CTNND1- tumor cells, CTNND1+ tumor cells showed the highest cell-cell communication strength with stromal cells (fibroblasts and endothelial cells) (Figure 4B-4C), and presented significantly more specific interactions with stromal cells (Supplementary Figure 2A-2B). For example, CTNND1+ tumor cells communicated with stromal cells via many ligand-receptor pairs, such as COL1A1_a1b1 complex, which involved in the activation and differentiation of endothelial and fibroblasts ^{50, 51}, suggesting that the accumulation of CTNND1+ tumor cells will be accompanied by the enrichment of stromal cells, leading to the formation of desert TME. Besides, our results also showed that some interactions associated with tumor stemness, chemoresistance, and angiogenesis, such as WNT7B-FZD4 and NOTHC1-JAG1 ^{52, 53}, were detected between CTNND1+ tumor cells and stromal cells (Supplementary Figure 3A-3B). Some studies proved that selectively blocking JAG/NOTCH can disrupt angiogenesis by unique mechanisms to inhibit tumor growth ^{54, 55}. Totally, compared with CTNND1- tumor cells, CTNND1+ tumor cells interacted with microenvironment subpopulations frequently, suggesting that the complex crosstalk between CTNND1+ tumor cells and other subpopulations will contribute to the formation of TME.

High expression of CTNND1 was related to immunosuppressive status of NSCLC

Next, we explored the effect of CTNND1 on the remodeling of immune microenvironment of NSCLC. As shown in Supplementary Figure 4A-4D, we found that CTNND1 expression showed no significant differences between patients with different clinic parameters, suggesting that clinical parameters did not influence the expression levels of CTNND1. Then, we calculated the tumor purity and immune infiltration score in each sample through ESTIMATE, and carried out correlation analysis with the expression of CTNND1. The results showed that the expression of CTNND1 was significantly positively correlated with the tumor purity, and significantly negatively with tumor purity (Figure 5A-B). The samples were divided into high and low expression groups based on the median expression of CTNND1, and the results showed that the immune infiltration score of the high expression patients was significantly lower than that of the low expression patients (Figure 5C-D). Furthermore, we investigated the relationship between CTNND1 and the degree of infiltration of various types of immune cells, which showed that the infiltration degree of immune cells was negatively correlated with CTNND1 expression (Figure 5E). Meanwhile, patients with high expression of CTNND1 showed resistance to immunotherapy, while 62.5% of patients with low expression showed sensitivity (Figure 5F). The results of the GSE135222 dataset also confirmed this conclusion (Figure 5G-H), and survival analysis showed that patients with high CTNND1 expression had poorer progression-free survival compared to those with low expression (Figure 5I). In addition, in other epithelial cell carcinomas, such as gastric carcinoma, breast cancer and melanoma, patients in the CTNND1-high group also showed poor immunotherapeutic response than those in the CTNND1-low group (Supplementary Figure 5A-5C).

The median expression value of CTNND1 in NSCLC cases in TCGA was used as the threshold to divide the samples into two groups: CTNND1 high and CTNND1 low, and then performed differential analysis. The volcano plot shows up-regulated and down-regulated genes. The results of functional enrichment analysis showed that the ANOIKIS-related pathway in the CTNND1 high groups were significantly upregulated (Figure 6A-B). The samples from TCGA also showed the feature of low immune infiltration in the CTNND1 high expression group (Figure 6C). The levels of immune receptor activation, antigen presentation, immune activation, and cytokine between high- and low-CTNND1 expression groups were then compared by ssGSEA scores, in order to further investigate the relationship between CTNND1 and immune function. The results showed that the immune function of the CTNND1 high group was significantly lower than patients of the low expression (Figure 6E). Comparing the activation levels of various types of immune cells between the two groups, we found that the abundance of all immune cells was lower in high CTNND1 group, especially CD8+T cells (Figure 6F, Supplementary Figure 6 and Figure 7). To further investigate the role of this marker in immune suppression, we investigated the relationship between CTNND1 and various immune checkpoints, and the results showed a negative correlation with most immune checkpoints, indicating that high expression patients were ineffective in immunotherapy (Figure 6G). Finally, we compared the enrichment scores of various steps in the cancer immune cycle between CTNND1 high and low groups, our finding suggested a significant negative correlation between CTNND1 expression and various processes of the immune cycle (Figure 4I). To validate our conclusion, we conducted further analysis on the validation set GSE42127 (Figure 7). Our results indicate that CTNND1 expression is significantly negatively correlated with immune cells abundance and immune-related signaling pathways. Additionally, we observed that the expression of immune checkpoints in the high expression group is significantly lower than that in the low expression group, which may indicate ineffective immunotherapy.

Last but not least, given of the importance of genomic variants in diagnosis and therapeutic guidance, we compared the mutations between the CTNND1-high and low groups (Table S3). Results showed that some genomic variants, such as TP53 ⁵⁶ and KMT2D ⁵⁷, associated with worse clinical outcomes and tumor invasion were enriched in CTNND1-high group, while CTNND1-low groups had higher mutant frequency of some gene associated with immune infiltration and favorable outcomes, such as CACNA1C ⁵⁸ and LMO7 ⁵⁹ (Supplementary Figure 8). These results were consistent with the phenotypes found at the transcriptional levels, and further indicated that for those patients in the CTNND1-high group, targeting some mutations, such as TP53 can benefit patients.

Validation of CTNND1 expression and immuno-correlation in NSCLC

We next validated the expression patterns of CTNND1 in the in-house cohort. We found that CTNND1 was highly expressed in tumor tissues compared with para-tumor tissues (Figure 8A-8B) and CTNND1 was specifically expressed in NSCLC (Figure 8C-8D). We also divided NSCLC as immuno-hot and immune cold tumors based on the CD8 score. Obviously, CTNND1 was decreased in immuno-hot tumors (Figure 8E-8F) and negatively correlated with the CD8 score (Figure 8G). Overall, the negative correlation between CTNND1 expression and T cell infiltration could be validated in the in-house cohort, which largely increases the credibility of results by public cohorts.