Skip to main content


Front. Oncol., 20 October 2021
Sec. Thoracic Oncology

Identification of Prognostic Model Based on Immune-Related LncRNAs in Stage I-III Non-Small Cell Lung Cancer

Qiaxuan Li1,2†, Lintong Yao1,2†, Zenan Lin3†, Fasheng Li1,4, Daipeng Xie1, Congsen Li1,2, Weijie Zhan1,4, Weihuan Lin4, Luyu Huang1,2, Shaowei Wu1,2 and Haiyu Zhou1,2,4,5*
  • 1Department of Thoracic Surgery, Guangdong Provincial People’s Hospital & Guangdong Academy of Medical Sciences, School of Medicine, South China University of Technology, Southern Medical University, Guangzhou, China
  • 2College of Medicine, Shantou University, Shantou, China
  • 3Guangdong Provincial People’s Hospital, Guangdong Academy of Medical Sciences, Guangzhou University of Chinese Medicine, Guangzhou, China
  • 4The Second School of Clinical Medicine, Southern Medical University, Guangzhou, China
  • 5Department of Thoracic Surgery, Jiangxi Cancer Hospital, Nanchang, China

Background: Long non-coding RNAs (lncRNAs) participate in the regulation of immune response and carcinogenesis, shaping tumor immune microenvironment, which could be utilized in the construction of prognostic signatures for non-small cell lung cancer (NSCLC) as supplements.

Methods: Data of patients with stage I-III NSCLC was downloaded from online databases. The least absolute shrinkage and selection operator was used to construct a lncRNA-based prognostic model. Differences in tumor immune microenvironments and pathways were explored for high-risk and low-risk groups, stratified by the model. We explored the potential association between the model and immunotherapy by the tumor immune dysfunction and exclusion algorithm.

Results: Our study extracted 15 immune-related lncRNAs to construct a prognostic model. Survival analysis suggested better survival probability in low-risk group in training and validation cohorts. The combination of tumor, node, and metastasis staging systems with immune-related lncRNA signatures presented higher prognostic efficacy than tumor, node, and metastasis staging systems. Single sample gene set enrichment analysis showed higher infiltration abundance in the low-risk group, including B cells (p<0.001), activated CD8+ T cells (p<0.01), CD4+ T cells (p<0.001), activated dendritic cells (p<0.01), and CD56+ Natural killer cells (p<0.01). Low-risk patients had significantly higher immune scores and estimated scores from the ESTIMATE algorithm. The predicted proportion of responders to immunotherapy was higher in the low-risk group. Critical pathways in the model were enriched in immune response and cytoskeleton.

Conclusions: Our immune-related lncRNA model could describe the immune contexture of tumor microenvironments and facilitate clinical therapeutic strategies by improving the prognostic efficacy of traditional tumor staging systems.


Lung cancer has the leading rates of incidence and mortality worldwide, with the highest estimated deaths and 119,100 newly diagnosed cases in USA from Cancer statistic, 2021 (1). Non-small cell lung cancer (NSCLC) is the most common subtype, accounting for 85% of lung cancer cases. Target therapy and immune checkpoint inhibitors are emerging therapeutic strategies for NSCLC, but the identification of potential responders remains critical (2).

Long non-coding RNAs (lncRNAs) are diverse repertoires of RNA transcripts that are over 200 nucleotides in length, which lack the capacity for direct protein coding but are involved in the regulation of critical biological processes and cellular behavior (3). LncRNAs influence gene expression and have essential roles in carcinogenesis by combining regulatory molecules with proteins or directly binding to nucleic acids, protein complexes, or transcription factors (4). The roles of lncRNAs in immune systems are extensively investigated, and it has been summarized that lncRNAs are involved in the differentiation, activation, and function of immune cells. Emerging research has emphasized the role of lncRNAs in the regulation of carcinogenesis, immunosurveillance, and antitumor immune responses. (5) Specific lncRNAs have been found overexpressed in tumor-associated macrophages, shaping the tumor immune microenvironment and inhibiting tumor apoptosis (6).

Tumor infiltrating immune cells, compromising dendritic cells, mast cells, natural killer (NK) cells, macrophages, and tumor infiltrating lymphocytes (TILs), are present in the tumor microenvironment such as the tumor center peritumor stroma, or invasive margin (7). Immune checkpoints such as LAG-3, CTLA-4, and PD-L1 expressed in TILs are associated with cancer prognosis and therapeutic response, especially immunotherapy (8). It is evident that CD8+ T cells in tumor microenvironments are critical to the immune response. Recent studies have emphasized the function of tumor infiltrating B cells in immune surveillance and regulation of immunotherapy (9).

The traditional staging system for non-small cell lung cancer is the tumor, node, and metastasis (TNM) classification, which stratifies tumor subtypes and predicts cancer prognosis based on tumor size and its invasiveness to lymph nodes or distant organs, but does not take account of the tumor molecular and immune characteristics (5). Immune characteristics could be good candidates to improve prognostication of the TNM staging system in NSCLC (10).

Based on the fact lncRNAs could shape tumor microenvironments and predict the characteristics of NSCLC, we focused on immune-related lncRNAs to filter effective prognostic signatures. In this study, we extracted significant immune-related lncRNAs in stage I-III NSCLC and constructed immune-related lncRNAs based prognostic model, exploring the immune characteristics in tumor microenvironments and the potential therapeutic response of immunotherapy.

Material and Methods

Data Acquisition

All the clinical information and RNA-sequencing data of patients with NSCLC were downloaded from the Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO). NSCLC patients with clinical stage I-III samples and complete follow-up information were included in this study. After sifting, 1357 cases with NSCLC from three data sets were incorporated into this study, including 970 patients from TCGA (, 226 patients in GSE31210 (, and 161 patients in GSE30219 ( The 970 patients from TCGA were randomly separated in a 7:3 ratio to form the training cohort (n=717) and the testing cohort (n=253). GSE31210 and GSE30219 were combined into another independent validation cohort. Detailed baseline clinical features of three datasets are shown in Table 1. Batch effects were removed by the “ComBat” package of R software. Meanwhile, immune-related lncRNAs were downloaded and extracted from the Immlnc dataset ( (11).


Table 1 Characteristics baseline of patients in cohorts.

Construction and Validation of the Immune-Related LncRNA Model

The immune-related lncRNA model was identified using the training cohort, and the validation cohort and GEO datasets were used to validate the accuracy and efficacy of the model. We selected the immune-related lncRNAs by taking the intersection of lncRNAs between and TCGA and GEO datasets. The least absolute shrinkage and selection operator (LASSO) was chosen to reduce overfitting and to analyze the optimal immune-related lncRNA signature for predicting the overall survival of NSCLC patients. The “glmnet” R package was used for LASSO regression analysis. We calculated the risk score of each sample as follows: risk score = expression value of lncRNA 1 * coefficient + expression value of lncRNA 2 * coefficient + … + expression value of lncRNA n * coefficient. Then, the NSCLC patients were sorted into high-risk and low-risk groups based on the optimal cut-off value of the risk score. The area under the curve (AUC), performed by the “Survival” package of R software, was used to validate the sensitivity and specificity of the immune-related signature.

Tumor Infiltrating Immune Cells Signature

Single sample gene set enrichment analysis (ssGSEA) was applied to quantify the immune infiltration level of 28 immune cell phenotypes. We obtained the gene set from previous studies, which included various immune cell phenotypes such as activated B cells, activated CD8 T cells, T follicular helper cell, and so on (12, 13). Also, CIBERSORT ( was applied to describe the abundance of 22 immune cell types in each NSCLC sample using R software. We evaluated the immune cell infiltration, stromal content, and tumor purity with the ESTIMATE algorithm. By comparing the results of ssGSEA, CIBERSORT, and ESTIMATE among high-risk and low-risk groups, we described the relationship between the tumor immune microenvironment and the immune-related lncRNA signature.

Prediction of Potential Immunotherapy Response

A T cell dysfunction and exclusion signature was constructed by the tumor immune dysfunction and exclusion (TIDE, algorithm to predict tumor immune escape, which could influence patients’ response to immunotherapy. We used TIDE score mapping to compare the potential clinical response to immune checkpoint inhibitors between high-risk and low-risk groups.

Potential Biological Mechanisms of the Immune-Related LncRNA Model

We downloaded the single nucleotide variation (SNV) data from TCGA to analyze the difference of SNV among high-risk and low-risk groups through Fisher test. To explore the potential biological mechanisms, we used the “Limma” package of R software to identify the differentially expressed genes (DEGs) between high-risk and low-risk groups in training cohort. The gene ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed to analyze the functional enrichment of DEGs, using R software. The protein-protein interaction (PPI) network of DEGs was described by using the STRING database (Version 11.0) and constructed by using Cytoscape (Version 3.8.2). Then, gene set enrichment analysis (GSEA) was performed on GSEA software ( by using the Molecular Signatures Database (MSigDB) Version 7.2 collections C2 (curated gene sets) (14).

Verification of the LncRNA Expression Between NSCLC Tissues and Adjacent Normal Tissues by qRT-PCR

We collected sixteen paired NSCLC and adjacent normal tissue samples from Jiangxi Cancer Hospital after gaining ethical approval from the Human Research Ethics Committee in Jiangxi Cancer Hospital. Total RNA was isolated using RNAiso Plus (Takara) according to the manufacturer’s instructions. The RNA was reverse-transcribed using the Primer Script Reverse Transcriptase Reagent Kit with gDNA Eraser (Takara, RR047A). Real-time PCR was performed using the TB Green™ Premix Ex Taq™ (Takara, RR420A) and analyzed using the Bio-Rad CFX96 thermal cycler. The primer sequences used for the investigated genes are listed in online Supplementary Table 3. GADPH was used to standardize the gene expression. In order to compare the expression levels of lncRNA in different samples, the 2-ΔΔCt method was adopted to calculate the expression levels of the immune-related lncRNA from the risk model.

Statistical Analysis

Univariate cox proportional hazard regression was applied to assess the prognostic value of immune-lncRNA signatures by evaluating the association between risk score and overall survival in the training cohort. The correlation of the overall survival with immune-related lncRNA signature and the clinicopathological characteristics was calculated using the Kaplan-Meier curve. The Wilcoxon rank sum test and t-test were conducted for the comparison between two groups. Two-tailed P value < 0.05 was considered significant statistically. All statistical analyses were performed in R software, version 4.0.1.


Construction and Verification of the Immune-Related LncRNA Model

We extracted 1034 immune-related lncRNAs by taking the intersection among the lncRNAs in the TCGA-LUAD dataset, GEO and Immlnc database. We filtrated 15 immune-related lncRNAs with LASSO, then constructed the prognostic immune-related lncRNA model. (Figures 1A, B) Multivariate cox regression analysis showed potential prognostic properties in these 15 immune-related lncRNAs, and the expression level of 9 of the lncRNAs were positively associated with overall survival. Reciprocally, 6 immune-related lncRNAs were correlated with worse prognosis. All 15 lncRNAs are either significantly protective or risk factors for survival of stage I-III LUAD. (Figure 1C) Survival analysis of these 15 lncRNAs suggest significant survival differences between high and low expression levels of lncRNAs. (Supplementary Figure 1) Equation for the risk model from 15 significant immune-related lncRNAs is exhibited as follows:


Figure 1 Construction and verification of the immune-related lncRNA prognostic model. (A, B) The least absolute shrinkage and selection operator was utilized to construct an immune-related lncRNA model. (C) The forest plot of fifteen immune-related lncRNAs was figured by multivariate cox regression analysis. (D–F) Survival analysis showed better survival among low-risk patients in training cohort, internal validation cohort, and external validation cohort. lncRNA, long non-coding RNA.

Risk score = 0.0136*expression value of LINC01116 + 0.1285*expression value of WWC2.AS2 + 0.0415*expression value of CASC15 + 0.0545*expression value of PRKG1.AS1 + 0.1052*expression value of HOTAIR + 0.0011*expression value of TMPO.AS1 - 0.0711*expression value of CDC42.IT1 - 0.1914*expression value of BANCR - 0.0866*expression value of COLCA1 - 0.1751*expression value of RPARP.AS1 - 0.1610*expression value of CRNDE - 0.0422*expression value of C20orf197 - 0.1677*expression value of CSNK1G2.AS1 - 0.1954*expression value of LINC00528 - 0.0693*expression value of LINC00896.

Training cohort was divided into a high-risk group (n=109) and a low-risk group (n=608) according to the optimal cut-off value which was most significantly associated with overall survival. The high-risk group showed higher mortality than the low-risk group. The Kaplan-Meier survival analysis of the immune-related lncRNA signatures suggested better a survival probability in the low-risk group of the training cohort. (Figure 1D) Similar results for score distribution and survival analysis were found in the validation cohort and GEO cohort (Figures 1E, F).

Clinical Value of Immune-Related LncRNA Model

Univariate cox regression analysis of the clinical characteristics and immune-related lncRNA model suggested that gender, AJCC T stage, AJCC N stage, AJCC TNM stage, and immune lncRNA model were significant prognostic factors for overall survival in NSCLC. All the significant candidates were included into multivariate cox regression analysis, indicating that AJCC T stage and immune-related lncRNA model were independent prognostic factors (Table 2).


Table 2 Univariate and multivariate analysis of clinical characteristics.

To explore the robustness of the prognostic effect of the immune-related lncRNA model, Kaplan-Meier survival analysis was performed and stratified by clinicopathological characteristics, including T stage (T1, T2, T3-T4), N stage (N0, ≥N1), and gender (male, female). Similar results were found in the high-risk group, which had worse overall survival than low-risk group, and were delivered from different gender, T stage, and N stage in the training cohort (Figures 2A–G).


Figure 2 Clinical value of immune-related lncRNA model. (A–G) Survival analysis showed favorable survival for low-risk patients in different gender, node stage, and tumor stage. (H, I) Receiver operating characteristic curves were used to compare the predictive efficacy of the immune-related lncRNA based model alone, AJCC TNM staging alone, and the combination model in training cohort and validation cohort. lncRNA, Long non-coding RNA; AJCC, American Joint Committee on Cancer; TNM staging, Tumor, Node, and metastasis staging system.

The area under the curve (AUC) of the combination between the AJCC TNM staging system and the immune-related lncRNA model was 0.90, while the AUC of AJCC TNM staging system alone was 0.88 in the training cohort for 3 years survival (Figure 2H). The combination of the AJCC TNM staging system and the immune-related lncRNA model presented a larger AUC than the AJCC TNM staging system alone in the internal validation cohort (AUC: 0.92 vs 0.89, Figure 2I) and external validation cohort (AUC: 0.64 vs 0.61, Supplementary Figure 2). Enhancement of prognostic accuracy of AJCC staging indicated a potential clinical application of the immune-related lncRNA model.

Exploration of Immune Landscape and Immune Response

Tumor infiltrating lymphocytes analysis was conducted through ssGSEA and CIBERSORT. Through ssGSEA, almost all of the infiltrating immune cells showed higher infiltration abundance in low-risk group, especially activated B cells (p<0.001), immature B cells (p<0.001), effector memory CD4 T cells (p<0.001), activated CD8 T cells (p<0.01), effector memory CD8 T cells (p<0.001), activated dendritic cells (p<0.01), immature dendritic cells (p<0.001), CD56+ Natural killer cells (p<0.01), mast cells (p<0.001), monocytes (p<0.01), and T follicular helper cells (p<0.001) (Figure 3A). A higher proportion of CD8+ tumor infiltrating lymphocytes could be detected in the low-score group with CIBERSORT (Supplementary Figure 3). Tumor purity as calculated by the ESTIMATE algorithm was not a significant prognostic factor for NSCLC based on univariate cox regression analysis in our study (Table 2).


Figure 3 Exploration of immune landscape and immune response. (A) Single sample gene set enrichment analysis suggested a higher proportion of multiple immune cells such as activated B cells, CD8+ T cells, CD4+ T cells, and dendritic cells. ns for p>0.001, * for p<0.001, ** for p<0.0001, *** for p<0.00001. (B–D) Higher immune score, estimate score, and lower tumor purity are analyzed by ESTIMATE algorithm. (E, F) TIDE analysis estimated T cell dysfunction and exclusion and predicted response of immunotherapy. (G) Low-risk patients showed significantly higher PD-1 and CTLA-4. lncRNA, Long non-coding RNA; ssGSEA, Single Sample Gene Set Enrichment Analysis; TIDE, Tumor Immune Dysfunction and Exclusion. ns for p>0.05, * for p<0.05, ** for p<0.01, *** for p<0.001.

To assess the immune infiltration in histological aspects, the ESTIMATE algorithm was used to infer the proportion of stromal cells and immune cells in tumor tissue. The low-risk group has higher immune score (p<0.001) and estimate score (p=0.004) significantly, which indicates more immune cell infiltration and lower tumor purity (p=0.004). (Figures 3B–D) ESTIMATE analysis verified the results of the tumor infiltrating lymphocyte evaluation with ssGSEA and CIBERSORT.

Observing the correlation between the immune-related lncRNA model and the tumor microenvironment, we further explored the association between this model and the response of immunotherapy. The TIDE score integrates the T cell dysfunction and exclusion signature to evaluate tumor immune escape. The high-risk group calculated from immune-related lncRNA signature was found to have a higher TIDE score, indicating potential tumor T cell dysfunction and exclusion. The predicted proportion of responders to immune checkpoint blockade was lower in the high-risk group as well. (Figures 3E, F) In addition, among patients in low-risk group, the gene expression of PDCD1 and CTLA-4 were significantly higher than those in high-risk group (Figure 3G).

Somatic Mutation and Pathway Analysis of the Immune-Related LncRNA Model

In order to compare the difference of somatic mutation of the immune-related lncRNA model, we used the “maftools” R package to calculate the SNV among the high-risk and low-risk groups. We found that the high-risk group had high mutation of TP53, TTN, MUC16, RYR2, CSMD3, XIRP3, USH2A, ZFHX4, LRP1B, and KEAP1 (Supplementary Figure 4A). The low-risk group was characterized by frequent mutation of TP53, TTN, MUC16, CSMD3, RYR2, LRP1B, USH2A, ZFHX4, KRAS, and FLG (Supplementary Figure 4B). Meanwhile, we analyzed the level of EGFR mutation in the two groups and we found no difference between the high-risk group and low-risk group (Fisher test, p=0.246, Supplementary Figure 4C).

We identified 19592 differentially expressed genes (DEGs) between the high-risk and low-risk groups, among which 72 DEGs were significant (p<0.05) and had fold change ≥2. GO functional enrichment analysis showed up-regulated DEGs for the low-risk group enriched in humoral immune response and channel activity, and down-regulated DEGs for the low-risk group enriched in skin development, keratinocyte differentiation, and cytoskeleton. (Figures 4A, B) KEGG pathway analysis indicated significant enrichment in the PI3K-AKT signaling pathway, Ras/Rap signaling pathway, MAPK signaling pathway, and regulation of actin cytoskeleton. (Supplementary Figure 5) GSEA showed significantly positive enrichment of the Fc epsilon RI pathway, asthma, vascular smooth muscle contraction, T cell receptor signaling, B cell receptor signaling, renin-angiotensin-aldosterone system, and VEGF signaling pathways in the low-risk group. (Figures 4C, D and Supplementary Table 1) Meanwhile, up-regulated pathways enriched in high-risk group were DNA replication, cell cycle, nucleotide excision repair, p53 signaling, homologous recombination, mismatch repair, pentose phosphate pathway, and tricarboxylic acid cycle. (Figures 4E, F; Supplementary Figure 6 and Supplementary Table 2).


Figure 4 Pathway analysis of immune-related lncRNA model. (A, B) Gene ontology analysis was used to explore the potential functional mechanism of immune-related lncRNA model, and the results are visualized in the low-risk group (A) and high-risk group (B). Immune-related lncRNA signaling pathway obtained by gene set enrichment analysis, including T cell receptor signaling pathway (C), Fc epsilon RI pathway (D), p53 signaling pathway (E), and cell cycle (F).

Expression Level of 15 Immune-Related LncRNAs Between NSCLC Tissues and Adjacent Normal Tissues by qRT-PCR

Finally, we measured the expression levels of 15 immune-related lncRNA from the risk model in the sixteen paired NSCLC and adjacent normal tissues using qRT-PCR. The result showed that three lncRNAs (CASC15, BANCR, and RPARP.AS1) had low expression in NSCLC tissues when compared with normal tissues (Supplementary Figure 7).


We constructed a 15 immune-related lncRNA based model using the LASSO algorithm for NSCLC in the TCGA dataset and validated it using the GEO datasets. Combined model demonstrated an enhanced prognostic efficacy compared to AJCC TNM staging system alone. A higher proportion of immune cell infiltration was detected in the survival-benefit group stratified by our model. The immune infiltration score described by the ESTIMATE algorithm is higher in the survival-benefit group, which is predicted with better immunotherapy response reasonably. GO functional analysis showed an enrichment in the humoral immune response, while KEGG showed enrichments in the PI3K-AKT signaling pathway, MAPK signaling pathway, and Ras/Rap signaling pathway. GSEA showed enrichments in the Fc epsilon RI pathway, asthma, T cell receptor signaling, and B cell receptor signaling.

Our study constructed and validated an immune-related lncRNA signature for stage I-III NSCLC with good prognostic AUC in all enrolled datasets. A total of 15 lncRNAs were identified to construct the immune-related lncRNA model, including 9 lncRNAs (CDC42.IT1, BANCR, COLCA1, RPARP.AS1, CRNDE, C20orf197, CSNK1G2.AS1, LINC00528, and LINC00896) associated with better overall survival and 6 lncRNAs (LINC01116, WWC2.AS2, CASC15, PRKG1.AS1, HOTAIR, and TMPO.AS1) correlated with worse prognosis. Some of these lncRNAs had been confirmed to be related to cancer progression and prognosis in previous studies. CASC15 has been reported to be upregulated in various types of tumor tissues (15), including NSCLC. As part of HIF-1α/CASC15/SOX4/β-catenin axis, CASC15 plays an essential role in cell proliferation, invasion, and tumor development in NSCLC (16). Meanwhile, CASC15 can promote lung cancer metastasis via miR-766-5p/KLK12 axis (17). Homebox (HOX) are vital in embryonic development and oncogenesis and the most studied HOX-lncRNAs is HOTAIR. HOTAIR is also significantly upregulated in NSCLC and is known for its association with higher TMN staging, lymphatic metastasis, and poor prognosis (18). The previous studies had found that HOTAIR could promote the level of miR-149-5p to facilitate the process of invasion, migration, and cell proliferation in NSCLC (19). Furthermore, HOTAIR has been linked to drug resistance in several types of tumor. Silencing HOTAIR expression can revert the gefitinib resistance of lung adenocarcinoma (20). Several studies have emphasized the potential value of LINC01116 as a prognostic marker or therapeutic target in various kinds of cancer. LINC01116 has been confirmed to accelerate tumor progression by regulating tumor-associated genes such as MYC (21) and p53 (22). In lung cancer, the upregulation of LINC01116 is an important reason for tumor cell proliferation and migration as it enhances the process of epithelial-mesenchymal transition (EMT) (23). In addition, as one of effectors in the IFN/STAT1 pathway, IFI44 is repressed by LINC01116, leading to acquired resistance of gefitinib in NSCLC (24). TMPO.AS1 performs its tumor-promoting function via activating the PI3K/AKT/mTOR pathway in gastric cancer and HCC (25). However, studies on BANCR and CRNDE function on tumor procession have shown conflicting results. For BANCR, some researchers have identified its influence in tumor invasion and migration (26). In contrast, overexpression of BANCR can control the content of N-cadherin, and vimentin and E-cadherin were shown to inhibit EMT in another study. (27) In our study, BANCR was used as a positive prognostic predictor for overall survival in NSCLC. As for CRNDE, most studies show that CRNDE can promote cell proliferation, invasion, and migration and inhibit apoptosis in colorectal cancer, lung cancer, glioma, and other cancers. CRNDE was also shown to affect cancer microenvironments and metabolism via the PI3K/AKT/mTOR and Raf/MAPK pathways. In conclusion, CRNDE might be a potential cancer promoter (28). However, recent research indicates the unique function of CRNDE in attenuating chemoresistance in gastric cancer by reducing the stability of SRSF6 (29). Meanwhile, our immune-related lncRNA signature also describes CRNDE as a favorable factor for overall survival in NSCLC. KEGG pathway analysis of our signature also shows enrichment of the PI3K/AKT and MAPK signaling pathways. The molecular mechanism of CRNED in cancer prognosis has not been completely investigated, and requires additional experiments to explore it in the future.

The TNM staging system is a cancer staging manual constructed by the American Joint Committee on Cancer on anatomic extent, in the pursuit of definitive prognoses and selecting the most beneficial therapeutic strategies. However, clinical outcomes vary among different patients in the same stage due to diverse biological behavior determined by molecular and genetic features (30). Immune contexture represents the results of dynamic interaction between tumor cells and the immune system in the tumor microenvironment. Prognostic immune parameters have been studied to predict the prognosis of cancer. Immunoscore is an unprecedented biomarker describing the proportion of immune cells in the tumor centre, invasive edge, and peritumor stroma (31). In colorectal cancer, immunoscore is used to predict prognosis, therapeutic effects, and disease relapse after immune checkpoint inhibitor therapy (32). With respect to NSCLC, a Norwegian study identified stromal CD8+ density and CD45RO+ memory T lymphocytes as independent prognostic factors for NSCLC regardless of endpoints, and proposed them as a supplement to the TNM-staging system (10). Prospective multicenter clinical trials are designed to evoke the attention of TNM-immunoscore for clinical application (33). After identifying prognostic immune-related lncRNA signatures, the differences of immune contexture including tumor infiltrating lymphocytes and immunescore were explored with our immune-related lncRNA model. An improvement of prognostic efficacy was found in the training cohort, internal validation cohort, and external validation cohort, and provides a novel immune indicator as a nonanatomic supplement of tumor features for the TNM staging system.

Infiltration of different types of immune cells is associated with cancer progression and patient survival in NSCLC (34). Our study identified differences in tumor infiltrating cells between high-risk and low-risk groups according to an immune-related lncRNA model, focusing on cells including dendritic cells, B cells, CD4 T cells, CD8 T cells, Natural killer cells, T follicular helper cells, and mast cells. Cytotoxic CD8+ T lymphocytes (CTLs) are essential immune cells against tumor cells (35). The priming of CTLs requires antigen presentation and co-interaction with mature dendritic cells, natural killer cells, and CD4+ T cells. Mature dendritic cells and natural killer cells, motivated and licensed by CD4+ T cells, secret costimulatory molecules and cytokines, and then CTL is priming (36). Infiltrating B cells are emphasized as active participants that orchestrate the antitumor immune response. B cells are involved in antigen presentation to T cells directly, or facilitate the antigen uptake of dendritic cells (37). T follicular helper cells rely on antigen-specific B cells, and reciprocally facilitate B cell proliferation and differentiation, which is crucial in humoral response (38). Significant associations have been verified by multiple studies between the clinical outcome of cancer patients and CD4+ T cells, CD8+ T cells, B cells, dendritic cells, natural killer cells, and mast cells (39, 40).

In the present study, we inferred that atopy may be relevant to the development and prognosis of lung cancer based on evidence that the low-risk group had a higher infiltration of mast cells, and enriched Fc epsilon RI pathways, asthma and vascular smooth muscle contractions from GSEA. An up-to-30 year prospective study with 37747 participants in Denmark observed a 10-fold higher IgE level in non-Hodgkin lymphoma, oral or pharynx cancer, lung cancer, and esophagus cancer. However, these results were non-significant after multivariable adjustment (41). Recent research based on the Surveillance, Epidemiology, and End Results (SEER) database in the United States suggested asthma was associated with reduced risk of liver cancer, which could be attributed to the activation of immunosurveillance from allergic response (42). The cellular mechanism of the antitumor function of IgE/FcϵRI combination could be explained by cross-presentation with dendritic cells to induce the priming of cytotoxic T lymphocytes (43). High-affinity IgE receptors and Fc epsilon RI signaling pathways constituted the most significant pathway of prognostic signature for lung adenocarcinoma. MS4A2, an IgE receptor related gene expressed in tumor-infiltrating mast cells, was an independent favorable prognostic biomarker (44). Ultra-low IgE is correlated with a higher risk of malignancy and could be a diagnostic and prognostic biomarker for lung adenocarcinoma (42). Hence, our results validate the critical role of IgE and allergic reactions in the antitumor response of low-risk patients, stratified based on immune-related lncRNA signature.

Immune checkpoint inhibitor (ICI) immunotherapy has emerged as an effective treatment for various kinds of cancers, including NSCLC (45). However, only 20%-30% patients with NSCLC respond to immunotherapy. Predictive biomarkers, such as PD-1 (programmed cell death protein 1) and CTLA-4 (cytotoxic T lymphocyte-associated antigen-4), are frequently used to assess the response of ICI in NSCLC. In our study, immune checkpoint genes such as PD-1 and CTLA-4 were significantly higher in the low-risk group than in high-risk group (46). The TIDE scores were also significantly higher in the low-risk group. Thus, patients in the low-risk group may have better response to ICI therapy. Tumor-infiltrating immune cells have also been regarded as a predictor for response to immunotherapy (46). For example, KEYNOTE-001 has found that a high percentage of CD8+ T cell infiltration showed a strong association with superior ICI treatment responses when treated with pembrolizumab. (47) We also found that the low-risk group showed higher immune cells infiltration abundance, such as activated B cells, immature B cells, CD8+ T cells, and so on. In conclusion, our 15-immune-related lncRNA signature was closely related to ICI response.

Some limitations of this study should be considered when interpreting the results. First, all of the results were completed based on retrospective studies and public datasets. The accuracy of our immune-related lncRNA signature should be further verified with a clinical real-world dataset. Second, laboratory explorations are needed to verify and illuminate the molecular mechanisms of these immune-related lncRNAs. Finally, as the predictive marker of ICI response, TIDE score is just verified in several datasets. Hence, available immunotherapy cohorts are warranted to confirm the clinical application of the immune-related lncRNA signature.

In conclusion, we established a 15-immune-related lncRNA prognostic model for NSCLC. This model could be applied in clinical situations as a supplement to the TNM staging system for an improvement in the predictive efficacy of cancer malignancy and prognosis. Low-risk patients stratified by our model have higher infiltration of immune cells such as dendritic cells, CD8+ T cells, CD4+ T cells, B cells, natural killer cells, and mast cells. Pathway analysis of our immune-related lncRNA signature might indicate an underlying mechanism associated with humoral immunity, cell-mediated immunity, and the regulation of the cell cycle.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Author Contributions

QL, LY, and HZ designed this work. ZL performed experiments such as qRT-PCR. FL and DX integrated the data and conducted the analyses. CL and WZ wrote this manuscript. WL, LH, and SW edited and revised themanuscript. All authors contributed to the article and approved the submitted version.


Guangzhou Municipal Science and Technology Project and Natural Science Foundation of Guangdong Province.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.


We appreciate the data obtained from the TCGA database and the GEO database. This study was funded by the Science and Technology Program of Guangzhou (Grant Number 201903010028) and Natural Science Foundation of Guangdong Province, China (Grant Number 2021A1515010838).

Supplementary Material

The Supplementary Material for this article can be found online at:


1. Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer Statistics, 2021. CA Cancer J Clin (2021) 71(1):7–33. doi: 10.3322/caac.21654

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Duma N, Santana-Davila R, Molina JR. Non-Small Cell Lung Cancer: Epidemiology, Screening, Diagnosis, and Treatment. Mayo Clin Proc (2019) 94(8):1623–40. doi: 10.1016/j.mayocp.2019.01.013

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Quinn JJ, Chang HY. Unique Features of Long Non-Coding RNA Biogenesis and Function. Nat Rev Genet (2016) 17(1):47–62. doi: 10.1038/nrg.2015.10

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Prensner JR, Chinnaiyan AM. The Emergence of Lncrnas in Cancer Biology. Cancer Discov (2011) 1(5):391–407. doi: 10.1158/2159-8290.CD-11-0209

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gibb EA, Brown CJ, Lam WL. The Functional Role of Long Non-Coding RNA in Human Carcinomas. Mol Cancer (2011) 10:38. doi: 10.1186/1476-4598-10-38

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Slack FJ, Chinnaiyan AM. The Role of Non-Coding Rnas in Oncology. Cell (2019) 179(5):1033–55. doi: 10.1016/j.cell.2019.10.017

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Bremnes RM, Busund LT, Kilvaer TL, Andersen S, Richardsen E, Paulsen EE, et al. The Role of Tumor-Infiltrating Lymphocytes in Development, Progression, and Prognosis of Non-Small Cell Lung Cancer. J Thorac Oncol Off Publ Int Assoc Study Lung Cancer (2016) 11(6):789–800. doi: 10.1016/j.jtho.2016.01.015

CrossRef Full Text | Google Scholar

8. He Y, Yu H, Rozeboom L, Rivard CJ, Ellison K, Dziadziuszko R, et al. Lag-3 Protein Expression in Non-Small Cell Lung Cancer and Its Relationship With Pd-1/Pd-L1 and Tumor-Infiltrating Lymphocytes. J Thorac Oncol Off Publ Int Assoc Study Lung Cancer (2017) 12(5):814–23. doi: 10.1016/j.jtho.2017.01.019

CrossRef Full Text | Google Scholar

9. Paijens ST, Vledder A, de Bruyn M, Nijman HW. Tumor-Infiltrating Lymphocytes in the Immunotherapy Era. Cell Mol Immunol (2020) 18(4):842–59. doi: 10.1038/s41423-020-00565-9

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Donnem T, Hald SM, Paulsen EE, Richardsen E, Al-Saad S, Kilvaer TK, et al. Stromal Cd8+ T-Cell Density-A Promising Supplement to Tnm Staging in Non-Small Cell Lung Cancer. Clin Cancer Res (2015) 21(11):2635–43. doi: 10.1158/1078-0432.CCR-14-1905

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Li Y, Jiang T, Zhou W, Li J, Li X, Wang Q, et al. Pan-Cancer Characterization of Immune-Related Lncrnas Identifies Potential Oncogenic Biomarkers. Nat Commun (2020) 11(1):1000. doi: 10.1038/s41467-020-14802-2

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Jia Q, Wu W, Wang Y, Alexander PB, Sun C, Gong Z, et al. Local Mutational Diversity Drives Intratumoral Immune Heterogeneity in Non-Small Cell Lung Cancer. Nat Commun (2018) 9(1):5361. doi: 10.1038/s41467-018-07767-w

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Xue Y, Tong L, LiuAnwei Liu F, Liu A, Zeng S, Xiong Q, et al. Tumorinfiltrating M2 Macrophages Driven by Specific Genomic Alterations Are Associated With Prognosis in Bladder Cancer. Oncol Rep (2019) 42(2):581–94. doi: 10.3892/or.2019.7196

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (Msigdb) Hallmark Gene Set Collection. Cell Syst (2015) 1(6):417–25. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Gu X, Chu Q, Zheng Q, Wang J, Zhu H. The Dual Functions of the Long Noncoding RNA Casc15 in Malignancy. Biomed Pharmacother = Biomed Pharmacother (2021) 135:111212. doi: 10.1016/j.biopha.2020.111212

CrossRef Full Text | Google Scholar

16. Sun J, Xiong Y, Jiang K, Xin B, Jiang T, Wei R, et al. Hypoxia-Sensitive Long Noncoding RNA Casc15 Promotes Lung Tumorigenesis by Regulating the Sox4/Beta-Catenin Axis. J Exp Clin Cancer Res (2021) 40(1):12. doi: 10.1186/s13046-020-01806-5

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Bai Y, Zhang G, Cheng R, Yang R, Chu H. Casc15 Contributes to Proliferation and Invasion Through Regulating Mir-766-5p/ Klk12 Axis in Lung Cancer. Cell Cycle (2019) 18(18):2323–31. doi: 10.1080/15384101.2019.1646562

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Li L, Wang Y, Song G, Zhang X, Gao S, Liu H. Hox Cluster-Embedded Antisense Long Non-Coding Rnas in Lung Cancer. Cancer Lett (2019) 450:14–21. doi: 10.1016/j.canlet.2019.02.036

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Li H, Cui Z, Lv X, Li J, Gao M, Yang Z, et al. Long Non-Coding RNA Hotair Function as a Competing Endogenous RNA for Mir-149-5p to Promote the Cell Growth, Migration, and Invasion in Non-Small Cell Lung Cancer. Front Oncol (2020) 10:528520. doi: 10.3389/fonc.2020.528520

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Liu Y, Jiang H, Zhou H, Ying X, Wang Z, Yang Y, et al. Lentivirus-Mediated Silencing of Hotair Lncrna Restores Gefitinib Sensitivity by Activating Bax/Caspase-3 and Suppressing Tgf-Alpha/Egfr Signaling in Lung Adenocarcinoma. Oncol Lett (2018) 15(3):2829–38. doi: 10.3892/ol.2017.7656

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Xing H, Sun H, Du W. Linc01116 Accelerates Nasopharyngeal Carcinoma Progression Based on Its Enhancement on Myc Transcription Activity. Cancer Med (2020) 9(1):269–77. doi: 10.1002/cam4.2624

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Zhang ZF, Xu HH, Hu WH, Hu TY, Wang XB. Linc01116 Promotes Proliferation, Invasion and Migration of Osteosarcoma Cells by Silencing P53 and Ezh2. Eur Rev Med Pharmacol Sci (2019) 23(16):6813–23. doi: 10.26355/eurrev_201908_18720

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zeng L, Lyu X, Yuan J, Wang W, Zhao N, Liu B, et al. Long Non-Coding RNA Linc01116 Is Overexpressed in Lung Adenocarcinoma and Promotes Tumor Proliferation and Metastasis. Am J Transl Res (2020) 12(8):4302–13.

PubMed Abstract | Google Scholar

24. Wang H, Lu B, Ren S, Wu F, Wang X, Yan C, et al. Long Noncoding RNA LINC01116 Contributes to Gefitinib Resistance in Non-Small Cell Lung Cancer Through Regulating Ifi44. Mol Ther Nucleic Acids (2020) 19:218–27. doi: 10.1016/j.omtn.2019.10.039

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Guo X, Wang Y. Lncrna Tmpo-As1 Promotes Hepatocellular Carcinoma Cell Proliferation, Migration and Invasion Through Sponging Mir-329-3p to Stimulate Foxk1-Mediated Akt/Mtor Signaling Pathway. Cancer Med (2020) 9(14):5235–46. doi: 10.1002/cam4.3046

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Fang S, Liu Z, Guo Q, Chen C, Ke X, Xu G. High Bancr Expression Is Associated With Worse Prognosis in Human Malignant Carcinomas: An Updated Systematic Review and Meta-Analysis. BMC Cancer (2020) 20(1):870. doi: 10.1186/s12885-020-07177-6

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Yu X, Zheng H, Chan MT, Wu WKK. BANCR: A Cancer-Related Long Non-Coding RNA. Am J Cancer Res (2017) 7(9):1779–87.

PubMed Abstract | Google Scholar

28. Lu Y, Sha H, Sun X, Zhang Y, Wu Y, Zhang J, et al. Crnde: An Oncogenic Long Non-Coding RNA in Cancers. Cancer Cell Int (2020) 20:162. doi: 10.1186/s12935-020-01246-3

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Zhang F, Wang H, Yu J, Yao X, Yang S, Li W, et al. LncRNA Crnde Attenuates Chemoresistance in Gastric Cancer via Srsf6-Regulated Alternative Splicing of Picalm. Mol Cancer (2021) 20(1):6. doi: 10.1186/s12943-020-01299-y

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Amin MB, Greene FL, Edge SB, Compton CC, Gershenwald JE, Brookland RK, et al. The Eighth Edition Ajcc Cancer Staging Manual: Continuing to Build a Bridge From a Population-Based to a More "Personalized" Approach to Cancer Staging. CA: Cancer J Clin (2017) 67(2):93–9. doi: 10.3322/caac.21388

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Bruni D, Angell HK, Galon J. The Immune Contexture and Immunoscore in Cancer Prognosis and Therapeutic Efficacy. Nat Rev Cancer (2020) 20(11):662–80. doi: 10.1038/s41568-020-0285-7

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Pagès F, Mlecnik B, Marliot F, Bindea G, Ou FS, Bifulco C, et al. International Validation of the Consensus Immunoscore for the Classification of Colon Cancer: A Prognostic and Accuracy Study. Lancet (2018) 391(10135):2128–39. doi: 10.1016/S0140-6736(18)30789-X

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Donnem T, Kilvaer TK, Andersen S, Richardsen E, Paulsen EE, Hald SM, et al. Strategies for Clinical Implementation of Tnm-Immunoscore in Resected Non-Small Cell Lung Cancer. Ann Oncol (2016) 27(2):225–32. doi: 10.1093/annonc/mdv560

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Bremnes RM, Al-Shibli K, Donnem T, Sirera R, Al-Saad S, Andersen S, et al. The Role of Tumor-Infiltrating Immune Cells and Chronic Inflammation at the Tumor Site on Cancer Development, Progression, and Prognosis. J Thorac Oncol (2011) 6(4):824–33. doi: 10.1097/JTO.0b013e3182037b76

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Farhood B, Najafi M, Mortezaee K. Cd8(+) Cytotoxic T Lymphocytes in Cancer Immunotherapy: A Review. J Cell Physiol (2019) 234(6):8509–21. doi: 10.1002/jcp.27782

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Borst J, Ahrends T, Babala N, Melief CJM, Kastenmuller W. Cd4(+) T Cell Help in Cancer Immunology and Immunotherapy. Nat Rev Immunol (2018) 18(10):635–47. doi: 10.1038/s41577-018-0044-0

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Sharonov GV, Serebrovskaya EO, Yuzhakova DV, Britanova OV, Chudakov DM. B Cells, Plasma Cells and Antibody Repertoires in the Tumor Microenvironment. Nat Rev Immunol (2020) 20(5):294–307. doi: 10.1038/s41577-019-0257-x

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Crotty S. T Follicular Helper Cell Biology: A Decade of Discovery and Diseases. Immunity (2019) 50(5):1132–48. doi: 10.1016/j.immuni.2019.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Bisheshar SK, De Ruiter EJ, Devriese LA, Willems SM. The Prognostic Role of Nk Cells and Their Ligands in Squamous Cell Carcinoma of the Head and Neck: A Systematic Review and Meta-Analysis. Oncoimmunology (2020) 9(1):1747345. doi: 10.1080/2162402X.2020.1747345

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Melaiu O, Chierici M, Lucarini V, Jurman G, Conti LA, De Vito R, et al. Cellular and Gene Signatures of Tumor-Infiltrating Dendritic Cells and Natural-Killer Cells Predict Prognosis of Neuroblastoma. Nat Commun (2020) 11(1):5992. doi: 10.1038/s41467-020-19781-y

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Helby J, Bojesen SE, Nielsen SF, Nordestgaard BG. Ige and Risk of Cancer in 37 747 Individuals From the General Population. Ann Oncol (2015) 26(8):1784–90. doi: 10.1093/annonc/mdv231

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ferastraoaru D, Bax HJ, Bergmann C, Capron M, Castells M, Dombrowicz D, et al. Allergooncology: Ultra-Low Ige, A Potential Novel Biomarker in Cancer-a Position Paper of the European Academy of Allergy and Clinical Immunology (Eaaci). Clin Transl Allergy (2020) 10:32. doi: 10.1186/s13601-020-00335-w

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Platzer B, Elpek KG, Cremasco V, Baker K, Stout MM, Schultz C, et al. Ige/fcepsilonri-Mediated Antigen Cross-Presentation by Dendritic Cells Enhances Anti-Tumor Immune Responses. Cell Rep (2015) 10(9):1487–95. doi: 10.1016/j.celrep.2015.02.015

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Ly D, Zhu CQ, Cabanero M, Tsao MS, Zhang L. Role for High-Affinity Ige Receptor in Prognosis of Lung Adenocarcinoma Patients. Cancer Immunol Res (2017) 5(9):821–9. doi: 10.1158/2326-6066.CIR-16-0392

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Havel JJ, Chowell D, Chan TA. The Evolving Landscape of Biomarkers for Checkpoint Inhibitor Immunotherapy. Nat Rev Cancer (2019) 19(3):133–50. doi: 10.1038/s41568-019-0116-x

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Doroshow DB, Sanmamed MF, Hastings K, Politi K, Rimm DL, Chen L, et al. Immunotherapy in Non-Small Cell Lung Cancer: Facts and Hopes. Clin Cancer Res (2019) 25(15):4592–602. doi: 10.1158/1078-0432.CCR-18-1538

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Garon EB, Hellmann MD, Rizvi NA, Carcereny E, Leighl NB, Ahn MJ, et al. Five-Year Overall Survival for Patients With Advanced Non-Small-Cell Lung Cancer Treated With Pembrolizumab Results From the Phase I Keynote-001 Study. J Clin Oncol (2019) 37(28):2518–27. doi: 10.1200/JCO.19.00934

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: non-small cell lung cancer, long non-coding RNA, immune, prognostic model, tumor microenvironment

Citation: Li Q, Yao L, Lin Z, Li F, Xie D, Li C, Zhan W, Lin W, Huang L, Wu S and Zhou H (2021) Identification of Prognostic Model Based on Immune-Related LncRNAs in Stage I-III Non-Small Cell Lung Cancer. Front. Oncol. 11:706616. doi: 10.3389/fonc.2021.706616

Received: 07 May 2021; Accepted: 24 September 2021;
Published: 20 October 2021.

Edited by:

Paul Takam Kamga, Université de Versailles Saint-Quentin-en-Yvelines, France

Reviewed by:

Jiajia Zhang, Johns Hopkins University, United States
Ryota Kurimoto, Tokyo Medical and Dental University, Japan
Tabatha Gutierrez Prieto, University of São Paulo, Brazil

Copyright © 2021 Li, Yao, Lin, Li, Xie, Li, Zhan, Lin, Huang, Wu and Zhou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Haiyu Zhou,

†The authors share first authorship