Identification of novel immunomodulators in lung squamous cell carcinoma based on transcriptomic data

1.   Introduction

Lung cancer is one of the most lethal cancers in the world, and it is classified into two kinds: small cell lung cancer (SCLC) and non-small cell lung cancer (NSCLC). NSCLC accounts for approximately 85% of lung cancer cases and is mainly divided into two categories based on etiology and histological pattern: lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) ^[1,2]. Lung cancer is treated with various treatments, including surgical resection, chemotherapy, radiotherapy, targeted therapy, and immunotherapy ^[3]. Surgical resection is used to treat patients with early-stage NSCLC. Because the tumor cannot be physically removed, the best therapeutic option for patients in advanced stages of NSCLC is targeted therapy or immunotherapy coupled with chemotherapy ^[4].

Since inhibitors of epidermal growth factor receptor (EGFR), anaplastic lymphoma kinase (ALK), and vascular endothelial growth factor (VEGF) became available in the early twenty-first century, patients with advanced LUAD had a better prognosis ^[5]. Patients with LUSC, on the other hand, did not have similar increases in survival; this might be because most lung cancer therapy breakthroughs in the last decade have improved the prognosis for adenocarcinoma, but not for LUSC ^[6]. The current slate of options for squamous lung cancer therapy is limited.

Understanding and researching possible molecular targets for LUSC might identify new methods for treating this kind of lung cancer. After immune checkpoint therapy (ICB), a type of therapy that utilizes checkpoint inhibitors to allow the immune system to recognize and attack tumor cells, some LUSC patients continue to develop further, and treatment results are poor ^[7]. Some cytokines, such as IL-6, enable tumor cells to evade immune monitoring by promoting the stemness of LUSC cells in the tumor microenvironment ^[8]. Evasion of immune surveillance and chronic inflammation are hallmarks of tumor growth in all tumor tissues, including LUSC ^[9]. In this study, we investigated differentially-expressed genes (DEGs) in LUSC. After analyzing these genes and their pathways, a prognosis-related core gene set was constructed, and survival models were built using multifactorial Cox regression. This resulted in a better understanding of occurrence and development processes of LUSC and its prognosis. We selected IGLC7, the survival model's key gene, for further investigation.

To date, the mechanisms underlying the oncogenic role of IGLC7 in lung cancer remain mostly unknown. In our study, we found that IGLC7 is involved in the control of several immune-related pathways, suggesting its possible role in the immune microenvironment. However, to date, no mechanism has been reported regarding the immunological role of IGLC7 in cancer. Here, we explored the association of IGLC7 with immune cells and analyzed how it is involved in immune regulatory mechanisms. Our findings demonstrate a role for IGLC7 in tumor immunity, suggesting that IGLC7 may be a potential immunotherapeutic target for LUSC.

2.   Materials and methods

2.1.   Data collection

The Cancer Genome Atlas (TCGA) (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga) and the University of California, SantaCruz (UCSC) Xena (https://xenabrowser.net/datapages/) databases were used to obtain gene expression profiles and clinical information for LUSC patients, including gender, age, histopathological type, and survival data. We used the ESTIMATE algorithm to calculate the immune/stromal scores and the ESTIMATE scores for LUSC patients (https://bioinformatics.mdanderson.org/estimate/rpackage.html). The immune cell marker gene set used for Gene Set Variation Analysis (GSVA) was obtained from Gabriela et al. ^[10].

2.2.   Identification of immune DEGs

DEG calculation was conducted using R software (version 3.6.0; https://www.r-project.org/) and the limma package ^[11]. The cutoffs for screening of immune DEGs were absolute fold change ≥ 1.5 and adjusted p-value < 0.05.

2.3.   Bioinformatics analysis

The Cluster Profiler R package was used for gene ontology (GO) analysis of genes associated with the immune response. KEGG (Kyoto Encyclopedia of Genes and Genomes) analysis was also performed to identify enriched pathways. The Microenvironment Cell Populations-counter (MCP-counter) method was used to calculate the absolute abundance of immunological cell populations. Scores of gene sets associated with immunological functions and inflammatory activities were calculated using GSVA.

2.4.   Statistical analysis

The R programming language was used for all statistical testing. The correlations between continuous variables were estimated using the Wilcoxon and Kruskal-Wallis tests. Cox proportional hazards model analysis was used to assess IGLC7's prognostic value. The R packages ggplot2, pheatmap, survival ROC, circlize, and corrplot were used to further statistical computations and graphical work, and a statistically significant difference was defined as p-value < 0.05.

2.5.   Survival model construction

We use the R package glmnet to perform lasso cox regression for the analysis. First, we analyze the trajectory of each independent variable, and in the training set we perform 1000 lasso regression analyses, followed by model construction using 10-fold cross-validation, and then aggregate the results of each downscaling, counting the number of occurrences of each probe in 100 times, where the combination of the maximum frequency of occurrence can be observed.

3.   Results

3.1.   Correlation of immune, stromal, and ESTIMATE scores with clinical subgroup information in lung squamous cell carcinoma

The TCGA and UCSC Xena (https://xenabrowser.net/datapages/) databases were used to obtain gene expression profiles, demographics, and clinical characteristics of 496 LUCS samples. After analysis, we found that the 496 LUSC cases had immunological scores distributed between-1646.821 and 3063.732, stromal scores between-2286.256 and 1799.134, and ESTIMATE scores between-3414.484 and 4445.312. Separating results based on tumor staging, the three scores for stages I-III became concentrated relative to stage IV, although it was not statistically significant (n.s., not significant; *p < 0.05, ***p < 0.01). For the T, N, and M stages, no significant differences were seen across stages. However, immune scores and stromal scores significantly differed between the sexes, suggesting the potential significance of immune, stromal, and ESTIMATE scores in differentiating clinical features (Figure 1).

In addition, we assessed the potential association between overall survival and immune score. Using the immune score obtained from the ESTIMATE algorithm, we divided the 496 LUSC cases with both immune score and survival data into two subgroups: high a (score ≥ 0) and low (score < 0). Kaplan-Meier curves showed that between these two groups, the prognosis of the high-scoring group was worse than that of the low-scoring group, and the difference was significant (p = 0.039); This result suggests that the immune score enables a certain degree of sample classification into high- and low-risk groups (Figure 2).

3.2.   Differential expression and functional enrichment of immune-microenvironment-related gene features in lung squamous cell carcinoma

To assess the relationship between whole gene expression profiles and immune scores, we compared high immune score subgroups with low immune score subgroups and screened for upregulated and downregulated differential genes. We defined the genes identified here as immune DEGs. The hierarchical clustering heat map shows the expression profiles of the top 15 upregulated and downregulated genes in the high- and low-immune score groups (Figure 3A). Among all DEGs (Figure 3A), 1020 genes showed significant upregulation, while 144 genes showed significant downregulation. We further performed follow-up analyses to elucidate the potential functions of immune-microenvironment-related DEGs (immune DEGs).

The relevance of certain functions was alluded to by the enrichment factor of GO and KEGG pathway analyses. Immune DEGs were particularly enriched in the following categories in the GO analysis: lymphocyte-mediated immunity, immune response-activating cell surface receptor signaling pathway), adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains, immunoglobulin complex, external side of plasma membrane, and antigen binding. Immune DEGs were significantly enriched in the hematopoietic cell lineage, cytokine-cytokine receptor interaction, Staphylococcus aureus infection, and cell adhesion molecules pathways in the KEGG pathway study (Figure 3B, C). GSEA analysis mainly showed significant differences in gene expression in cellular calcium ion homeostasis, receptor-ligand activity, regulation of GTPase activity, calcium signaling pathway, and cytokine-cytokine receptor interaction (Figure 3D, E).

3.3.   Identification and analysis of immune DEGs with prognostic value

We used a cross-set analysis of immune DEGs and survival significance genes to find more immune cell-related genes with prognostic value. Of the 1164 previously screened immune DEGs and the 3990 survival significance genes culled from the bulk survival data analysis (Figure 4A), 57 immune DEGs were found to be present in both and thus potentially be of prognostic value in patients with LUSC (Table 1). To further understand the mechanism of action of this gene set, we focused on the potential functions of the 57 intersection genes. Pathway analysis showed that these were mainly enriched in hematopoietic cell lineage, cytokine-cytokine receptor interaction and Staphylococcus aureus infection. In the GO report, the genes were mainly enriched in the acute inflammatory response, platelet degranulation, and positive regulation of protein secretion (Figure 4B, D, Table 2).

3.4.   Prognostic implication of IGLC7-associated immunomodulators in lung cancer

We generated a 9-gene LUSC prognostic profile. We then split the samples 1:1 to form a training set and a test set, and tested this risk model in the training set versus the overall data set (Figure 5A). The model was evaluated using Kaplan-Meier survival curves with cox p = 0.0017 in the training set and cox p < 0.0001 in the overall dataset (Figure 5B). According to the ROC curves, the ROC values were greater than 0.6 in both the training set as well as the overall dataset, indicating good predictive performance (Figure 5C).

To further explore the impact of this survival model on the prognosis of LUSC patients, we created scatter plots of gene expression and corresponding survival times in different samples. The death rate gradually increased as the risk score increased, indicating that the survival model could predict the survival of LUSC patients (Figure 6A). As shown in Figure 6, risk score was significantly associated with survival in LUSC in the univariate COX regression model [hazard ratio (HR) = 2.1, 95% confidence interval (CI) = 1.6-2.9, p < 0.001]. In addition, multivariate Cox regression showed that risk score remained an independent predictor of prognosis in LUSC after adjusting for gender, T, M, and N stage (HR = 2.22, 95% CI = 1.68-2.9, p < 0.001). In addition, multivariate Cox regression showed that IGLC7, also known as immunoglobulin lambda constant 7, was an independent predictor of prognosis in LUSC after adjusting for individual gene expression profiles (HR = 0.87, 95% CI = 0.81-0.94, p < 0.001) (Figure 6B, C).

Based on the C-index values, a nomogram integrating the risk score, gender, and TNM stage was constructed. Total points were calculated by adding the points of the risk score, age, and TNM stage. The calibration curves for predicting 3- and 5-year OS (Overall Survival) showed that the survival rates predicted by the nomograms closely correlated with the actual survival outcomes (Figure 6D, E). In addition, our analysis showed that the gene IGLC7 had the most significant effect on patient prognosis in this model. There are multiple pathways by which genes can influence patient prognosis, and immune-related genetic prognostic biomarkers are a potentially effective prognostic classifier. The effects of IGLC7 on immune cells and the synergistic effects of IGLC7 with other immune modulators are still largely unknown. To fill this gap, we explored the relevance of IGLC7 to immune cell populations, the synergistic effects of IGLC7 with immune checkpoints, and the relevance of IGLC7 to specific cellular immune and inflammatory responses.

3.5.   Relationship between IGLC7 expression and immune cell populations

In previous studies, only the physical location of human immunoglobulin lambda-like (IGLL) was demonstrated to be located on this gene, and its corresponding functional description is still unclear ^[12]. In patients with LUSC, the MCP-counter algorithm was used to calculate the absolute abundance of immune cell populations to clarify the immune manipulation function of IGLC7 further. Thus, as the expression of the IGLC7 gene increases, the expression abundance of immune cell populations likewise increases, indicating that IGLC7 expression plays a role in the immune response. Subsequently, we further calculated the correlation between IGLC7 and different immune cell populations. The marker genes for 24 immune cell species were used to predictthe abundance patterns of these cell populations in LUSC (Figure 7A). IGLC7 expression was strongly correlated with the immune cell population fractions of B cells (cor.p = 0.69), T cells (cor.p = 0.53), cytotoxic cells (cor.p = 0.48), regulatory T-cells (cor.p = 0.46), macrophages (cor.p = 0.44), T helper type 1cells (cor.p = 0.45) and T follicular helper cells (cor.p = 0.46), but less correlated with T helper 17cells, plasmacytoid dendritic cells, gamma delta T cells, and central memory T cells (Figure 7B). These results show that IGLC7 functions in T-cell and B-cell immune processes and is also involved in other cellular immune processes. This suggests that IGLC7 may possess multiple regulatory modes. Interestingly, other investigators have not revealed this result and may be a new insight for immunotherapy.

3.6.   Association between IGLC7 expression and specific cell immune responses

In LUSC, the precise immunomodulatory function of IGLC7 is unclear. In order to investigate the immunomodulatory mechanisms of IGLC7, after analyzing the correlation between IGLC7 and immune cells, we performed a GSVA analysis of its correlation with specific immune pathways. Consistent with the above results, we found that IGLC7 was closely and positively associated with multiple immune pathways (Figure 8A), particularly the activation pathway of B cells and the positive regulatory pathway of T cell proliferation. This suggests a detailed immunological function of IGCL7 and not just a physical location overlap. This result indicates the potential of IGLC7 in regulating the immune response and could be a target for future studies.

3.7.   Relationship between IGLC7 and inflammation

To provide further insight into the mechanisms by which IGLC7 mediates inflammatory activity, we used the GSVA algorithm to derive 104 genes from 7 clusters and defined them as metagenes associated with different inflammation and immune functions ^[13] (Table 3). These metagenes distinguished the prognosis of tumors in previous studies, and with this analysis, we identified a correlation between IGLC7 and the inflammatory response ^[14]. IGLC7 showed some degree of positive correlation with the STAT1, LCK, HCK, MHC-I, MHC-II, IgG, and interferon gene clusters (Figure 8B). Among these seven clusters, IGLC7 showed the strongest correlation with LCK. Previous studies have shown LCK to be an important mediator of B-cell receptor signaling in other diseases ^[15], in addition to studies demonstrating that LCK can act as a novel prognostic-related gene ^[16]. Thus, IGLC7 might regulate various immune cells, which will help in improving the effectiveness of immunotherapy for LUSC.

3.8.   Synergy of IGLC7 with other immune checkpoint members

We next evaluated the association of IGLC7 with other immune checkpoint members to further explore the synergistic role of IGLC7 in LUSC-induced immune responses, inflammatory pathways related to LUSC collected by Liu et al. ^[12,14] (Figure 8C). Detailed R and P values for the correlations between IGLC7 and other checkpoint members are listed in Table 4. We found that IGLC7 was closely associated with several checkpoint members, including CTLA4, BTLA, CD27, and CD48. Additionally, IGLC7 showed a positive correlation with all checkpoint members involved in the experiment, consistent with the correlation trend between IGLC7 and immune and inflammatory responses.

4.   Discussion

In this study, we developed a survival model that can accurately predict the survival of LUSC patients through TCGA database mining and experimental validation. We identified IGLC7 as having the most significant effect on overall survival in LUSC patients, potentially becoming a new tumor microenvironment regulator. We assessed the correlation of IGLC7 with immune cell populations, synergy with other immune checkpoints, and correlation with specific cellular immune and inflammatory responses.

First, we generated scores for 496 LUSC patients, calculated the score distribution based on clinical information, and divided the patients into high and low immune score groups according to their immune scores. We found a significant difference between the high and low immune score groups in the Kaplan-Meier curves between the high and low score groups (p = 0.039). According to the survival differences between the different groups, we demonstrated that the level of immunity, as implied by the immune scores that had an impact on patient prognosis. By analyzing differential gene expression between the groups, we identified 1164 genes associated with immune infiltration in LUSC.

The 1164 DEGs obtained were subjected to GO and KEGG pathway analysis as well as GSEA analysis. Most of the genes, such as MS4A1, ADH7, and GPX2, were involved in regulating the tumor microenvironment ^[17,18]. GO terminology consistently showed that these DEGs interact with immune cells in the immune microenvironment of lung squamous carcinoma, with enriched terms including lymphocyte-mediated immunity, immune response-activating cell surface receptor signaling pathway, and adaptive immune response based on somatic recombination of immune receptors built from immunoglobulin superfamily domains. Similarly, KEGG pathway analysis also suggests that the genes are related to an immune response, with pathways involving autoimmune thyroid disease and intestinal immune network for IgA production being enriched. Consistent with GO and KEGG analysis, GSEA analysis also showed that this gene set is involved in regulating immune processes.

Then, using batch overall survival analysis, we performed an intersection set analysis and identified a total of 57 survival gene markers (Table 1). The GO and KEGG analysis results confirmed that the 57 genes associated with the prognosis of the patients were largely involved in immunomodulatory functions.

Of the 57 genes, we identified a novel survival-related immunomodulator, IGLC7, with independent predictive power based on our constructed model. The immune microenvironment has been shown to play an important role in the field of tumour biology and numerous therapeutic approaches, including immune checkpoint inhibitors, immunotherapy and gene therapy, play a powerful role in the treatment of cancer. Given the critical role of immune-related gene prognostic biomarkers in disease diagnosis and immunotherapy, we performed a multifaceted immune function analysis of IGLC7. To further elucidate the immune function of IGLC7 in LUSC, the abundance of immune cell populations was calculated using a microenvironmental cell population counting algorithm. The results indicate that IGLC7 is potentially associated with various immunomodulators and immune cells, demonstrating that IGLC7 may be involved in complex immune regulation and can become a new immunomodulator.

The immune activation of both T cells and B cells is considered to evaluate immune-related functions ^[16,17]. We found that IGLC7 is closely associated with T cell immunity and B cell immunity, implying a significant role in immune regulation. It was noted in previous studies that Memory CD8⁺ T cells are critical in the immune process, but the development of vaccines targeting T cells remains problematic, mainly due to the limited knowledge of CD8⁺ T cells. In our study, we found a high correlation between IGLC7 and T cells, and by targeting IGLC7 treatment may fill the current gap and provide a new idea. There is still a need for more in-depth studies on how other immune cells are involved in the immune process of tumors ^[18], but IGLC7 is shown here to act on the regulation of other immune cells. This implies a multimodal regulatory role for IGLC7 in LUSC and that IGLC7 plays an important immune and inflammatory regulatory function in LUSC.

5.   Conclusions

In conclusion, IGLC7 expression levels are related to the malignancy of LUSC based on the constructed prognostic model and can thus be a therapeutic target for patients with LUSC. Furthermore, IGLC7 may work in concert with other immune checkpoint members to regulate the immune microenvironment of LUSC. These discoveries might lead to a fresh understanding of the complicated interactions between cancer cells and the tumor microenvironment, particularly the population of immune cells, and a novel approach to future immunotherapeutic treatments for patients with LUSC.

Acknowledgments

We would like to thank you for Jiangsu University medical clinical science technology development fund (JLY2021061) to support our work.

Conflict of interest

The authors declare there is no conflict of interest.