Tumor Microenvironment-Associated Immune-Related Genes for the Prognosis of Malignant Pleural Mesothelioma

Malignant pleural mesothelioma (MPM) is a rare but highly aggressive thoracic malignancy. ESTIMATE algorithm-derived immune scores are commonly used to quantify the immune and stromal components in tumors. Thus, this algorithm may help determine the tumor microenvironment (TME)-related gene expression profile associated with tumor immunity. This study aimed at mining public databases to determine a potential correlation between differentially expressed genes (DEGs) and survival in patients with MPM. We categorized patients from the Gene Expression Omnibus database according to their immune/stromal scores into high- and low-score groups. Functional enrichment analysis and the construction of protein–protein interaction networks showed that the DEGs identified were primarily involved in the TME. Furthermore, we validated these genes in an independent cohort of patients with MPM from The Cancer Genome Atlas database. DEG analysis showed that 29 DEGs were cancer driver genes. Subsequently, 14 TME-related genes, which have been previously neglected, were shown to exhibit significant prognostic potential in MPM. In conclusion, immune/stromal scores are novel predictors of a poor prognosis in patients with MPM. We identified DEGs that are involved in immunity against MPM and may contribute to patient survival. Owing to their potential as prognostic factors for MPM, these 14 TME-related genes need to be studied in detail in the future.


INTRODUCTION
Malignant pleural mesothelioma (MPM) is a highly aggressive thoracic malignancy, with three major histological subtypes, namely, epithelioid, sarcomatoid, and biphasic mesotheliomas. MPM is a rare cancer, occurring in 30 out of a million individuals, although its incidence has slightly increased over the last decade (1). Chinese individuals are predominantly affected by MPM, owing to their exposure to high levels of asbestos (2,3). Multimodality treatment (4), including surgical resection, radiation therapy, chemotherapy, and immunotherapy, is preferred for MPM because of the potential to increase efficacy and prognostic value. Thus, a multidisciplinary therapy-based model for the diagnosis and treatment of patients with MPM is extremely important (5)(6)(7).
Immunotherapy has become increasingly important in the management of cancers. Immune checkpoint inhibitors targeting programmed cell death protein-1 (PD-1), its ligand (PD-L1), and cytotoxic T-lymphocyte antigen 4 (CTLA-4) are commonly used as first-line therapeutics against various cancers (5,(8)(9)(10). Several clinical trials have evaluated the efficacy and safety of immunotherapy in patients with MPM (1,3,11,12). Thus, the KEYNOTE-028 trial (1) showed potential clinical benefits of pembrolizumab, an anti-PD-L1 monoclonal antibody, in patients with MPM showing positive PD-L1 expression. Several ongoing phase II trials, including the KEYNOTE-158 basket trial (ClinicalTrials.gov Identifier: NCT02628067) and KEYNOTE-139 study (ClinicalTrials.gov Identifier: NCT02399371), are assessing the activity of pembrolizumab as a second-line therapy for MPM. The potential therapeutic benefit of combination immunotherapy (CTLA-4 and PD-L1 blockade) against MPM is under investigation, and thus far, the data from the CheckMate-743 (ClinicalTrials.gov Identifier: NCT02899299) and NIBIT-MESO-1 trials seem promising (3). There are only a few ongoing trials using a combination of immunotherapy and chemotherapy or surgery, or both, in patients with MPM. The DREAM trial is a multicenter, single-arm, open-label phase II study that aims to determine the effects of durvalumab in combination with chemotherapy for the treatment of MPM. Based on the data presented at the 2018 American Society of Clinical Oncology Conference, the therapeutic efficacy of this regimen seems promising.
The tumor immune microenvironment, including endothelial, stromal, and immune cells, plays a vital role in tumor surveillance and antitumor effects (2,7). Understanding the correlation between tumor immunity and antitumor effectors in the MPM immune microenvironment is critical for enhancing the potential efficacy of immunotherapy (2,7). Increasing evidence suggests that analysis of gene expression or copy numbers in cancer samples helps understand immune cell infiltration into the tumor microenvironment (TME). The ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) algorithm (5) can be used to analyze transcription profiles in tumor samples and their association with tumor cellularity and different infiltration characteristics. This algorithm generates scores to predict the levels of infiltrating stromal and immune cells and tumor heterogeneity (5).
In this study, we used ESTIMATE algorithm-derived immune scores of patients with MPM to generate a list of TME-associated genes based on the Gene Expression Omnibus (GEO) database. Subsequently, genes associated with poor outcomes were validated in a cohort of patients with MPM available from The Cancer Genome Atlas (TCGA) dataset.

Microarray Data Analysis
Gene expression profiles of 55 MPM tumor samples and paired normal tissues (for 41 tumor samples) were obtained from the GEO database under accession number GSE51024 (13). Gene expression was analyzed using the Affymetrix Human Genome U133 Plus 2.0 Array (Affymetrix, Inc., Santa Clara, CA, USA). The raw data were processed using the Robust Multi-Array Average method and the "Oligo" package from BioConductor (http://www.bioconductor.org) to normalize the data and annotate probe information.

The Cancer Genome Atlas Data Analysis
Gene expression was analyzed using the Illumina HiSeq 2000 RNA sequencing platform at the University of North Carolina Genome Characterization Center. Level 3 data were obtained from the TCGA Data Coordinating Center. This dataset showed the predicted transcriptome profiles, presented as log 2 (x + 1)transformed RSEM-normalized counts. Genes were mapped to human genome coordinates using HUGO probeMap. A method description is available from the University of North Carolina TCGA Genome Characterization Center. Demographic and clinical data, such as sex, age, histological type, survival, and outcome, were also obtained from TCGA.

Initial Data Processing and Identification of Differentially Expressed Genes
Normalized signal intensities of the data were imported into BRBArrayTools (v.4.5; http://linus.nci.nih.gov/BRB-ArrayTools. html) for the initial processing and identification of differentially expressed genes (DEGs). A false discovery rate (FDR) of 0.05 was considered significant. DEGs were selected based on a fold change of ≥1.5 and FDR of <0.05. Immune and stromal scores were calculated using the ESTIMATE algorithm (5). Heatmaps were generated, and clustering analyses were performed using the open-source web tool ClustVis (14).

Enrichment Analysis of Differentially Expressed Genes
Functional enrichment analysis of DEGs was performed by DAVID (The Database for Annotation, Visualization and Integrated Discovery) to identify Gene Ontology (GO) categories by their biological processes (BPs), molecular functions (MFs), or cellular components (CCs). The DAVID database was also used to perform pathway enrichment analysis with reference from Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. FDR < 0.05 was used as the cut-off.

Identification of Network-Based Functional Modules
The Microarray Data Analysis tool from ReactomeFIViz was used for network-based functional analysis (24). After using DEGs as the input, we filtered the network modules using a cutoff of ≥2.

Overall Survival Curve
Kaplan-Meier plots were generated to illustrate the relationship between patients' overall survival (OS) and gene expression levels of DEGs. The relationship was tested by log-rank test. p < 0.05 was used as the cut-off.

Study Design and Workflow
In this study, we identified, using a public database, TMErelated genes that might affect the development and progression of MPM. Comparison of patients with high and low immune or stromal scores led to the identification of 74 stromal and immune cell-associated DEGs in the GEO database. GO and KEGG analyses showed that all 74 genes were associated with the TME. Analysis of the DEGs (tumor vs. normal) indicated that 29 DEGs were cancer driver genes, and survival analysis showed that 14 of these 29 genes were associated with a poor prognosis in patients with MPM. The prognostic values of the 14 TME-related genes were validated in patients with MPM from the TCGA database (Figure 1). Furthermore, a PPI network was constructed, in which seven network modules comprising 32 genes were obtained, to better understand the interactions among these DEGs.

Correlation Between Differentially Expressed Genes and Immune and Stromal Scores Using the Gene Expression Omnibus Database
To determine a correlation between DEG profiles and immune and/or stromal scores, we analyzed the data from a transcription microarray of 55 tumor tissues with paired normal tissues (for 41 tumor samples) from the GEO cohort. With the use of the BRBArrayTools software, 2,134 DEGs were identified based on the stromal scores, and 334 DEGs were identified based on the immune scores (Figures 2A-C). A total of 74 DEGs were validated using the Venn algorithm ( Figure 2D).

Significant Association of Immune and Stromal Scores With the Survival of Patients With Malignant Pleural Mesothelioma
The initial pathological diagnoses for 87 patients with MPM, for whom we obtained the gene expression profiles and clinical information from the TCGA database, were made between 1999 and 2013. Demographic and clinicopathological information, including sex, age, tumor location, history of asbestos exposure, histological classification, differentiation grade, pathological T, N, and M stages, and survival, was also retrieved from the  database. However, we observed no correlations between the clinicopathological characteristics and immune and/or stromal scores (Supplementary Figure 1), except for the pathological T stage.
Kaplan-Meier survival curves were used to evaluate the correlation between the immune and/or stromal scores and the prognosis of patients with MPM (Figure 3). The median OS of patients with low stromal scores (43 cases) was higher than that of patients with high stromal scores (43 cases) (719 vs. 414 days, respectively; p = 0.006, by a log-rank test). Similarly, the median relapse-free survival (RFS) of 40 cases with low stromal scores was higher than that of 44 cases with high stromal scores (620 vs. 414 days, respectively; p = 0.022, by a log-rank test). By contrast, the median OS of patients from the high-immune-score group (46 cases) was higher than that of patients from the lowimmune-score group (38 cases) (620 vs. 434 days, respectively; p = 0.015, by a log-rank test). Meanwhile, the median RFS of patients from the high-immune-score group (55 cases) was higher than that of patients from the low-immune-score group (29 cases) (629 vs. 365 days, respectively; p = 0.011, by a logrank test).

Significant Association of Immune Function-Related Genes With the Survival of Patients With Malignant Pleural Mesothelioma
Survival analysis showed the presence of 29 cancer-specific DEGs. Of these, 14 immune-related genes were associated with a poor OS (Figure 4) and RFS (Figure 5), while 15 genes had no significant association with the survival of patients with MPM (Supplementary Figure 2). The 14 immune-related genes ( Table 1) were NT5E, ETS1, C16orf54, WIPF1, FLI1, ARHGEF6, HLA-DRA, PRKCH, C5AR1, SPATA13, RCAN3, NAPSB, SAMHD1, and STX11; and their expression levels are shown in Figure 6A. The relationships between the immune/stromal scores and abnormal expression of these 14 immune-related genes are shown in Figures 6B,C.

Protein-Protein Interaction and Hub Gene Identification
To explore the interplay among the 74 overlapping DEGs, a PPI network, which contained 32 nodes and 42 edges, was

.10E-08 Wilcoxon
The p-value is for the differences in the expression of this gene in the GEO database. Genes in bold is also essential to in protein interaction analysis.
constructed using the Cytoscape software. In the PPI network, HLA-DRA, CD28, ITK, PTPRC, CXCR4, SYK, and PIK3CG were significant nodes because they had the largest numbers of connections with other nodes. Finally, the following nine TMErelated hub genes were identified: HLA-DRA, SYK, CXCR4, ITK, PTPRC, HLA-DPB1, PIK3CG, ETS1, and HCLS1. Functional enrichment analysis using ReactomeFIViz showed that the genes were involved in the C-MYB transcription factor network, MHC class II antigen presentation (R), antigen processing and presentation (K), cell adhesion molecule (CAM) interactions, Fc-epsilon receptor I signaling in mast cells (N), Ctype lectin receptor signaling pathway (K), endothelin signaling pathway (P), T-cell activation (P), and T-cell receptor signaling pathway (K), most of which are related to immunity (Figure 7).

DISCUSSION
It has been previously shown that molecular features of tumors can affect immune responses and the TME (5, 28). Therefore, identifying TME-related genes that determine the prognosis of MPM is crucial to preventing cancer progression and managing the potential immune response. In this study, we used the ESTIMATE algorithm (5) to screen for tumor immune-related genes in the GEO database and validate the prognostic values of genes in the TCGA database. To the best of our knowledge, this is the first report to demonstrate the involvement of 14 immune-related genes in MPM.
The composition of the TME may influence cancer development (29). Interactions between cancer stem cells and immune cells are important in the carcinogenesis of MPM (30). The TME may also be associated with the prognosis of MPM and may determine the potential efficacy of anticancer therapies (31). Immunotherapy exploits these mechanisms and imparts antitumor effects (32). There are multiple reports on the development of novel therapeutic strategies based on targeting the TME (29,33). Recently, encouraging results of a first-line dual-immune phase III clinical trial on MPM were released at the 2020 World Conference on Lung Cancer, which are expected to break the 15-year deadlock without new drugs. The CheckMate-743 phase III clinical study showed that nivolumab injection, combined with ipilimumab, could significantly improve the OS of patients with previously untreated, unresectable MPM. At the shortest follow-up of 22 months, nivolumab, combined with ipilimumab, reduced the patient's risk of death by 26%. The patient's median OS was 18.1 months, compared with 14.1 months in the chemotherapy group (hazard ratio: 0.74; 96.6% confidence interval: 0.60-0.91; p = 0.002). The 2-year survival rate of patients in the nivolumab plus ipilimumab group was 41%, compared with 27% in the chemotherapy group (34). This study suggested that checkpoint inhibitors might be a potential effective treatment for MPM. However, more clinical trials, with OS as a primary end point, investigate the efficacy of immunotherapy compare with the best treatment (pemetrexed + cisplatin ± bevacizumab) for MPM is needed.
In our study, we aimed to explore whether genes related to the tumor immune microenvironment might predict the prognosis of patients with MPM. These genes might help to select the potential beneficiaries of immunotherapy. Five of the 14 immune-related genes (ETS1, FLI1, WIPF1, HLA-DRA, and PRKCH) represented essential elements in the PPI network. ETS1 (also known as p54 and EWSR2) is a member of the ETS family of transcription factors, which are involved in stem cell development, cell senescence and death, and tumorigenesis. ETS1 functions as an oncogene and is a crucial regulator of phenomena involved in tumorigenesis, such as the mesenchymal phenotype in various tumors, including head and neck squamous cell carcinoma (35), breast cancer (36), prostate cancer (37), and glioblastoma (38). ETS1 also functions in several immune-related pathways, including the C-MYB transcription factor network, BCR signaling, HIF-1α transcription factor network, HIF-2α transcription factor network, and IL4-mediated signaling (https://pubchem.ncbi. nlm.nih.gov/source/Pathway%20Interaction%20Database). FLI1 (also known as SIC-1 and BDPLT21) is a transcription factor containing an ETS DNA-binding domain; it is involved in gene fusions (39). FLI1 is associated with a poor prognosis in multiple tumors, including non-small-cell lung cancer (40), breast cancer (41), and acute myeloid leukemia (42). WIPF1 (also known as PRPL-2, WAS2, WASPIP, and WIP) is vital to the organization of the actin cytoskeleton. WIPF1 is involved in Fc-epsilon receptor I signaling, which is an important immune-related pathway, in mast cells. HLA-DRA, an HLA class II alpha-chain paralog, is involved in several immune-related pathways, such as STAT4-mediated IL12 signaling, CXCR4-mediated signaling, TCR signaling in naïve CD4 + T cells, and cytokine inflammatory responses. PRKCH is a protein kinase C and a member of the PKC family of proteins. PRKCH phosphorylates a wide variety of targets and is involved in diverse signaling pathways, including calcium regulation in cardiac cells, G-protein signaling pathways, the Wnt signaling pathway, and pluripotency. Taken together, these genes show prognostic potential and deserve further exploration to determine their potential as therapeutic targets in MPM.
The significance of this study is that we identified 14 TMErelated genes. These genes are involved in immune responses, may predict the survival of patients with MPM, and may also play a role as biomarkers of the sensitivity to immunotherapy. Moreover, we validated these results in two databases. However, because of the low incidence of mesothelioma, it is very difficult to collect fresh clinical samples to verify our results.
Taken together, these identified genes show promise as prognostic markers and immune effectors. To the best of our knowledge, this is the first study to report that immune and stromal scores can predict the prognosis of MPM. High-immune and low-stromal scores are predictive of a good prognosis of MPM in patients. However, most clinicopathological characteristics, including the history of asbestos exposure, histological classification, differentiation grade, and pathological N stage, have no impact on the prognosis of MPM.

CONCLUSION
In summary, our findings indicate the importance of TMErelated genes that are involved in immune-related pathways in the prognosis of MPM and serve as potential predictors to improve the efficacy of precision immunotherapy. Based on the ESTIMATE algorithm, immune and stromal scores were calculated and found to be useful in determining the prognosis of patients with MPM. The DEGs were validated in two independent cohorts, from the GEO and TCGA databases. As the involvement of these genes in MPM was not confirmed in clinical samples from China, further studies are warranted in the future to confirm their roles in MPM.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
XX and WM guaranteed the integrity of the manuscript and were involved in the conceptualization and design of the study. YF and LC contributed to the data collection, analysis, and interpretation. All the authors contributed to writing and revising the manuscript. The final manuscript has been read and approved by all authors.