Identification of genes with therapeutic and prognostic values in lung adenocarcinoma microenvironment

As the most diagnosed malignancy, lung cancer is also the primary cause of cancer death in the entire world. The therapy of lung adenocarcinoma (LUAD), which is the most prevalent subtype of lung cancer, draw researchers’ increasing attentions. This research aimed to investigate the tumor microenvironment (TME)-related hub genes which might be novel targets for treatment. LUAD-associated data packages, including RNA-Seq information and clinical data of 522 patients, were obtained from The Cancer Genome Atlas (TCGA). For better evaluation of stromal and immune cell components, immune scores, stromal scores and estimate scores were obtained with ESTIMATE algorithm based on gene expression levels in tumors. The R package heatmap and clustering analysis were used to explore interested genes. Differentially expressed genes (DEGs) were identified by Venn diagram. Protein-protein interaction (PPI) network was applied to explore intrinsic connections of DEGs. Kaplan-Meier (K-M) survival curves, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis were applied to investigate the prognostic values and intricate biological functions of DEGs. The relationships between 4 survival-related hub genes and 6 types of immune cells were examined using TIMER database. The LinkedOmics database was applied to look for kinase targets of hub genes. A total 702 589 and 113 were identified. GO and KEGG analysis showed the suggested that the top 8 nodes were FPR2, C3AR1, MCHR1, CCR5, FPR1, CCL19, CCR2 and CXCL10. K-M survival curves indicated that FPR2, C3AR1, MCHR1 and CCR5, as hub genes, were significantly correlated with the overall survival (OS) of LUAD patients. The expression levels of C3AR1 and CCR5 were positively correlated with immune cell infiltration. LYN, LCK and SYK were the targeted kinases of these hub genes. FPR2, C3AR1, MCHR1 and CCR5 were TME-related genes and potential biomarkers for the therapy and prognosis of LUAD.


Introduction
Lung cancer is the most frequent malignant carcinoma, and also the major cause of cancer death [1]. Base on histopathologic characteristics, non-small cell lung cancer (NSCLC) accounts for approximately 85% of all lung cancer, and lung adenocarcinoma (LUAD) is the major component of NSCLC [2]. While LUAD patients benefited from multiple treatments, including radiotherapy, chemotherapy and immunotherapy, the 5-year survival rate was still under 20% [3]. Therefore, it is eager to identify novel biomarkers to raise the curative rates, reduce the death rates and improve survival quality.
Tumor microenvironment (TME) composes of cells (endothelial cells, neuroendocrine cell and immune cells, etc.) and extracellular components (cytokines, chemokines and extracellular matrix, etc.). TME has significant effects on the initiation, development, metastasis, and recurrence of carcinoma [4,5]. It was reported that TME affected gene expression, such as SDF1 [6] and CXCR4 [7]. Emerging evidence showed that regulation of TME facilitated chemotherapy and immunotherapy to improve the prognosis of patients [8].
It is vital to find TME-related genes to increase the efficiency of therapies.
As the essential ingredients of TME, immune and stromal cells play substantial effects on signal transduction, immune surveillance and immune scape [9,10]. To accurately evaluate the percentage of immune and stromal cells in tumors, Yoshihara et al. obtained scores for tumor purity, the infiltration levels of immune cells and stromal cells in tumor tissues based on gene expression levels using ESTIMATE algorithm [11]. Precious studies also suggested that TME-related biomarkers were helpful in histological diagnosis and prognostic prediction using big-data based ESTIMATE algorithm [12,13].
As a publicly available database, The Cancer Genome Atlas (TCGA) includes abundant cancer-causing genomic alterations data of more than 30 human tumor types [14]. In the current study, we obtained transcription profiles and clinical information of LUAD patients, including ages, genders, tumor stages, survival states and prognosis from TCGA database.
The immune and stromal scores were calculated, and specific gene expression signatures of immune and stromal cells were analyzed. TME-related prognostic genes were identified, and the correlations between hub genes and immune cell infiltration were explored. The kinase targets of hub genes were also screened.

Data collection and preprocessing
The transcription profiles and clinical data of 522 LUAD patients were obtained from TCGA GDC website (http://portal.gdc.cancer.gov/) [14]. Normalization process was carried out with the R package limma [15]. Using ESTIMATE algorithm, the immune scores, stromal scores and estimate scores were calculated.

Screening of differentially expressed genes
Bioconductor limma package was applied to screen differentially expressed genes (DEGs) from the immune and stromal score groups. |log(FC)| > 1 and False Discovery Rate (FDR) < 0.05 were set as the cutoffs to screen DEGs [15].
Heatmap and clustering analysis R package heatmap and clustering analysis were used to obtain up-and down-regulated genes in the immune and stromal score groups, respectively [16]. Venn diagram package in R software was used to obtain the up-and down-regulated genes in both the immune score group and stromal score group [17]. The total overlapping DEGs were retained for subsequent analysis.

Functional and signal pathways enrichment analysis
The GO and KEGG analysis were applied to explore the function of DEGs. GO analysis included 3 parts: biological processes (BPs), molecular functions (MFs) and cellular components (CCs). To make data visible, 4 R packages, clusterProfiler, org.Hs.eg.db, enrichplot and ggplot2, were applied [18], and the top 10 enrichment clusters were shown in the bar charts. KEGG analysis was shown in the bubble chart. P < 0.05 and q < 0.05 were defined as cut-offs criterion.
Construction of protein-protein interaction network and overall survival curve STRING (https://string-db.org/) was applied to construct protein-protein interaction (PPI) network [19], with the consideration of 0.9 as minimum interaction score. Overall survival (OS)-related hub genes were analyzed with K-M Plotter (http://kmplot.com/analysis/) [20].
Identified the kinase targets of hub genes LinkedOmics database (http://www.linkedomics.org/) was used to explore the kinase targets of hub genes [22]. P < 0.05 was considered as statistical significance.

Data source and clinical characteristics
The 522 LUAD patient's clinical data and profiles of 60483 mRNA expression were downloaded from TCGA database in February 2020. The clinical characteristics of LUAD patients were exhibited in Table 1.  1A). The lower immune scores (P = 0.003) and ESTIMATE scores (P = 0.028), the larger tumor sizes (Fig. 1B). However, there was no statistically significant difference between these scores and lymph node metastases ( Fig. 1C). LUAD patients with distant metastases had lower stromal scores (P = 0.007) and ESTIMATE scores (P = 0.016, Fig. 1D).
According to median immune scores, LUAD patients were separated into low immune score cluster (n = 267) and high immune score cluster (n = 268). Based on median stromal scores, LUAD patients were also separated into low stromal score cluster (n = 267) and high stromal score cluster (n = 268). K-M survival curves showed that the high immune score group had longer OS than the low immune group (P = 0.021, Fig. 2A), and that the high stromal score group had longer OS than the low stromal group (P = 0.103, Fig. 2B). In addition, the high ESTIMATE score group had longer OS than the low ESTIMATE group (P = 0.034, Fig. 2C). These results exhibited that the immune and stromal scores were significantly associated with OS.
Screening DEGs based on immune/stromal scores down-regulated (Fig. 3D). These 702 genes were defined as DEGs (Table S1)  Integrated the interrelation of DEGs and identified significance genes STRING method (https://string-db.org/) was used to construct a PPI network of the 702 DEGs, including 366 nodes and 637 edges, with the average node degree of 3.48 (Fig. 4C).

Different roles of the hub genes in immune cell infiltration
The expression of FPR2, C3AR1, MCHR1 and CCR5 were significantly associated with the OS of LUAD patients (Fig. 5). These 4 genes were defined as hub genes. FPR1, CCL19, CCR2 and CXCL10 were not significantly associated with the OS of LUAD patients ( Figure   S1).
The online program TIMER was used to investigate the correlations between the hub genes and immune cells. Our results suggested positive correlations between the expression of C3AR1 and B cell, CD4 + T cell, CD8 + T cell, macrophage, neutrophil, and dendritic cell infiltration. The expression levels of CCR5 were also significantly correlated with immune cell infiltration. However, FPR2 and MCHR1 were only partially related with immune cell infiltration (Fig. 6).

Kinase targets of hub genes
LinkedOmics database includes miRNA targets, kinase targets and clinical features of 32 cancer types. It was used to explore the kinase targets of hub genes. The top 5 kinase targets of each hub gene were shown in Table 2. LYN, LCK and SYK were mutual kinase targets of hub genes ( Figure S2).

Discussion
Lung cancer draws more attention due to its substantial morbidity and mortality. The noteworthy features of lung cancer are fast development, highly malignant, low rates of early diagnosis, and adverse prognosis, seriously harming human health and wellbeing [1].
LUAD is the most prevalent subtype of lung cancer. During last decades, the discoveries of many therapy targets substantially improved the efficacy, including KRAS, EGFR, BRAF, ALK [23][24][25]. However, LUAD prognosis are still not satisfied. Our research aimed to screen and identify novel targets related to immune infiltration in TME.
The analysis showed that DEGs were mainly enriched in signaling pathways related to B/T cell receptors, CAMs, and Chemokines. Some DEGs were also involved in the transductions of NF-κB signals and Toll-like receptor signals [26]. These results accorded with previous reports on the functions of immune and stromal cells in the TME of LUAD patients [27][28][29].
To analyze the inherent relations of DGEs, PPI analysis was performed by STRING tool, and the results showed that most DEGs participated in immune or inflammatory responses. We obtained the top 8 nodes, including FPR2, C3AR1, MCHR1, CCR5, FPR1, CCL19, CCR2 and CXCL10. Previous study showed that FPR2 deficiency increased secretion of CCL2 and facilitated M2 polarization of macrophages, promoting tumor progression [30]. In addition,

FPR2 enhanced invasion and migration of cancer cells in gastric cancer and colorectal
cancer [31,32]. Mathern et al. found that T cell-expressed C3AR1 amplified the expression of antigen presenting cell costimulatory molecules and innate cytokines, promoting CD8 + T cell proliferation and expansion [33]. Moreover, C3A/C3AR1 signaling facilitated the lung metastasis via activating PI3K/AKT signaling pathway in breast cancer [34]. In gastric cancer, a combined treatment of anti-CCR5 and anti-PD1 improved CD4 + and CD8 + T cell infiltration and reduced the load of mice bearing tumors [34]. The eutopic expression of CCR5 activated calcium signaling and enhanced Treg cell differentiation and infiltration [35]. Previous studies also showed that FPR1 deficient dendritic cells could not elicit anti-tumor T cell immunity, and that FPR1 was required in chemotherapy-induced anti-tumor immunity [36]. As a significant immune-related cytokine, CCL19 enhanced the development of an immune-stimulating intertumoral niche via increasing anti-tumorous CD8 + T cell accumulation in tumor tissues [37]. In hepatocellular carcinoma, knocking out of CCR2 suppressed the tumor progression via inhibiting infiltration and M2 polarization of macrophages, activating CD8 + T cells in TME [38]. CD8 + T cells and macrophages played important roles in the TME and tumor prognosis. In pancreatic tumor tissues, the higher proportions of CD4 + and CD8 + T cells, the better long-term prognosis for patients [40]. In breast carcinoma, annexin1 promoted polarization of M1 macrophages and induced pro-inflammatory genes [41]. Nevertheless, it was reported that the higher numbers of T or B cells, the worse capsular and perineural invasion in prostate cancer [42].
We also proved that LYN, LCK and SYK were the mutual kinase targets of the hub genes. It was reported that down-regulation of LYN inhibited the proliferation, migration and invasion of melanoma cell lines via suppressing PI3K/AKT signaling pathway [43].
Combined treatment of SYK inhibitor and PI3K inhibitor induced anti-tumor immune responses via promoting pro-inflammatory macrophage phenotype and activating CD8 + T cells [44]. LCK was also reported to activate T cell via facilitating costimulatory molecule CD28-mediated signal transduction and LAT phosphorylation [45].
There are several limitations in our current study. First, this study was mainly based on publicly accessible big-data applying the computational-biology algorithm. Second, the intrinsic mechanisms of between these hub-genes and immune cell infiltration remained unclear.
Their kinase targets might be used for immune-combined therapy. Our studies provided new perspectives on LUAD immunotherapy.

Conclusions
Four hub genes, FPR2, C3AR1, CCR5 and MCHR1, were significantly associated with the prognosis and immune cell infiltration in LUAD. LYN, LCK and SYK were the mutual kinase targets of these hub genes. Our researches enhanced TME-related gene features and provided novel immune targets for clinical researches.

Ethics approval and consent to participate
Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
The datasets during and/or analyzed during the current study available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests. Additional Files Figure S1. Association between the OS and FPR1, CCL19, CCR2 and CXCL10 in LUAD patients. Figure S2. Identification of mutual kinase targets of hub genes by Venn diagram. Table S1. Identification of the up-and downregulated DEGs in LUAD. Figure 1 The  The relationships between the immune/stromal scores and OS of LUAD patients.    and CCR5) were showed. P < 0.05 was considered as statistical significance.