Fibroblast-derived LPP as a biomarker for treatment response and therapeutic target in gastric cancer.

Association of tumor microenvironment and immune checkpoint (e.g., PD-L1) is important for immune escape, impacting chemotherapy and immunotherapy efficacy. We aimed to investigate biomarkers and therapeutic targets against treatment resistance in gastric cancer. Abundances of tumor-infiltrating immune cells were estimated in multiple datasets. Three patient subgroups (A, B, and C) were identified based on seven types of PD-L1- and IFN-γ-associated immune cells. Patients yielded increased prognosis from subgroup A to C (p = 0.027). Subgroup A was characterized by high activated CD4+ memory T cell infiltration, while more resting CD4+ memory T cells were in subgroup C. Further, a risk score was developed for prognostication. Lipoma preferred partner (LPP), as the hub gene in subgroup-related regulatory network, was upregulated (p < 0.01) and was associated with high risk score (p < 0.001) and poor survival (p < 0.05). Bioinformatics analyses and experiments found that LPP expressed restrictively in fibroblasts and associated with activated CD4+ memory T cell infiltration and tumor growth. High-LPP patients yielded fewer benefits from chemotherapy or immunotherapy, compared with the low-LPP group. We finally identified 28 compounds as sensitive drugs for high-LPP patients. Our findings suggested LPP might be a biomarker for treatment response and therapeutic target in gastric cancer.

Introduction

Recently, tumor immune microenvironment has been validated to play increasingly critical roles in the pathogenesis of multiple malignances, making remarkable influences on patient managements and therapy strategies in the clinical setting.1, 2, 3 Application of chemotherapy, with combining immune checkpoint blockade such as anti-programmed cell death protein 1 (PD-1) and anti-programmed cell death ligand 1 (PD-L1), has been demonstrated with groundbreaking efficacies in multiple tumor entities,, including gastric cancer (GC)., However, the response rates and the benefits vary among GC patients receiving adjuvant chemotherapy and immunotherapy. Thus, it is urgently needed to investigate the principle of resistance to adjuvant treatment and to find effective therapeutic targets to improve the benefits from chemotherapy and immunotherapy.

Several studies have reported that the heterogeneity of the immune microenvironment within and among individual tumors is greatly responsible for such variations of treatment responses and patient outcomes., For example, the tumor environment context observed at initially diagnosis could reflect the immune responses and chemotherapy benefits during subsequent clinical managements, and changes in the abundances of tumor-infiltrating immune cells, such as CD4+ T cells, CD8+ T cells, and macrophages, were significantly associated with clinical outcomes in various tumor types, including GC., These evidences imply that characterizing the tumor immune microenvironment and unraveling the intrinsic heterogeneity for abundances of tumor-infiltrating immune cells could make an important opportunity to advance understanding the underlying link between the efficacy of adjuvant treatment and the principle of immune surveillance and escape.

The major objective in this study was to identify prognostic biomarkers in GC for predicting response to chemotherapy and immunotherapy through uncovering the association between the heterogeneity of tumor immune microenvironment and the complexity of genomic alterations, as well as their impacts on the regulatory network underlying the pathology of GC. We assumed that integrating multi-omics profiles with tumor microenvironment characteristics could not only help to highlight the molecular principles regarding tumor cells escape immune surveillance in the presence of multiple types of tumor-infiltrating immune cells, but also identify the underlying relevant biomarker for treatment response. Our main deliverable of the biomarker lipoma preferred partner (LPP) might serve as a prognostic predictor for individual selection for chemotherapy and immunotherapy and could provide an underlying therapeutic target regarding adjuvant treatment resistance in GC.

Results

PD-L1- and IFN-γ-associated subgroups in GC

Abundances of 22 tumor-infiltrating immune cell types were estimated for each tumor sample from a total of 892 GC patients in five cohorts, including TCGA-STAD (n = 350), GSE26899 (n = 93), GSE26901 (n = 109), GSE66229 (n = 300), and GSE28541 (n = 40). We found that seven types of immune cells (i.e., CD8+ T cells, activated CD4+ memory T cells, follicular T helper cells, gamma delta T cells, activated NK cells, M1 macrophages, and neutrophils) had a significantly positive correlation, and four types of immune cells (i.e., plasma cells, resting CD4+ memory T cells, regulatory T cells, and M2 macrophages) were negative correlated with PD-L1 expression level, in all of five cohorts (Figure 1A). Similar analyses were also performed for IFN-γ, a well-known inducer of PD-L1 transcription that is secreted by activated T and NK cells, and the results also revealed consistent correlation patterns with that of PD-L1 in the above five cohorts (Figure 1B).

Figure 1

PD-L1 and IFN-γ associated patient subgroups in GC

(A) Dot plot showing Spearman's correlation between PD-L1 expression with abundances of 22 types of tumor-infiltrating immune cells in five independent GC cohorts. (B) Dot plot showing Spearman's correlation between PD-L1 expression with abundances of 22 types of tumor-infiltrating immune cells. (C) Unsupervised hierarchical clustering analysis revealing three clusters based on seven types of LASSO-selected, PD-L1- and IFN-γ-associated tumor-infiltrating immune cells in TCGA-STAD cohort. (D and E) Box plots showing distinct PD-L1 and IFN-γ transcript levels across subgroups A, B, and C in TCGA-STAD cohort. The boxes represent the median ± 1 quartile, with whiskers extending to the most extreme point within 1.5 interquartile range from the box boundaries. ∗∗p <0.01; ∗∗∗∗p <0.0001. (F) Kaplan-Meier plot showing distinct overall survival outcomes for GC patients in three subgroups.

We then performed least absolute shrinkage and selection operator (LASSO) regressions to prioritize the relevant tumor-infiltrating immune cell types with PD-L1 and IFN-γ expression level in these cohorts. And seven types of the most relevant immune cells, including resting CD4+ memory T cells, activated CD4+ memory T cells, gamma delta T cells, monocytes, M1 macrophages, resting mast cells, and neutrophils, were identified. The tumor-infiltrating abundances of seven immune cell types were used for unsupervised hierarchical clustering analysis. The results showed that patients were clustered into three distinct subgroups, defined as subgroup A, B, and C respectively (Figure 1C). In particular, samples in subgroup A yielded a higher proportion of activated CD4+ memory T cells, while more resting CD4+ memory T cells were found in subgroup C. Interestingly, both levels of PD-L1 and IFN-γ expression were significantly distinct among three subgroups, gradually decreasing from subgroup A to C (Figures 1D and 1E). Survival comparison analysis showed patients in subgroup A yielded the worst prognosis, whereas those in subgroup C yielded the best survival outcomes (Figure 1F).

Genomic characteristics of PD-L1- and IFN-γ-associated subgroups in GC

Landscapes of mutation and copy number variation, and expression profiles were analyzed across three subgroups using the TCGA-STAD cohort. The level of tumor mutational burden was significantly different among distinct subgroups, with a relatively high and low level in subgroup A and C, respectively (Figure 2A). We further examined ten genes with highest mutation frequency in each subgroup (Figure 2B). Five genes (TTN, TP53, MUC16, CSMD3, and LRP1B) simultaneously shared high mutation frequency in three subgroups. The oncoplots display mutational profiles for those genes with frequent mutation in each subgroup (Figures 2C–2E). In all three subgroups, TTN and TP53 were commonly observed with frequent missense mutation. In specific, subgroup A had more ARID1A frame shift mutation and PIK3CA missense mutation than the other two subgroups. Unlike subgroup A, MUC16 missense mutation tended to be more widespread in subgroups B and C. Interestingly, we found more mutational co-occurrence in subgroup A, such as PIK3CA and ARID1A mutations (Figure 2F), whereas more mutually exclusive mutation events were observed in subgroup B (e.g., ARID1A and TP53) (Figure 2G); and in subgroup C, most of the mutational events were relatively independent without interactions (Figure 2H).

Figure 2

Mutation landscape among subgroups A, B, and C

(A) Box plot showing distinct levels of tumor mutational burden across subgroups A, B, and C. The boxes represent the median ± 1 quartile, with whiskers extending to the most extreme point within 1.5 interquartile range from the box boundaries. ∗p <0.05; ∗∗p <0.01. (B) Venn diagram showing the distribution of top 10 genes with highest mutation frequency in each subgroup. (C–E) Oncoplot displaying mutation profile of a union set of the top 10 genes with highest mutation frequency in subgroups A, B, and C, respectively. (F–H) Graphs showing mutually exclusive and co-occurrence of mutational events in subgroups A, B, and C, respectively.

Copy number variation and expression profiles were compared between subgroup A and C, in which PD-L1 and IFN-γ expression levels, tumor mutational burdens, and survival outcomes were extremely distinct (Figure S1A and Table S1). Significant amplification for genomic region 17q21.31 (BRCA1) and 4p16.1 (WDR1) and deletion of region 3p14.2 (FHIT) were found in subgroup A. In subgroup C, amplification of region 1p31.1 (NEGR1) and 3q29 (ATP13A4) and deletion of region 17p11.2 (FAM106A, USP32P2, and CCDC144B) were observed. After comparison of expression profile between the two subgroups, a total of 339 and 375 genes were detected with significant upregulation in subgroup A and C, respectively (Figure S1B and Table S2). Functional enrichment analyses showed genes upregulated in subgroup A were enriched in pathways related to immunological activities such as antigen processing and presentation, and immune cell differentiation, while genes upregulated in subgroup C were enriched in pathways associated with digestive metabolisms including mineral absorption and gastric acid secretion (Figure S1C and Table S3).

PD-L1- and IFN-γ-associated risk score for prognostication in GC

Among the aforementioned differentially expressed genes between subgroups A and C, 45 genes were identified as prognostic genes by using Cox regression analyses, including 14 favorable genes and 31 risk genes (Table S4). Then, among them, nine favorable and five risk genes were selected as optimal signature genes by using LASSO-Cox algorithm (Figure 3A). The final risk score was developed based on the expression level of signature genes and corresponding coefficients (see materials and methods for details). The optimal cutoff for risk score was identified as cutoffrisk-score = 0.5358 for patient stratification, whereby patients were classified into low- and high-risk groups with distinct survival outcomes (Figures 3B and 3C). Survival comparison showed the low-risk group yielded significantly better survival outcomes than the high-risk group (log rank p < 0.0001). Stratification efficacy of the risk score was also validated by using another independent dataset GSE26899, and similar results were found (Figures 3D and 3E).

Figure 3

Risk score for prognosis stratification in GC

(A) Graph presenting hazard ratios and corresponding coefficients for 14 LASSO-selected signature genes using TCGA-STAD cohort. (B) Graph displaying the distribution of risk score and survival status (top) and expression profile of signature genes (bottom) for GC patients in TCGA-STAD cohort. (C) Kaplan-Meier plot showing distinct overall survival outcomes between the high- and low-risk groups in TCGA-STAD cohort. (D) Graph showing the distribution for risk score and survival status (top) and expression profile of signature genes (bottom) for GC patients in validation cohort GSE26899. (E) Kaplan-Meier plot showing distinct overall survival outcomes between the high- and low-risk groups in GSE26899 cohort.

LPP as a hub node in PD-L1- and IFN-γ-associated lncRNA-miRNA-gene regulatory network

We further constructed a PD-L1- and IFN-γ-associated lncRNA-miRNA-gene regulatory network to identify the critical component from the complex mechanism network that established the underlying foundation for the distinct genetic characteristics and outcomes among the subgroups (see materials and methods for details). A total of 59 miRNAs and 19 lncRNAs were found with dysregulation between subgroups A and C (Figure 4A). Figure 4B displays the predicted regulatory pairs between dysregulated miRNAs and genes (miRNA-gene) from three databases (i.e., miRTarBase, miRDB, and TargetScan). After merging the lncRNA-miRNA interactions from the StarBase database, the final lncRNA-miRNA-gene regulatory network was constructed (Figure 4C). Based on the network topology, we identified the gene LPP as the hub node in this network. The LPP-cored subnetwork is illustrated in Figure 4D. Interestingly, we found the LPP expression gradually decreased from subgroup A to subgroup C (Figure 4E), which was similar with PD-L1 and IFN-γ expression. Notably, patients in the high-risk group harbored higher LPP expression than those in the low-risk group (Figure 4F). These findings indicated that LPP might play crucial roles in the oncogenesis and development of GC.

Figure 4

LPP as a hub node in lncRNA-miRNA-gene regulatory network

(A) Volcano plot displaying miRNAs and lnRNAs dysregulated in subgroups A and C. (B) Venn diagram summarizing the prediction pairs from miRTarBase, miRDB, and TargetScan for upregulated miRNAs in subgroup A and C. (C) Visualization of the whole lncRNA-miRNA-gene regulatory network by using the Cytoscape software. Nodes in yellow, blue, and red indicate lncRNAs, miRNAs, and genes, respectively. Edges between the nodes indicate significant regulatory interactions for related molecules. (D) Graph displaying the LPP-cored subnetwork derived from the whole regulatory network based on the hub node identification. (E) Box plot showing distinct LPP levels among subgroups A, B, and C in TCGA-STAD cohort. (F) Box plot showing distinct LPP expression between high- and low-risk groups in TCGA-STAD cohort. The boxes represent the median ± 1 quartile, with whiskers extending to the most extreme point within 1.5 interquartile range from the box boundaries. ∗∗p <0.01; ∗∗∗∗p <0.0001.

High expression of LPP is an independent risk factor for poor prognosis in GC

Transcript level of LPP was compared between tumorous and normal tissues in 31 cancer types from the TCGA project (Figure 5A). The results showed LPP was significantly upregulated in four cancer types and downregulated in nine malignances. In particular, LPP was observed with significant overexpression in GC. Consistent results were also found in GC patients from independent datasets GSE66229 and GSE29272 (Figures 5B and 5C). By employing a previously developed method, maximally selected rank statistics, which allowed to optimally classify patients with the most distinct prognoses, we divided patients into two groups according to their LPP expression levels (Figure 5D). Survival comparison found the high-LPP group yielded worse survival than the low-LPP group (log rank p = 0.026) (Figure 5E). Similar results were also found in another two datasets GSE62254 and GSE84437 (Figures 5F–5G). Multivariate Cox regression analyses showed high expression of LPP was an independent risk factor for poor prognosis of GC patients in the above three datasets (Figure S2).

Figure 5

Overexpression of LPP is associated with poor prognosis in GC

(A) Graph presenting LPP expression differences between cancerous and normal tissues for 31 TCGA cancer types. Colors of the abbreviation of TCGA cancer type on the top indicate LPP dysregulation in tumors (red, upregulated; green, downregulated). Expression comparison in GC is marked by a blue box. (B and C) Box plot for comparison of LPP expression in normal and tumor tissues using datasets GSE66229 and GSE29272. The boxes represent the median ± 1 quartile, with whiskers extending to the most extreme point within 1.5 interquartile range from the box boundaries. (D) dot plot showing standardized log rank statistics for all candidate cutoffs of LPP expression in TCGA-STAD dataset. Dots in dark and red indicate patients classified in low- and high-LPP groups respectively, under the classification by the optimal cutoff. (E) Kaplan-Meier plot showing overall survival differences between low- and high-LPP groups in TCGA-STAD cohort. (F and G) Kaplan-Meier plot showing distinct overall survival outcomes between low- and high-LPP groups in GSE62254 and GSE84437 cohorts.

Fibroblast-derived LPP associates with tumor progression and resistance to chemotherapy and immunotherapy in GC

We used the immunohistochemistry from The Human Protein Atlas portal to observe the protein expression of LPP in GC (Figure S3A). The staining showed the LPP protein was strongly and restrictively expressed in stromal fibroblasts in tumor entity, while it was nearly not detected in normal stomach tissue. We also found high LPP expression was associated with increased fibroblast infiltration in the tumor entity (Figure S3B). Consistently, LPP expression was positively correlated with six well-known markers of fibroblasts that are widely used for fibroblast identification and quantification with sequencing data (Figure S3C). As mentioned earlier, LPP expression increased in subgroup A, which was characterized by poor prognosis with high infiltration of activated CD4+ memory T cells, and decreased in subgroup C, in which a better prognosis with high infiltration of resting CD4+ memory T cells was found. Herein, we consistently found that high expression of LPP was associated with more infiltration of activated CD4+ memory T cells and less infiltration of resting CD4+ memory T cells (Figures S3D and S3E).

To further validate the expression of fibroblast-derived LPP and its roles in the pathology of GC, we collected 15 GC samples coupled with corresponding normal tissues and used immunohistochemistry staining to analyze the protein expression of LPP. We observed similar results that the LPP protein was highly expressed in tumor tissues and was restrictively enriched in stromal fibroblasts (Figures 6A and 6B). The fibroblast-derived LPP was also confirmed by using a single-cell RNA-seq dataset of 26 primary GC samples containing 119,862 single cells (see materials and methods for details). The findings also revealed that LPP was predominantly enriched in fibroblasts, although low expression was observed in stem cells and endothelial cells (Figures 6C and 6D). Staining on serial tissue slides showed LPP protein expressed restrictively in fibroblasts, and it showed distinctly different distribution against those markers of endothelial cells, such as CD34 (Figure 6E). To further evaluate the role of fibroblast-derived LPP in GC development, the stable LPP knockdown fibroblast cell line was established for in vivo tumor growth experiments (Figure 6F). The results showed that the mice injected with GC cells mixed with LPP knockdown fibroblasts harbored smaller tumor volume than those injected with cancer cells mixed with the control (Figures 6G and 6H). These findings indicated fibroblast-derived LPP was associated with the tumor growth of GC.

Figure 6

Fibroblast-derived LPP involves in GC pathology

(A) H&E staining and immunohistochemistry for analyzing LPP expression in GC versus stomach tissue. (B) Box plot showing distinct levels of LPP protein expression in GC and normal tissues. The boxes represent the median ± 1 quartile, with whiskers extending to the most extreme point within 1.5 interquartile range from the box boundaries. (C) Graph showing cell clusters identified by using single-cell RNA-seq dataset GSE183904. (D) Violin plot presenting LPP expression level in each type of cell clusters. (E) H&E staining and immunohistochemistry of serial tissue slides for comparing the expression of LPP protein and the endothelial cell marker CD34. (F) Western blot for confirming the efficacy of LPP knockdown in BJ fibroblast cell line. (G) Subcutaneous GC tumor in mouse from the LPP knockdown group and the control group. For tumor observation, the GC cells were stably transfected with the vector of green fluorescent protein before injection. (H) Bar plot comparing tumor volume in the LPP knockdown group and the control group. Error bar indicates mean + standard deviation.

To investigate the effect of fibroblast-derived LPP expression on chemotherapy response, we evaluated the half-maximal inhibitory concentration (IC50) value of four clinic-used chemotherapy compounds (i.e., cisplatin, docetaxel, gemcitabine, and paclitaxel) in the GDSC database using the TCGA-STAD dataset. The results showed patients with high LPP expression had higher IC50 values of four compounds than those with low LPP expression, suggesting lower sensitivity for these compounds based on chemotherapy in patients with high LPP expression (Figures 7A–7D). By using another dataset, GSE62254, we validated that patients with high LPP expression yielded fewer benefits than those with low LPP expression, whether under the condition of chemotherapy or not (Figures 7E–7H). Overall, these results suggested that GC patients with high LPP expression responded less to current chemotherapy regimens.

Figure 7

High LPP expression associates with fewer benefits from chemotherapy for GC patient

(A–D) Box plot comparing the IC50 value of chemotherapy drugs cisplatin, docetaxel, gemcitabine, and paclitaxel, between the high- and low-LPP patient groups. The boxes represent the median ± 1 quartile, with whiskers extending to the most extreme point within 1.5 interquartile range from the box boundaries. (E and F) Kaplan-Meier plot showing distinct overall survival outcomes between high- and low-LPP patient groups, under condition of receiving chemotherapy or not, by using GSE62254 dataset. (G and H) Kaplan-Meier plot showing distinct disease-free survival outcomes between high- and low-LPP patient groups, under condition of receiving chemotherapy or not, by using GSE62254 dataset.

We then examined the association of LPP expression with immunotherapy response, by using a dataset containing 45 GC patients receiving anti-PD1 therapy (see materials and methods for details). The results showed LPP expression was lower in responsive patients than in those non-responders (Figures 8A and 8B). And patients with high LPP expression yielded a lower response rate than those with low LPP expression (Figure 8C). We also performed similar analyses using another cohort containing 26 melanoma patients receiving anti-PD1 therapy, and similar results were also obtained, although the statistical significance was not reached due to the relatively small sample size (Figures 8D–8F). Further, we calculated the tumor immune dysfunction and exclusion (TIDE) score to predict response to immune checkpoint blockade for GC patients in TCGA-STAD cohort (see materials and methods for details). The results showed that patients with high LPP expression yielded higher TIDE scores than those with low LPP expression (Figure 8G). According to their TIDE scores, patients were then divided into two groups, i.e., responder and non-responder (Figure 8H). We found that response rate in the high-LPP group was lower than that in the low-LPP group (Figure 8I). These findings indicated high expression of LPP could increase the potential of tumor immune escape and resistance to immunotherapy for GC patients.

Figure 8

High LPP expression is associated with resistance to immunotherapy

(A) Box plot showing distinct LPP expression between responder and non-responder after anti-PD1 therapy in 45 GC patients from Kim et al. study. (B) Graph showing distribution of LPP expression for responsive and non-responsive patients from Kim et al. study. (C) Bar plot showing distinct response rates between the high- and low-LPP groups in GC patients from Kim et al. study. (D) Box plot showing distinct LPP expression between responder and non-responder after anti-PD1 therapy in 26 melanoma patients from Hugo et al. study. (E) Graph showing distribution of LPP expression for responsive and non-responsive patients from Hugo et al. study. (F) Bar plot showing differences in response rate between the high- and low-LPP groups using dataset from Hugo et al. study. (G) Box plot showing distinct levels of TIDE score between the low- and high-LPP groups in GC patients from TCGA-STAD cohort. (H) Graph showing distribution of TIDE score for the responder and non-responder groups in TCGA-STAD cohort. (I) Bar plot showing distinct response rates between the high- and low-LPP groups in GC patients from TCGA-STAD cohort. The boxes represent the median ± 1 quartile, with whiskers extending to the most extreme point within 1.5 interquartile range from the box boundaries.

Underlying chemical compounds sensitive to GC patients with high LPP expression

We further explored the underlying chemical compounds sensitive to GC patients with high LPP expression, by using the GDSC database. Based on the IC50 value of 138 compounds in each TCGA-STAD sample, a total of 28 compounds were identified as sensitive drugs for GC patients with high LPP expression (Figure S4). These compounds targeted canonical pathways involved in oncogenesis, such as IGF1R, JNK and p38, PI3K/MTOR, cell cycle, and apoptosis regulation signaling pathway. Our findings could help to explore effective drugs and to establish targeted treatment strategy regarding combination with chemotherapy and immunotherapy against drug resistance.

Discussion

A better understanding of immune microenvironment in tumor and a well-designed tool to evaluate the comprehensive signature of immune microenvironment are of critical importance and urgent need for prognostication of GC patients and for guiding more effective therapy strategies in the clinical setting. In this study, through investigating the association of PD-L1 and IFN-γ expression with tumor microenvironment, we finally identified the fibroblast-derived LPP as an effective predictor for response to chemotherapy and immunotherapy. Our data demonstrated that the LPP expression was an independent risk factor for poor prognosis, and patients with high LPP expression yielded fewer benefits from current chemotherapy and immunotherapy in GC. Finally, we identified 28 compounds as underlying sensitive drugs for those GC patients with high LPP expression.

Although the importance of molecules involving cytoimmunity pathways and bioprocesses such as PD-L1 and IFN-γ has been well-established in predicting immune status and response, the underlying principles for immune escape and its link with tumor-infiltrating immune cells and treatment outcomes remain elusive. Several studies have reported that the impact of tumor-infiltrating immune cells abundances on immunotherapy and chemotherapy responses for patients cannot be underestimated in various malignances.,, For example, it was found that high proportion of CD8+ T cells in tumor tissues could predict a better therapy response and survival outcome for GC patients. In our analyses, we found that GC patients with poor prognosis were characterized by strong expression of PD-L1 and IFN-γ, as well as high infiltrating proportion of activated CD4+ memory T cells in tumor microenvironment, while weak expression of PD-L1/IFN-γ and more resting CD4+ memory T cells were found in those with better prognosis. What is more, patients with high LPP expression whose tumor had more activated CD4+ memory T cell infiltration yielded fewer benefits from chemotherapy and immunotherapy than those with low LPP expression and more resting CD4+ memory T cell infiltration. These findings suggest that chemotherapy and immunotherapy resistance had a closed relation with PD-L1/IFN-γ expression and infiltration by CD4+ memory T cells in tumor microenvironment. The immune escape mediated by immune checkpoint might be partially responsible for this issue. And depicting the principles underlying the complex relationship among the expression of PD-L1/IFN-γ, infiltration and activation of CD4+ memory T cells, as well as immune response might provide a new insight into targeting against therapy resistance and selecting sensitive patients for chemotherapy and immunotherapy.

Due to the complex heterogeneity within and among individual patients, the challenges still remain in aspects of personalized prognostication and management, and follow-up scheduling for GC patients. Attempts have been made to link the molecular characteristics and genetic profiles to specific molecular subtypes of GC, although the elusive heterogeneity has rendered the ensuring data largely inclusive., In the present study, we tried to investigate distinct subtypes of GC from a new insight into tumor-infiltrating immune cells abundances, differing from most previous works based on genetic characteristics for specific molecule type, such as genes, miRNAs, and lncRNAs.14, 15, 16 By integrating seven types of PD-L1- and IFN-γ-associated immune cells, we identified three distinct GC subgroups, which had distinct PD-L1 and IFN-γ expression and survival outcomes. In subsequent analyses, significant differences in genomic characteristics were found across the three subgroups, including expression profile, mutational and epigenetic landscape, tumor mutational burden, as well as interactions of co-occurrence and exclusive mutation events. These findings support previous attempts about integrating the expression of immune checkpoint molecules with those important genetic characteristics for prognostication in multiple malignances. For example, non-small-cell lung carcinoma patients with tumors simultaneously harboring high level of both PD-L1 expression and tumor mutational burden yielded better clinical outcomes than those whose tumor had only one of these features, although there were also inconsistent conclusions from other reports., Through integrating associations between PD-L1/IFN-γ and tumor immune microenvironment, we identified a promising biomarker LPP for patient prognostication and response prediction for chemotherapy and immunotherapy in GC. We found high LPP expression was associated with resistance to chemotherapy and immunotherapy, as well as unfavorable survival outcomes for GC patients. These interesting findings imply that LPP might provide additional information for current clinical decision-making and might be an underlying therapeutic target against resistance to chemotherapy and immunotherapy.

Uncovering the complex crosstalk among various biomolecules underlying tumorigenesis pathology remains an unsettled challenge for researchers. Many systematically and comprehensively computational strategies have been developed to investigate such complex biological interactions. Regulatory network, as one of the effective tools for unraveling regulatory mechanisms underlying the biological process of oncogenesis and cancer progression, has been validated and applied in various malignance entities, including GC. One of the major deliverables from the regulatory network was the hub molecular, defined as the node served as a central role in the network. Hub node identification has been demonstrated as an effective strategy to explore potential candidates for extracting critical components and underlying therapeutic targets from complex regulatory mechanisms. In the present study, by systematically inferring interactions among three types of biomolecules (i.e., lncRNA, miRNA, and gene) whose expressions were significantly dysregulated among the tumor infiltration- and prognosis-distinctive patient subgroups, we constructed an lncRNA-miRNA-gene regulatory network in GC. Based on the network, we finally identified the hub molecular of LPP. LPP has been known as a member of the zyxin family of proteins that regulates cytoskeletal organization and cell adhesion and motility. Several studies have demonstrated that LPP served as a proto-oncogene and its overexpression could promote cell invasion and metastasis in multiple malignances., These previous results were consistent with our findings that LPP was overexpressed in tumor, and high LPP expression was associated with poor prognosis in GC patients. However, to date, the roles of LPP, especially regarding principles of immune response and immune escape in GC, are still unclear. Cecilia et al. found LPP as an endothelial adhesion protein in endothelial cells, which can be regulated by cancer-associated fibroblasts and promote chemoresistance in ovarian cancer. However, these results in ovarian cancer are quite different from our current findings in GC. Our immunohistochemistry data and single-cell RNA-seq analysis consistently revealed that LPP was overexpressed in tumor and was restrictively enriched in the stromal fibroblasts of GC. Further in vivo experiments validated that such fibroblast-derived LPP could promote the growth of subcutaneous tumor in mice. Notably, our subsequent findings suggested that high LPP expression was associated with increased resistance to chemotherapy and immunotherapy for GC patients. How does such fibroblast-derived LPP regulate the drug sensitivity? Our further analyses showed tumors with high LPP expression tended to be infiltrated with more activated and less resting CD4+ memory T cells in the tumor microenvironment. As mentioned earlier, patients with strong PD-L1 and IFN-γ expression had more activated CD4+ memory T cells in the tumor microenvironment and yielded unfavorable prognoses. Collectively, these consistent results demonstrated an underlying mechanism hypothesis that fibroblast-derived LPP associates with the infiltration of activated CD4+ memory T cells in the tumor microenvironment for patients with high PD-L1/IFN-γ expression, which might be contributed by complex intercellular crosstalk, such as chemokines secretion, resulting in immunosuppression. Further efforts are awaited to validate this interesting mechanism, and more evidence is needed to elucidate the complete role of LPP in regulating immune response to chemotherapy and immunotherapy in GC.

There also exist several limitations in this study. First, similar to other analyses applying in silico approaches to infer tumor-infiltrating components, although 22 types of immune cells with their activation states were analyzed in this study, our analyses failed to interpret those immune cell types that were not assessed by the applied algorithm but modulated the immune microenvironment properties of GC. Second, since the immune infiltration and its impact on tumor microenvironment are different in distinct regions of the tumor entity, such as the core of tumor and the invasive margin, it is appropriate to evaluate the variation from different tumor regions accordingly. However, data used for our immune infiltration analysis were derived from bulk tumor tissues, making it impossible to take the location information into account.

Materials and methods

Samples and datasets

Analyzed samples were retrieved from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/) portal and Gene Expression Omnibus (GEO, https://www.ncbi.nlm.gov/geo/) database. We gathered a total of 350 GC samples with coupled transcript expression, somatic mutation, copy number variation, and survival information from the TCGA-STAD project. Four GC patient cohorts were collected from the GEO database, including GSE26899 (n = 93), GSE26901 (n = 109), GSE66229 (n = 300), and GSE28541 (n = 40). A dataset of cancerous and normal samples for TCGA pan-cancer types was collected from the TCGA and Genotype-Tissue Expression Project (GTEx, available at: https://gtexportal.org/). Three additional GEO datasets including GSE29272 (n = 268), GSE62254 (n = 300), and GSE84437 (n = 433) were used for validation analyses. To examine the association of LPP expression with immunotherapy response, cohorts of GC patients (n = 45) and melanoma patients (n = 26), who were receiving anti-PD1 therapy, were collected from Kim et al. and Hugo et al., respectively. Immunohistochemistry results from The Human Protein Atlas database were available at https://www.proteinatlas.org/. The robust multi-array average algorithm was employed to process and normalize the raw data of microarray datasets from Affymetrix platform by using the Affy software package. And those datasets from Illumina were processed using the lumi package. This study was approved by the Ethics Committee of Nanfang Hospital of Southern Medical University.

Estimating abundance of tumor-infiltrating cells

To infer the infiltrating abundance of immune cells in GC samples, the CIBERSORT algorithm was employed based on transcriptomic profiles, running the algorithm using the LM22 signature with 1,000 permutations. LM22 is a gene signature consisting of 547 genes that accurately distinguish 22 mature human hematopoietic populations and activation states, including seven T cell types, naïve and memory B cells, plasma cells, NK cells, and myeloid subsets. Abundances of fibroblasts were estimated by using the Microenvironment Cell Populations-counter.

Identifying PD-L1- and IFN-γ-associated patient subgroups

Spearman correlation was employed to evaluate associations between either PD-L1 (CD274) or IFN-γ (IFNG) expression level and abundances of tumor-infiltrating immune cells. The LASSO, implemented in R package “glmnet,” was applied to prioritize the most relevant immune cell types with PD-L1 or IFN-γ expression level. Then, the abundances of selected relevant immune cells were ln(x+1)-transformed, followed by unsupervised hierarchical clustering analyses. The results were visualized by using R package “pheatmap.”

Somatic mutation analysis

Mutation events analyzed included missense mutations, nonsense mutations, splice-site mutations, frameshift insertions/deletions, and in-frame insertions/deletions. Tumor mutational burden was calculated as the total number of somatic gene coding errors, base substitutions, gene insertion or deletion errors detected per million bases in the whole genome, by using R package “maftools.” The oncoplot was used to visualize the mutation profile by using R package “ComplexHeatmap.” Paired Fisher's exact test was used to identify mutually exclusive and co-occurrence of mutational events among clusters, and a p value of <0.05 was considered significant.

Copy number variation analysis

We employed the GISTIC algorithm to detect the genomic region with significant amplification or deletion for GC samples from distinct patient subgroups. Parameter thresholds for the GISTIC algorithm were set as follows: the amplification or deletion length greater than 98% length of the chromosome arm, with q < 0.05, and other parameters were set as default.

Differential expression analyses of gene, lncRNA, and miRNA

Expression differences of gene, lncRNA, and miRNA were analyzed by using R package “limma.” The statistical threshold was set as the absolute value of log-transformed fold change (|log2FC|) >log21.5, with a false discovery rate (FDR) adjusted p value < 0.05. KEGG pathway enrichment analyses were performed by using R package “Clusterprofiler,” and an adjusted p < 0.05 was considered significant.

Construction of lncRNA-miRNA-gene regulatory network

StarBase database (http://starbase.sysu.edu.cn/) was used to retrieve lncRNA-miRNA interactions. For miRNA-gene interaction, we initially searched all miRNA-gene pairs from three databases, including miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index.php), TargetScan (http://www.targetscan.org/vert_72/), and miRDB (http://mirdb.org/), and then we obtained those conserved in humans at least in two of these three sources. Then, based on these two interaction types (i.e., lncRNA-miRNA and miRNA-gene), three-node regulatory paths of lncRNA-miRNA-gene were constructed. Finally, we constructed the final lncRNA-miRNA-gene regulatory network through extracting those differentially expressed lncRNAs, miRNAs, and genes among distinct subgroups, as well as their interactions. The network topology was explored and the hub gene was identified as the node with the largest value of degree in the network, as proposed by Yu et al.

Risk score generation for prognostication

Univariate Cox regression analysis was employed to select prognostic genes from those with differential expression among distinct subgroups. Those genes yielding a hazard ratio (HR) > 1.2 or < 0.8, with p value < 0.05, were identified as risk or favorable genes, respectively. Then we employed the LASSO-Cox regression algorithm to prioritize the most optimally prognostic candidate genes to establish risk score, by using R package “glmnet.” The final risk score was generated as the sum of the products of the signature gene expression and corresponding LASSO coefficient. The “survminer” package was utilized to determine the optimal cutoff of risk score to stratify patients into low- and high-risk groups.

Single-cell analysis

The used single-cell RNA-seq dataset was available at GSE183904, and 26 primary GC samples containing 119,862 single cells were included. The analysis was performed using R package Seurat 4.0.4. Genes expressed in at least three cells and cells that contained 500 to 6,000 features with less than 20% mitochondria genes were eligible for further analysis. Principal component analysis was performed, and JackStraw package was applied to determine the statistically significant principal component. Cell clusters were identified by a shared nearest neighbor modularity optimization-based clustering algorithm, FindClusters, by setting the resolution to 0.5. Clusters were visualized by Uniform Manifold Approximation and Projection and were annotated by SingleR following manual curation.

Treatment response prediction

Based on the publically available pharmacogenomics database, Genomics of Cancer Drug Sensitivity (GDSC, https://www.cancerrxgene.org/), we used the TCGA-STAD transcriptomic profile to estimate therapeutic response of chemical compounds for each GC patient. The procedure was performed by using the R package “pRRophetic,” in which the IC50 value of the compound was estimated using ridge regression. A lower IC50 value indicated a higher efficacy for the compound to inhibit or kill cancer cells. The thresholds used to identify compounds sensitive to GC patients with high LPP expression were as follows: IC50high–LPP – IC50low–LPP < –1.2, with an FDR-adjusted p < 0.05.

TIDE was a newly developed computational framework to evaluate the potential of tumor immune escape by using transcriptomic profiles., The TIDE score computed for each tumor sample can serve as a surrogate biomarker to predict tumor response to immune checkpoint blockade. We calculated the TIDE score for GC patients, in which high TIDE score indicated increased potential for tumor immune escape and resistance to immunotherapy. An average score was used as cutoff to identify the response to immunotherapy for a GC patient.

Hematoxylin-eosin staining and immunohistochemistry

Hematoxylin-eosin staining and immunohistochemistry were performed as previously described., In brief, for incubation with primary antibodies, tissue slides were incubated at 4°C overnight with anti-rabbit antibody (LPP, 1:250, Abcam, ab126608; CD34, 1:250, Abcam, ab81289). Negative controls were treated identically, but without the primary antibody. Both goat anti-rabbit IgG fluorescent secondary antibody (Aliexa Flour488, 1:400, GB25303) and anti-rabbit horseradish peroxidase secondary antibody (HRP, 1:200, GB23303) were used in our experiments. Percentage positivity (continuous scale) in tissue area was scored in ten random fields, and the mean percentage positivity was used to access the expression level.

In vivo tumor growth assay

Human GC cell line AGS and human BJ fibroblast cell line were purchased from the Cell Resource Center, Shanghai Institute of Biochemistry and Cell Biology at the Chinese Academy of Sciences. The shRNA plasmid of LPP was purchased from the GeneCopoeia (Shanghai, China). AGS cell line with stably expressed green fluorescent protein (AGS-GFP) and BJ fibroblast cell line with stable LPP knockdown as well as its control (LPP-shRNA and Vector) were established. Five-week-old male athymic BALB/c nu/nu mice were purchased from the Central Laboratory of Animal Science at Southern Medical University (Guangzhou, China). The mice were maintained at the Laboratory Animal Centre of Nanfang hospital in a specific pathogen-free environment. For in vivo tumor growth assays, AGS-GFP cells and BJ fibroblasts (LPP-shRNA and Vector) were mixed and subcutaneously injected into the right flank of mice (n = 3 per group). Tumors were measured with calipers every 3 days after injection, and the tumor volumes were calculated as previously reported, V = 0.5 × length ×width2. After 30 days observing, the mice were sacrificed to harvest the tumor for further analysis.

Statistical analysis

All statistical analyses were conducted with R version 3.6.1. Networks were visualized by Cytoscape version 2.8. Kaplan-Meier estimator and log rank test were used for survival comparison. Continuous variables were compared by using Wilcox test, and categorical variables were compared by χ2 (or Fisher's exact test, if appropriate) test. Univariate and multivariate survival analyses were performed by using Cox proportional hazard regression with R package “survival.” A two-tailed p of <0.05 was considered statistically significant.