Characterization of Tumor Microenvironment in Lung Adenocarcinoma Identifies Immune Signatures to Predict Clinical Outcomes and Therapeutic Responses

Background and Objective Increasing evidence has elucidated the clinicopathological significance of individual TME component in predicting outcomes and immunotherapeutic efficacy in lung adenocarcinoma (LUAD). Therefore, we aimed to investigate whether comprehensive TME-based signatures could predict patient survival and therapeutic responses in LUAD, and to assess the associations among TME signatures, single nucleotide variations and clinicopathological characteristics. Methods In this study, we comprehensively estimated the TME infiltration patterns of 493 LUAD patients and systematically correlated the TME phenotypes with genomic characteristics and clinicopathological features of LUADs using two proposed computational algorithms. A TMEscore was then developed based on the TME signature genes, and its prognostic value was validated in different datasets. Bioinformatics analysis was used to evaluate the efficacy of the TMEscore in predicting responses to immunotherapy and chemotherapy. Results Three TME subtypes were identified with no prognostic significance exhibited. Among them, naïve B cells accounted for the majority in TMEcluster1, while M2 TAMs and M0 TAMs took the largest proportion in TMEcluster2 and TMEcluster3, respectively. A total of 3395 DEGs among the three TME clusters were determined, among which 217 TME signature genes were identified. Interestingly, these signature genes were mainly involved in T cell activation, lymphocyte proliferation and mononuclear cell proliferation. With somatic variations and tumor mutation burden (TMB) of the LUAD samples characterized, a genomic landscape of the LUADs was thereby established to visualize the relationships among the TMEscore, mutation spectra and clinicopathological profiles. In addition, the TMEscore was identified as not only a prognosticator for long-term survival in different datasets, but also a predictive biomarker for the responses to immune checkpoint blockade (ICB) and chemotherapeutic agents. Furthermore, the TMEscore exhibited greater accuracy than other conventional biomarkers including TMB and microsatellite instability in predicting immunotherapeutic response (p < 0.001). Conclusion In conclusion, our present study depicted a comprehensive landscape of the TME signatures in LUADs. Meanwhile, the TMEscore was proved to be a promising predictor of patient survival and therapeutic responses in LUADs, which might be helpful to the future administration of personalized adjuvant therapy.


INTRODUCTION
Lung adenocarcinoma (LUAD) is the commonest histological type of all lung cancers, accounting for approximately 50% of them (1,2). Nowadays, surgical resection remains as the standard treatment for early-stage LUADs. Meanwhile, chemotherapy plays an important role in LUAD patients at all stages of the disease. Recently, blockade of immune checkpoint has been proved as a promising therapy for patients with LUADs (3,4) which seems to be an alternative to conventional chemotherapy. Nevertheless, the fact that most patients do not derive any benefit from PD-1/PD-L1 blockade combined with the risk of serious immune-related adverse events (5,6) and the significant up-front costs, underscore the need for developing accurate tools for predicting therapeutic response to chemotherapy and immune checkpoint blockade (ICB) (7).
Nowadays, an increasing body of literature suggests a crucial role for the tumor microenvironment (TME) in cancer progression and therapeutic responses (8)(9)(10)(11). The TME context in LUADs has been reported not only to reflect the potential benefits from treatment (12)(13)(14), but also to predict patient survival (15,16), which includes tumor-infiltrating lymphocytes (e.g., CD8+ T cells, CD4+ T cells), tumorassociated macrophages (TAMs), cancer-associated fibroblasts (CAFs), and other cell types. With the introduction of computational methods to assess the abundance of cells infiltrating in the TME, several studies using these methodologies have explored the clinical utility of TME context (13,(17)(18)(19). However, the comprehensive landscape of TME infiltrates and its predictive power for therapeutic responses in LUADs have not been fully investigated.
In our present study, two previously proposed computational algorithms (20,21) were employed to estimate the fractions of 23 immune and stromal cells based on LUAD gene expression profiles from The Cancer Genome Atlas (TCGA) database. The TME infiltrating patterns of LUAD samples were investigated and correlated with both transcriptomic characteristics and clinicopathological features. Unsupervised clustering was applied to quantify the TME infiltrating patterns which were calculated as a TMEscore. Consequently, the TMEscore was proved to be a promising prognostic biomarker and a robust predictive factor for the therapeutic responses in LUADs.

LUAD Datasets and Preprocessing
The transcriptomic dataset of LUAD from TCGA database was downloaded from the UCSC Xena browser (https://xenabrowser. net/datapages/). LUAD patients without survival information were removed from further evaluation, among whom 499 were available to construct the TMEscore. Data of somatic mutations (MuSE Variant Aggregation and Masking) were downloaded from TCGA database, which included 567 LUAD specimens. Somatic mutation data, transcriptomic data and survival information were available in 493 of the 499 specimens, of which clinical characteristics were accessible in 474. The raw data in TCGA dataset generated from Illumina were processed using the lumi software package according to a previous study (8).
The microarray data (GSE68465) generated by Affymetrix were obtained from the Gene Expression Omnibus (https://www.ncbi. nlm.nih.gov/geo/) (GEO) as a validation dataset, which included 422 LUAD specimens with expression profiles and clinical outcomes available. The raw data for the dataset from Affymetrix were processed using the Range Migration Algorithm (RMA) for background adjustment in the Affy software package (22). The RMA was used to perform background adjustment, quantile normalization, and final summarization of oligonucleotides per transcript using the median polish algorithm. The detailed information of the LUAD datasets is listed in Supplementary  Table 1.

Assessment of Infiltrating Cells in TME
To calculate the proportions of immune cells in the LUAD samples, we used the CIBERSORT (cell type identification by estimating relative subset of known RNA transcripts) algorithm (20) and the LM22 gene signature, which allows for highly sensitive and specific discrimination of 22 human immune cell phenotypes (8). CIBERSORT is a deconvolution algorithm that uses a set of reference gene expression values (a signature with 547 genes) considered a minimal representation for each cell type and, based on those values, infers different cell type proportions in tumor samples using support vector regression (8). Gene expression profiles were prepared using standard annotation files, and data were uploaded to the CIBERSORT web portal (http://cibersort.stanford.edu/), with the algorithm run using the LM22 signature and 1,000 permutations. Proportions of stromal cells were estimated by applying the Microenvironment Cell Populations (MCP)-counter method, which allows for robust quantification of the absolute abundance of immune and stromal cell populations in heterogeneous tissues from transcriptomic data (21).

Consensus Clustering for TME-Infiltrating Cells
Hierarchical agglomerative clustering (based on Euclidean distance and Ward's linkage) was employed to group the samples with qualitatively different TME cell infiltration patterns. Unsupervised clustering methods (K-means) (23) for dataset analysis were used to identify TME patterns and classify patients for further analysis. A consensus clustering algorithm was applied to determine the number of clusters (8), which was repeated 1000 times to ensure the stability of classification using the ConsensuClusterPlus R package (24).

Generation and Analysis of TME Gene Signatures
To identify genes associated with TME cell infiltrating patterns, we classified patients into three TMEclusters. Differentially expressed genes (DEGs) among these three groups were determined using the R package limma (25). DEGs among different TME patterns were determined by significance criteria (adjusted p value < 0.05; |logFC|>0.58). An unsupervised clustering method (K-means) for analysis of DEGs was employed to classify patients into two groups for further analysis. TME signature genes were then obtained using the random forest classification algorithm to screen redundant genes. Gene-annotation enrichment analysis using the clusterProfiler R package (26) was performed on TME signature genes. Gene Ontology (GO) terms were identified with a strict cutoff of p < 0.01 and false discovery rate (FDR) of less than 0.05.

Establishment of TME Scores
Cox regression model was applied to assess the prognostic value of each signature gene which was classified according to its Cox coefficient. A method similar to gene expression grade index (24) was used to define the TMEscore of each patient: where X is the expression level of genes whose Cox coefficient is positive, and Y is the expression level of genes whose Cox coefficient is negative. The cut-off values of each dataset were evaluated based on the association between patient overall survival (OS) and TMEscore in each separate dataset using the survminer package (8). The R package MaxStat (27) which iteratively tests all possible cut points to find the one achieving the maximum rank statistic, was used to dichotomize TMEscore, and patients were then divided into low-and high-TMEscore subgroups.

Analysis of Tumor Mutation Profiles
Somatic mutation data were obtained from the publicly available TCGA database. Notably, one sample with merely silent mutation was excluded from the aforementioned 493 samples in our analysis. We prepared the Mutation Annotation Format (MAF) of somatic variants, and implemented the R package Maftools (https://bioconductor.org/packages/release/bioc/html/ maftools.html) which provides a multiple of analysis modules to perform the visualization process (28) to display somatic landscape. In addition, the R package SomaticSignatures (https://bioconductor.org/packages/release/bioc/html/ SomaticSignatures.html) was used to characterize the mutation signatures of the LUAD samples (29). Mutational signatures were extracted using 96 nonnegative components (single-base somatic substitutions and their immediate sequence context) and compared to the validated consensus mutational signatures in the Catalogue Of Somatic Mutations In Cancer (COSMIC) (30), version 2 (https://cancer.sanger.ac.uk/cosmic/signatures_v2) to identify the set of COSMIC mutational signatures in TCGA datasets (31). Moreover, the estimation of TMB in LUAD samples was conducted according to a previous study (32).

Predictive Value of TMEscore to Estimate Therapeutic Effect
Tumor Immune Dysfunction and Exclusion (TIDE) (http://tide. dfci.harvard.edu/) (33), a computational method to predict ICB response based on melanoma patients who underwent anti-PD-1 or anti-CTLA-4 agent, was used to investigate the predictive value of TMEscore for immunotherapy. TIDE uses a set of gene expression markers to estimate two distinct mechanisms of tumor immune evasion, including dysfunction of tumor infiltrating cytotoxic T lymphocytes (CTL) and exclusion of CTL by immunosuppressive factors (34). Patients with higher TIDE score have a higher chance of antitumor immune escape, thus exhibiting lower response rate of ICB treatment (33). The TIDE score was shown to have a higher accuracy than PD-L1 expression level and tumor mutation burden (TMB) in predicting survival outcome of cancer patients treated with ICB agents (34)(35)(36)(37). The R package MaxStat (27) was also employed to dichotomize the TMB level. The R package pRRophetic (38) was used to determine whether TMEscore could accurately predict clinical chemotherapeutic responses.

Immunohistochemistry Staining
LUAD samples resected from a cohort of chemo-and/or radionaïve patients (Supplementary Table 2) were obtained from the Second Affiliated Hospital of Soochow University, which was approved by the Institutional Review Board (IRB NO.JD-HG-2020-09). The chairperson of the ethics committee waived the need for patient consent. The sections of tumor tissues were firstly deparaffinized and rehydrated. Endogenous peroxidase was then quenched using 10% H2O2 for 10 min at room temperature. Subsequently, nonspecific proteins were blocked with 10% goat serum for 1 h. Afterwards, the sections were rinsed and incubated with anti-BTK (YM0083, Immunoway; diluted 1: 400) overnight at 4°C. The DAB Horseradish Peroxidase Color Development Kit (Beyotime, China) was used for color development. Finally, the sections were counterstained with hematoxylin and mounted.

Statistical Analysis
For comparisons of two subgroups, unpaired Student t tests was used to estimate statistical significance for normally distributed variables, and Wilcoxon rank-sum test was used for analyzing nonnormally distributed variables. To identify significant genes in the DEG analysis, Benjamini-Hochberg method was applied to converting the p values to FDRs (40). The Kaplan-Meier method was used to generate survival curves for the subgroups in each dataset, and the Log-rank test was used to determine the statistical significance of differences. The hazard ratios for univariate analyses were calculated using a univariate Cox proportional hazards regression model. A multivariate Cox regression model was used to determine independent prognostic factors. The R package pROC (41) was used to plot and visualize receiver operating characteristic (ROC) curves to calculate the area under the curve (AUC) and confidence intervals to evaluate the diagnostic accuracy of TMB and TMEscore. For comparison of AUCs, likelihood ratio test for two correlated ROC curves was used. All statistical analyses were conducted using R (https://www.r-project.org/) or SPSS software (version 25.0). A two-tailed p-value < 0.05 was considered statistically significant.

Characterization of TME in LUADs and Distinct Patterns of TME Subtypes
The general flowchart of our study is shown in Supplementary Figure  their transcriptomic data were included in our initial analysis. By applying CIBERSORT algorithm and MCP-counter method, we obtained the proportions of 23 different immune and stromal cells in the 499 samples (Supplementary Figure 1B). Meanwhile, a TME cell network was used to depict the comprehensive landscape of tumorimmune cell interactions and their effects on the OS of patients with LUADs ( Figure 1A and Supplementary Table 3).
Unsupervised learning using K-means algorithm was used on the dataset to perform group clustering, which identified K=3 according to the elbow method and gap statistic ( Supplementary Figures 2A, B). To identify the aforementioned optimal cluster number, we assessed clustering stability using the ConsensusClusterPlus R package which displayed the clustering stability using 1,000 iterations of hierarchical clustering. The consensus matrix supported the existence of three robust clusters of LUADs (Supplementary Figure 2C). As shown in Figures 1B, C, the proportions of infiltrating immune cells and histologic subtypes differ significantly among the three TME subtypes. We found that naïve B cells accounted for the majority in TMEcluster1, while M2 TAMs and M0 TAMs took the largest proportion in TMEcluster2 and TMEcluster3, respectively. However, log-rank test revealed no significant difference in survival among different TMEclusters (p = 0.45) (Supplementary Figure 2D).

Construction and Validation of the TMEscore in Different LUAD Datasets
A total of 3395 DEGs among the three TME clusters were determined by significance criteria (adjusted p value < 0.05; |logFC| > 0.58) as implemented in the R package limma (Supplementary Figure S3A). An unsupervised clustering method (K-means) for analysis of DEGs was then employed to classify the patients into two groups (Supplementary Figure  3B). Among these DEGs, 217 TME signature genes were obtained using the random forest classification algorithm, on which Gene-annotation enrichment analysis using the clusterProfiler R package (26) was performed. Consequently, the data indicated that these TME signature genes significantly enriched in pathways associated with T cell activation, lymphocyte proliferation and mononuclear cell proliferation (Supplementary Figure 3C).
Cox regression model was used to assess the prognostic value of each signature genes according to the Cox coefficient, by which the TMEscore was established for each patient. As shown in Figure 2A, the TMEscore could effectively distinguish significantly different OS in the entire cohort. Notably, patients with high TMEscore (n = 120) had significantly better survival than those with low TMEscore (n = 354) (p < 0.0001). Additionally, we visualized the relationships of TME subtypes and TMEscore subgroups as well as patient outcomes using an alluvial diagram ( Figure 2B). Moreover, TMEscore remained efficient in stratifying the patients with early-stage disease (stage I-II) into different OS (Supplementary Figures 4A, B), which exhibited its value in the external validation cohort from the GEO database as well (Supplementary Figure 4C). To be noted, even for smokers in the validation cohort, significantly different OS was also observed among the subgroups stratified by the TMEscore (Supplementary Figure 4D). A forest plot was used to summarize the predictive value of TMEscore in different patient cohorts ( Figure 2C).

Association Between TMEscore and Cancer Somatic Genomes
We analyzed the somatic variations of the 492 LUAD samples in TCGA database which revealed missense mutation as the leading type of single nucleotide variant (SNV) (Supplementary Figure  5A). Meanwhile, SNV was identified as the major variant in LUADs that occurred more frequently than insertion or deletion (Supplementary Figure 5B). In addition, it was observed that C>A was the predominant SNV type in LUADs (Supplementary Figure 5C). Besides, we calculated the TMB and showed the mutation type with different colors in LUAD samples (Supplementary Figures 5D, E), as well as the top 10 mutated genes in LUADs with ranked percentages (Supplementary Figure 5F). We also mapped the landscape of mutation profiles among the three TMEclusters, which characterized the mutation types of the frequently mutated genes (Supplementary Figure 6). Waterfall plots were then used to exhibit the mutation profiles of patients with high/low-TMEscore in which various colors with annotations at the bottom represented the different mutation types and TMB levels ( Figures 3A, B). Meanwhile, boxplots were applied to showing the mutation frequency of each frequently mutated genes in high-and low-TMEscore subgroups, respectively ( Figures 3C, D). classification (30). This 96 substitution classification is particularly useful for distinguishing mutational signatures that cause the same substitutions but in different sequence contexts (30). A previous study (30) applied this approach to the 30 cancer types and revealed 21 distinct validated mutational signatures. Each mutational signature was characterized by different substitutions and was associated with epidemiological and biological features of particular cancer types. Therefore the frequency distribution of the 96 mutations based on the six classes of base substitution were analyzed in both TMEscore subgroups ( Figures 4A, B). Moreover, by comparing the extracted signatures from our samples with those in COSMIC, we identified different mutational signatures between the two TMEscore subgroups. The identified signatures in both groups similarly showed a strong resemblance to COSMIC signature 13 and COSMIC signature 4 ( Figures 4C, D). A high similarity to COSMIC signature 7 was also found in the high-TMEscore subgroup ( Figure 4C). Meanwhile, a mutational signature similar to COSMIC signature 5 was present in the low-TMEscore subgroup ( Figure 4D). As reported previously (30), COSMIC signature 13 and COSMIC signature 4 were associated with APOBEC and smoking, respectively, while COSMIC signature 7 was mainly associated with ultraviolet light. A genomic landscape of the LUADs was thereby plotted by integrating the TMB information and clinical characteristics with TME profiles including the 217 TME signature genes and TMEscore ( Figure 4E).
Predictive Value of TMEscore as a Biomarker for Therapeutic Effect TIDE to evaluate TMEscore as a predictor of immunotherapy, interestingly, no significant difference (p = 0.8) was observed in TIDE score between the high-(n = 120) and low-TMEscore (n = 354) subgroups ( Figure 5A). However, we observed significantly different PD-L1 expression levels between the high-and low-TMEscore subgroups (p < 0.001) ( Figure 5B). Since microsatellite instability (MSI), the spontaneous loss or gain of nucleotides from repetitive DNA tracts, is a promising predictive biomarker for patient survival and response to immunotherapy (42,43), TMEscore was compared between the two subgroups stratified by MSIscore which was calculated by TIDE (44). As shown in Figure 5C, TMEscore in the high-MSI subgroup (n=314) was significantly higher than that in the low-MSI one (n = 160) (p < 0.001). Meanwhile, TMEscore in the high-TMB subgroup (n = 257) was significantly lower than that in the low-TMB one (n = 217) (p < 0.001) ( Figure 5D). Furthermore, the AUC indicated that TMEscore was superior to TMB alone in predicting response to ICB (p < 0.001) which exhibited greater accuracy in combination with TMB ( Figure  5E). Additionally, comparison of the 50% inhibitory concentration (IC50) of chemotherapy drugs indicated that the low-TMEscore subgroup had higher sensitivity to cisplatin (p = 0.019) ( Figure 6A) while the high-TMEscore subgroup was prone to the benefits from gemcitabine (p = 0.00042) ( Figure 6B).

DISCUSSION
Our findings indicated that assessment of the immune and stromal statuses via the TME signature provided a potent predictor of survival in early-stage patients with LUAD and a promising biomarker for therapeutic responses as well. Based on the DEGs and GO enrichment analysis, we observed that the TME signature genes significantly enriched in pathways mainly associated with lymphocyte activation and proliferation. Moreover, missense mutation was identified as the leading type of SNV which was identified as the major variant in LUADs. Meanwhile, the established TMEscore could stratify the patient cohort from TCGA database into two subgroups with distinct mutation profiles and COSMIC signatures. A genomic landscape of the LUADs was also characterized by integrating the TMB information and clinical characteristics with TME profiles. To our best knowledge, the present study is the first bioinformaticsbased study to comprehensively investigate the associations among clinicopathological features, the mutation spectra and the TME profiles of LUADs, which also developed a computational algorithm and demonstrated its predictive value for both ICB and chemotherapeutic agents. Previously, Yue et al. (45) identified a prognostic gene signature associated with TME, and validated its predictive accuracy for OS in LUAD patients. However, they omitted to analyze whether the signature could predict therapeutic responses to ICB or chemotherapy. In their study, univariate Cox regression analysis was initially employed to screen out 23 prognostic TME-related DEGs. Afterwards, Yue et al. (45) used selection operator (LASSO) and multivariate Cox regression analyses to identify three key genes for constructing a prognostic model, in which the analytical methods were different from ours. Interestingly, there are five common DEGs between their 23 TME-related DEGs and our 217 signature genes, including BTK, CCDC69, CD33, CD52, and LY86. Among the five common DEGs, only BTK was selected out as a signature gene in both their and our studies. We therefore validated the prognostic impacts of the five DEGs using the LUAD cohort from TCGA database (Supplementary Figures 7A-E). Additionally, the role of BTK as a prognostic factor was also validated in our domestic cohort using IHC (Supplementary Figure 7F). Similar to our data, Bi et al. (46) who performed a study based on TCGA mining also found that BTK was an immune-related gene and a promising prognostic factor for LUAD. Recently, Tan et al. (47) performed a bioinformatics study to characterize the immune landscape of LUADs, in which they divided the patients with LUAD into two immunophenotypes based on the tumor microenvironment. The two immunophenotypes were denoted as the Active Immune Class and Exhausted Immune Class. The former showed significant IFN, T-cells, M1 macrophage signatures, and better prognosis, while the latter presented an exhausted immune response with activated stromal enrichment, M2 macrophage signatures, and immunosuppressive factors such as WNT/ transforming growth factor-b. In addition, Tan and his colleagues (47) identified their developed Immune Class as a useful tool to predict the response to ICB. However, merely 32 patients in a metastatic melanoma cohort were included in their prediction of PD-L1 inhibitor response. In addition, the relationships of clinicopathological characteristics and immunophenotypes were not comprehensively characterized, neither were the somatic landscape and COSMIC signatures. Additionally, Huang et al. (48) identified two LUAD subtypes with specific immune and metabolic state based on the bioinformatic analyses of TME, who constructed a TME score to predict TME phenotypes in LUADs. Interestingly, they constructed the TME score using the principal component analysis algorithm based on the TME signature genes, which could represent the signature of the two TME clusters in their study (48). Differently, their TME score was applied to evaluating the expression patterns of immune-associated genes in LUADs, which merely exhibited diagnostic value. It remained unknown whether the TME score raised by Huang et al (48) could accurately predict patient survival and therapeutic responses to ICB or chemotherapy. Meanwhile, the data of somatic variations in LUADs were unavailable in their study. Hitherto, several bioinformatics studies (34,47,49,50) have identified the immune-related signatures as a prognostic biomarker in LUADs. However, the roles of immunophenotypederived signatures in predicting response to chemotherapy have not been fully clarified. Notably, a number of clinical investigations (51)(52)(53) have highlighted adjuvant chemotherapy as a prognostic factor for improved survival in patients with stage IB LUADs. Under the circumstances, our study not only proved TMEscore as a survival-related predictor, but also identified its potential in stratifying patients with distinct sensitivity to different chemotherapy regimens. According to our analysis, low-TMEscore subgroup had higher sensitivity to cisplatin while high-TMEscore group was more likely to respond to gemcitabine, which offered insights into the administration of personalized adjuvant therapy.
There are some limitations that should not be ignored in our study: 1) lack of domestic sequencing data to validate the associations between the TME infiltrating patterns and clinicopathological characteristics as well as mutation spectra; 2) lack of external validation cohort to confirm the roles of TMEscore in predicting therapeutic responses; 3) possessing predictive power though, the chosen set of signature genes are not necessarily valid in terms of their biological significance. For instance, the differences in IC50 scores between high-and low-TMEscore subgroups in Figure 6 do not seem to be obvious although statistically significant differences were obtained. It raises the possibility that the observed differences may not be biologically important; however, they are statistically significant due to a large number of samples included in the analyses. More prospective clinical trials are warranted to verify the potentials of TMEscore in predicting patient outcomes and response to both chemotherapeutic agents and ICB.
In conclusion, our present study depicted a comprehensive landscape of the TME signatures in LUADs. Meanwhile, the TMEscore was proved to be a promising predictor of patient survival and therapeutic responses in LUADs, which might be helpful to the future administration of personalized adjuvant therapy.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
Ethical review and approval were not required for the study on human participants in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required for this study in accordance with the national legislation and the institutional requirements.

AUTHOR CONTRIBUTIONS
DC, YW, XZ, and QD: conceptualization, formal analysis, and writing original draft. XW, YX, and WW: data curation, methodology, and formal analysis. YM: software, methodology, supervision, and reviewing the draft. YC and CC: supervision, reviewing, and editing the draft. All authors contributed to the article and approved the submitted version.