- 1Department of Thoracic Surgery, The Second Xiangya Hospital, Central South University, Changsha, China
- 2Hunan Key Laboratory of Early Diagnosis and Precision Therapy of Lung Cancer, The Second Xiangya Hospital, Central South University, Changsha, China
Background: Emerging scientific evidence has shown that long non-coding RNAs (lncRNAs) exert critical roles in genomic instability (GI), which is considered a hallmark of cancer. To date, the prognostic value of GI-associated lncRNAs (GI-lncRNAs) remains largely unexplored in lung adenocarcinoma (LUAC). The aims of this study were to identify GI-lncRNAs associated with the survival of LUAC patients, and to develop a novel GI-lncRNA-based prognostic model (GI-lncRNA model) for LUAC.
Methods: Clinicopathological data of LUAC patients, and their expression profiles of lncRNAs and somatic mutations were obtained from The Cancer Genome Atlas database. Pearson correlation analysis was conducted to identify the co-expressed mRNAs of GI-lncRNAs. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes enrichment analyses were conducted to determine the main biological function and molecular pathways of the differentially expressed GI-lncRNAs. Univariate and multivariate Cox proportional hazard regression analyses were performed to identify GI-lncRNAs significantly related to overall survival (OS) for construction of the GI-lncRNA model. Kaplan–Meier survival analysis and receiver operating characteristic curve analysis were performed to evaluate the predictive accuracy. The performance of the newly developed GI-lncRNA model was compared with the recently published lncRNA-based prognostic index models.
Results: A total of 19 GI-lncRNAs were found to be significantly associated with OS, of which 9 were identified by multivariate analysis to construct the GI-lncRNA model. Notably, the GI-lncRNA model showed a prognostic value independent of key clinical characteristics. Further performance evaluation indicated that the area under the curve (AUC) of the GI-lncRNA model was 0.771, which was greater than that of the TP53 mutation status and three existing lncRNA-based models in predicting the prognosis of patients with LUAC. In addition, the GI-lncRNA model was highly correlated with programed death ligand 1 (PD-L1) expression and tumor mutational burden in immunotherapy for LUAC.
Conclusion: The GI-lncRNA model was established and its performance was found to be superior to existing lncRNA-based models. As such, the GI-lncRNA model holds promise as a more accurate prognostic tool for the prediction of prognosis and response to immunotherapy in patients with LUAC.
Introduction
Lung cancer is one of the most common malignancies and causes the largest number of cancer-related deaths globally (Bray et al., 2018; Siegel et al., 2019). Lung cancer is histologically heterogeneous with lung adenocarcinoma (LUAC) as the most common pathological subtype, accounting for approximately 40% of all lung cancer cases (Denisenko et al., 2018). Although substantial progress and advances have been made in both the diagnosis and treatment of LUAC (e.g., surgical resection, immunotherapy, chemotherapy, targeted therapy), which have greatly improved the clinical outcome of LUAC patients, the prognosis of LUAC is still far from satisfactory with a 5-years survival rate as low as 21% (Macheleidt et al., 2018). Therefore, the development of a reliable prognostic tool is needed to precisely identify high-risk patients, and the implementation of optimal interventions is of great significance to improve patient prognosis in LUAC.
Genomic instability (GI), usually referred to as the high frequency of genetic mutations, is a hallmark of cancer, and these mutations allow cancer cells to adapt to environmental stress and drive the development of more aggressive cancer cells (Lengauer et al., 1998; Negrini et al., 2010). Abnormal transcriptional or post-transcriptional regulation potentially leads to gene mutations and chromosomal aberrations such as cell cycle checkpoints, DNA replication, DNA repair, mitosis, and epigenetic regulation (Soca-Chafre et al., 2019; Tam et al., 2019). For example, the aberrant expression of cell cycle-associated genes, such as cyclins and cyclin-dependent kinases, causes chromosomal change and promotes tumor progression (Broustas and Lieberman, 2014). A previous study revealed that nearly 2.5% of cancers can be attributed to mutations in DNA repair genes (Parry et al., 2017). Notably, lung cancer has the second highest frequency of somatic mutations (Kandoth et al., 2013), with an abnormal number of chromosomes, namely aneuploidy, which is detected in more than 60% of patients with non-small cell lung cancer, and genomic duplication occurs in more than 40% of patients with lung cancer (Burgess et al., 2020). These previous findings indicate the pivotal role of GI in the development and progression of lung cancer.
Long non-coding RNAs (lncRNAs), an emerging class of ncRNAs (Mattick and Rinn, 2015), exert regulatory roles in both genomic stability and GI (Liu, 2016; Nair et al., 2020). For example, the lncRNA GUARDIN promotes the expression of telomeric repeat-binding factor 2 by competitively binding to microRNA-23a, thus maintaining genomic stability (Hu et al., 2018). LncRNA LINC00657 inhibits mitosis, DNA repair, and DNA replication via binding to PUMILIO protein, which is essential for the maintenance of genomic stability (Elguindy et al., 2019). In contrast, lncRNA CCAT2 promotes carcinogenesis and GI (Chen et al., 2020). Given that lncRNAs have unique roles in maintaining genomic stability and promoting GI, we hypothesized that GI-associated lncRNAs may have prognostic value in LUAC.
In this study, we comprehensively analyzed the gene expression profiles, somatic mutations, and corresponding clinical data of LUCA patients from The Cancer Genome Atlas (TCGA) database, with the aim of developing a novel GI-lncRNA prognostic model to better predict the clinical outcomes of LUAC patients.
Materials and Methods
Data Acquisition and Processing
The level 3 transcriptome profiles of 535 LUAC tissues and 59 histologically normal tissues, somatic mutation profiles of 561 LUAC samples, and corresponding clinicopathological data of 522LUAC cases were acquired from TCGA database. For transcriptome profiles, mRNA data and lncRNA data were separated into a mRNA expression matrix and lncRNA expression matrix. For somatic mutation profiles, total mutation frequency of each case and frequency of the mutant gene were computed in all LUAC cases. After the LUAC cases with a survival time less than 30 days or incomplete follow-up information were excluded from further analyses, 490 patients with LUAC were randomly allocated into two cohorts: a training cohort (n = 246) and testing cohort (n = 244). The training cohort was used to identify GI-lncRNAs independently associated with overall survival (OS) for development of the prognostic index model, while the testing and total cohorts were used for validation of the newly developed model. The clinicopathological characteristics of patients in the training and testing cohorts are presented in Table 1.
Identification of GI-Associated lncRNAs in LUAC
To identify GI-lncRNAs in LUAC, we initially computed the total somatic mutations of each case and integrated them with the lncRNA expression matrix by the sample names. According to the mutator hypothesis-derived computational frame as previously described (Bao et al., 2019), we ranked the LUAC patients by total somatic mutations, and defined the top 25% of cases as the genome unstable-like (GU-like) group and the last 25% cases as the genome stable-like (GS-like) group. Then we performed differential expression analyses between the GU-like group and GS-like group using the “limma” R package, and lncRNAs with a false discovery rate (FDR) less than 0.05 and |logFC| > 1 were identified as GI-lncRNAs. A heatmap and volcano plot were constructed to visualize these differentially expressed lncRNAs using the “igraph” R package.
Hierarchical Clustering Analysis of the Differentially Expressed lncRNAs
We performed hierarchical clustering analysis of GI-lncRNAs in the LUAC cases. Briefly, “sparcl,” “pheatmap,” and “limma” R packages were used to compute Euclidean distances, and LUAC cases were stratified into two clusters. Then the clusters matrix was integrated with total somatic mutations count matrix. By comparing the median mutation counts of two clusters, the cluster with lower mutation count was defined as the GS-like cluster, while the other was defined as the GU-like cluster. Then the difference in total somatic mutation counts between two clusters was explored using the “limma” R package. Given that ubiquilin-4 (UBQLN4) has been demonstrated to be a driver gene of GI and is overexpressed in malignant tumors (Jachimowicz et al., 2019), we compared the expression levels of UBQLN4 between the two clusters in this study.
Gene Co-expression Network and Functional Enrichment Analysis
To reveal the potential biological function and molecular pathways of GI-lncRNAs, we performed gene co-expression analysis to identify their co-expressed mRNAs. Pearson correlation analysis was performed to identify mRNAs that were co-expressed with the GI-lncRNAs using the “limma” R package by the threshold of correlation coefficient >0.3 and p < 0.05. The top 10 co-expressed mRNAs were selected for subsequent analyses of the gene co-expression network using the “igraph” R package. Then Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed for the co-expressed mRNAs by the “clusterProfiler” R package. Molecular pathways with p < 0.05 were considered significantly enriched.
Construction and Validation of the GI-Associated lncRNA Signature
The training cohort was used to identify GI-lncRNAs independently associated with OS for construction of the prognostic index model. In brief, univariate Cox regression analysis was performed to identify GI-lncRNAs significantly related to OS (p < 0.05) using the “survival” R package. These GI-lncRNAs were subsequently subjected to multivariate Cox proportional hazard regression analysis to construct the optimal prognostic model (termed GI-lncRNAsSig) using the “survival” and “survminer” R packages. Patients’ risk score of the GI-lncSig was calculated as follows: risk score = Σ (ExpmRNAn × βmRNAn). According to the median risk score, patients were stratified into high- and low-risk groups. Kaplan–Meier survival analysis and receiver operating characteristic (ROC) analysis were used to evaluate the performance of the GI-lncRNAsSig for the prediction of OS using the “survminer” and “survivalROC” R packages. Finally, the testing cohort was used to validate the prognostic performance of the developed GI-lncRNAsSig in patients with LUAC.
Clinical Risk Stratification and Independent Prognostic Value of the GI-lncRNA-Based Model
To explore the applicability of the GI-lncRNA-based prognostic model, we performed clinical risk stratification in the total TCGA cohort. In brief, LUAC patients were stratified into subgroups based on clinical characteristics including age (≤65 and >65), gender (female and male), tumor stage (I–II and III–IV), pathologic T classification (T1-2 and T3–4), pathologic N classification (N0 and N1–3), and pathologic M classification (M0 and M1). Each subgroup was further stratified into high- and low-risk groups according to the median risk score of the newly developed GI-lncRNA-based model. Kaplan–Meier survival analysis was performed to explore the survival difference between the high- and low-risk groups. Multivariate Cox regression analysis was performed to determine whether the GI-lncRNA-based model could have independent prognostic value in the training, testing, and total cohorts.
Statistical Analyses
Statistical analyses were conducted using R software version 4.0.2. The Mann–Whitney test was used to compare quantitative data between different groups, and the chi-squared test was used to compare categorical data between groups. p < 0.05 (if not specified) was considered statistically significant.
Results
Identification of GI-lncRNAs in LUAC
We initially performed comprehensive analyses of the somatic mutation profiles of 561 LUAC cases, and the resulting data are presented in Supplementary Tables S1, S2 and Supplementary Figure S1. Somatic mutations were detected in the majority of LUAC tumor tissues (449/561, 88.95%) affecting 18,498 genes. The tumor protein p53 (TP53) (n = 272), giant-muscle filament titin (n = 265), mucoprotein-16 (n = 244), ryanodine receptor 2 (n = 228), and CUB and sushi multiple domains 3 (n = 215) were the top five most frequently mutated genes. After integrating the lncRNA expression matrix with the somatic mutations by the sample names, the top 25% cases with the highest mutation frequency were defined as the GU-like group, and the last 25% cases with the lowest mutation frequency were defined as the GS-like group. Differential analyses with the heatmap and volcano plot revealed that the expression levels of 138 lncRNAs were significantly different between the GU-like group and GS-like group, in which 59 were significantly upregulated and 79 were significantly downregulated (FDR <0.05, |logFC| > 1) (Figure 1). Based on the 138 differentially expressed lncRNAs, 490 LUAC patients were stratified into two clusters by hierarchical clustering analysis (Figure 2A). The cluster with the higher somatic mutation count was defined as the GU-like cluster, while the other was defined as the GS-like cluster. As shown in Figure 2B, the total somatic mutation count of the GU-like cluster was significantly higher than that of the GS-like cluster (p < 0.001). Notably, the expression of UBQLN4, a driver gene of GI, was significantly upregulated in the GU-like cluster (p < 0.001; Figure 2C). These results suggest that the 138 differentially expressed lncRNAs could be considered GI-lncRANs in LUAC.
FIGURE 1. Heatmap and volcano plot of differentially expressed lncRNAs in LUAC. (A) Heatmap of differentially expressed lncRNAs in the GS-like group versus the GU-like group. The upregulated lncRNAs are denoted in red, and the downregulated lncRNAs are indicated in green; (B) Volcano plot of differentially expressed lncRNAs in the GS-like group versus the GU-like group. The upregulated lncRNAs are denoted in red, the downregulated lncRNAs are indicated in green, and the unaltered lncRNAs are denoted in black.
FIGURE 2. Functional annotation of GI-lncRNAs in LUAC patients. (A) Hierarchical clustering analysis of differentially expressed lncRNAs in the two clusters: GS-like cluster and GU-like cluster; (B) Somatic mutation counts in the GS-like cluster and GU-like cluster; (C) Comparison of expression levels of UBQLN4, a driver gene of GI, between the GS-like cluster and the GU-like cluster; (D) Network analysis of the relationship between GI-lncRNAs and mRNAs; (E) GO analysis; (F) KEGG pathway analysis.
Functional Annotation and Molecular Pathway Analysis of GI-LncRNAs in LUAC
To reveal the potential biological functions and disrupted molecular pathways of the 138 GI-lncRNAs, we performed GO and KEGG enrichment analyses for their co-expressed mRNAs, which were identified by Pearson correlation analysis using a threshold of correlation coefficient >0.3 and p < 0.05. As shown in Figure 2D, a lncRNA-mRNA co-expression network was conducted using the GI-lncRNAs and their top 10 co-expressed mRNAs by the “igraph” R package. GO enrichment analysis indicated that the biological processes of the co-expressed mRNAs were mainly related to the formation of GI including cell cycle regulation, double-strand break repair, and DNA integrity checkpoint (Figure 2E). KEGG enrichment analysis indicated that these co-expressed mRNAs were enriched in molecular pathways associated with GI including cell cycle regulation, p53 signalling pathway, transcriptional dysregulation, nucleotide excision repair, and base excision repair (Figure 2F). Interestingly, programmed death-ligand 1 (PD-L1) expression and the programmed cell death protein 1 (PD-1) checkpoint pathway was enriched in these co-expressed mRNAs. Together, these results suggest that these 138 lncRNAs are closely associated with GI, and their aberrant expression may disrupt the lncRNA-mRNA regulatory network, thus promoting GI.
Construction and Validation of the GI-lncRNA-Based Prognostic Model
To identify GI-lncRNAs with prognostic value, univariate Cox regression analysis was performed in the training cohort. Among the 138 GI-lncRNAs, 19 were significantly associated with OS, including 9 risk factors for OS and 10 protective factors (Figure 3). Furthermore, multivariate Cox proportional hazard regression analysis identified nine GI-lncRNAs to construct the optimal prognostic model. The risk score to predict OS in LUAC patients was calculated as follows: risk score = (−0.192 * expression of LINC02159) + (−0.059 * expression of AC025154.2) + (0.049 * expression of LINC01671) + (0.051 * expression of FAM83A-AS1) + (−0.113 * expression of AC125603.1) + (0.423 * expression of AC021218.1) + (−0.409 * expression of AF131215.5) + −0.126 * expression of RHOXF1-AS1) + (0.151 * expression of LINC01116).
FIGURE 3. GI-lncRNAs associated with the prognosis of LUAC patients in the univariate Cox regression analysis. A total of 19 GI-lncRNAs were significantly associated with OS, including 9 risk factors for OS and 10 protective factors.
The performance of the newly constructed prognostic model was initially evaluated in the training cohort. After LUAC patients were stratified into high- and low-risk groups based on the median risk score of the prognostic model (Supplementary Table S3), Kaplan–Meier survival analysis showed a significant difference in OS between the groups (Figure 4A), with significantly poorer prognosis found in the high-risk group compared to the low-risk group (p < 0.001; Figure 4B). ROC curve analysis revealed that the area under the curve (AUC) value of the newly developed model was 0.771, indicating good performance in predicting the prognosis of patients with LUAC (Figure 4C). The prognostic model was validated using the testing cohort and total cohort in which LUAC patients were stratified into high- and low-risk groups in accordance with the median risk score of the prognostic model (Supplementary Tables S4, S5). Kaplan–Meier curve analysis indicated a significantly worse OS in the high-risk group (p = 0.024; Figure 4D) with the poorer OS correlated with higher risk scores (Figure 4E). The AUC value was 0.747. Similar results were obtained in the total cohort (Figures 4G–I). These results indicated that the performance of the newly develop GI-lncRNA-based prognostic model was well validated for the prediction of OS in patients with LUAC.
FIGURE 4. Performance evaluation and validation of the newly developed GI-lncRNA-based model for the prediction of OS in LUAC patients. Kaplan–Meier analysis of OS between the high-risk group and low-risk group in (A) the training cohort, (D) the testing cohort, and (G) the total cohort; The relationship between OS and risk-risk scores of the new prognostic model in the (B) the training cohort, (E) the testing cohort, and (H) the total cohort; The receiver operating curves for the new prognostic model in (C) the training cohort, (F) the testing cohort, and (I) the total cohort.
Determination of the Independent Prognostic Value of the GI-lncRNA-Based Model for LUAC
To determine if the newly developed GI-lncRNA model could independently predict prognosis, we initially performed risk stratification analysis based on key clinical characteristics of LUAC patients and risk scores of the GI-lncRNA model. Patients were stratified into different subgroups by key clinical characteristics, including age (≤65 and >65), gender (female and male), tumor stage (I-II and III-IV), pathologic T classification (T1-2 and T3-4), pathologic N classification (N0 and N1-3), and pathologic M classification (M0 and M1). Patients in each subgroup were further stratified into the high- and low-risk groups according to the risk scores of the GI-lncRNA model. As shown in Figure 5, Kaplan–Meier analysis showed that the OS of patients in the high-risk group was significantly worse than that of the low-risk group in all subgroups (p < 0.05), except for patients in the M1 subgroup and patients in the racial group of African American, which could be due to the relatively small sample size of the two subgroups. However, the OS of the low-risk group was better than that of the high-risk group in both subgroups, whereas there was no significant difference. To further determine whether the GI-lncRNA-based prognostic model could independently predict OS, we performed univariate Cox regression and multivariate Cox regression analyses in the training cohort, testing cohort, and total cohort. As illustrated in Table 2, the GI-lncRNA-based prognostic model was independent of key clinical characteristics to predict OS in LUAC patients in the training, testing, and total cohorts (all p < 0.05).
FIGURE 5. Kaplan–Meier analysis of OS between the high-risk and low-risk groups in LUAC patients with different clinical characteristics. Kaplan–Meier curves of LUAC patients with different clinical characteristics, including (A) age (≤65/>65); (B) gender (female/male); (C) race (Caucasian American/African American); (D) stage of lung cancer (Stage I-II/III-IV); (E) pathologic T stage (T1-2/T3-4); (F) pathologic N stage (N0/N1-3); (G) pathologic M stage (M0/M1).
TABLE 2. Univariate and multivariate analyses of risk factors associated with the prognosis of LUAC patients.
Performance Comparison of the Newly Developed GI-lncRNA Model With TP53 Mutation Status and Existing LncRNA-Based Models for Predicting the OS of LUAC Patients
In this study, TP53 was found to be the most frequently mutated gene in LUAC (Supplementary Figure S1), therefore, we analyzed the TP53 mutation pattern between the high- and low-risk groups as stratified by the newly developed GI-lncRNA model. As shown in Figures 6A–C, the frequency of TP53 mutation was different in the high- and low-risk groups. In the training cohort, TP53 mutation was detected in 60.7% of patients in the high-risk group, which was significantly higher than the 41.7% in the low-risk group (p = 0.005; Figure 6A). Similar results were obtained in the testing cohort (Figure 6B) and total cohort (Figure 6C). Given that TP53 mutation status is closely related to GI and it has been proposed as a biomarker with prognostic value in lung cancer (Zhou et al., 2019; Freudenstein et al., 2020), we further evaluated the performance between the newly developed GI-lncRNAs and TP53 mutation status in predicting the OS of patients with LUAC. According to the TP53 mutation status and risk score of the GI-lncRNA model, LUAC patients were further classified into four groups: TP53 Wild/high risk group, TP53 Wild/low risk group, TP53, Mutation/high risk group, and TP53 Mutation/low risk group. As illustrated in Figure 6D, Kaplan–Meier analysis indicated the OS was significantly different among these four groups (p < 0.001). In the TP53 wild-type (WT) group, patients with a high risk score (referred to as TP53 Wild/high risk) had a worse OS than those with low risk score (referred to as Wild/low risk). In the TP53 mutation group, patients with a high risk score (referred to as TP53 Mutation/high risk) also had a worse OS than those with a low risk score (referred to as TP53 Mutation/low risk). In the low-risk group, patients with WT TP53 (TP53 Wild/high risk) had a better OS than those with TP53 mutation status (TP53 Mutation/high risk). However, in the high-risk group, the OS of patients with TP53 wild-type (TP53 Wild/high risk) was similar to those with TP53 mutation type (TP53 Mutation/high risk), indicating that TP53 mutation status failed to discriminate the OS of patients in the high-risk group. Interestingly, the OS of the TP53 Mutation/low risk group was better than that of the TP53 Wild/high risk group. Thus, these findings indicate the overall better prognostic value of the GI-lncRNA model than the TP53 mutation status.
FIGURE 6. Comparison of the GI-lncRNA-based model and TP53 status for the prediction of OS in LUAC. Comparison of TP53 mutation status between the high-risk group and low-risk group stratified by the risk scores of the GI-lncRNA-based model in (A) the training cohort, (B) the testing cohort, and (C) the total cohort; (D) Kaplan–Meier analyses of LUAC patients in the following four groups stratified by the TP53 mutation status and risk score of the GI-lncRNA-based model: TP53 Wild/high risk group, TP53 Wild/low risk group, TP53 Mutation/high risk group, and TP53 Mutation/low risk group.
We further conducted a performance comparison of the GI-lncRNAs with three recently published lncRNA-based prognostic models, including 7 lncRNA signatures from Li et al. (referred to as LiLncSig) (Li et al., 2020), 7 lncRNA signatures from Zhou et al. (referred to as JinLncSig) (Jin et al., 2020), and 13 lncRNA signatures from Zhou et al. (referred to ZhouLncSig) (Zhou et al., 2020). As shown in Figure 7, the AUC for the newly developed GI-lncRNA model was 0.757, which was significantly greater than that of LiLncSig (AUC, 0.653), JinLncSig (AUC, 0.696), and ZhouLncSig (AUC, 0.689). These data indicate that the newly developed GI-lncRNA model is superior to the three published lncRNA-based models.
FIGURE 7. Performance comparison of the GI-lncRNA-based model with existing lncRNA-based prognostic models for LUAC. The receiver operating curves for the new prognostic model and the recently published lncRNA-based models for predicting the prognosis of LUAC patients.
Value of the GI-lncRNA Model in Predicting Response to Immunotherapy
KEGG enrichment analysis showed that PD-L1 expression and the PD-1 checkpoint pathway was enriched in the co-expressed mRNAs of GI-lncRNAs; thus, we further examined the expression pattern of PD-L1 and PD-1 in patients in different risk groups to explore the potential value of the newly developed GI-lncRNA model in predicting response to immunotherapy. As shown in Figures 8A,B, despite the significantly worse OS of patients in the high-risk and high PD-1 (PDCD1) group (p < 0.001), there was no significant difference in PD-1 expression between the high- and low-risk groups (p = 0.338), whereas the expression level of PD-L1 (CD274) in the high-risk group was significantly higher than that in the low-risk group (p = 0.006; Figure 8C) and the OS of the high-risk and high PD-L1 groups was significantly worse than that of the low-risk and PD-L1 groups (p < 0.001; Figure 8D). These findings indicate that the GI-lncRNA model can be used to predict the response to anti-PD-L1 immunotherapy. Combining tumor mutational burden (TMB) and PD-L1 expression greatly enhances the predictive power of response to immunotherapy efficacy (Hellmann et al., 2018). Given the hypothesis that patients with increased GI may have a higher frequency of somatic mutations and TMB, we further compared the difference in TMB between the high- and low-risk groups. Interestingly, the TMB of the high-risk group was significantly higher than that of the low-risk group (p < 0.001; Figure 8E), and Kaplan–Meier analysis indicated a significantly worse OS in patients with high TMB and in the high-risk group (p < 0.001; Figure 8F). In addition, we assessed the correlation between the GI-lncRNA model and PD-L1/TMB in the training cohort and testing cohort. As shown in Supplementary Figure S1, PD-L1/TMB was significantly higher in both cohorts. Therefore, LUAC patients in the high-risk group may benefit more from immunotherapy compared with those in the low-risk group, suggesting the value of the newly developed GI-lncRNA model for predicting the response to immunotherapy in patients with LUAC.
FIGURE 8. Value of the GI-lncRNA model for predicting the response to immunotherapy. Comparison of expression levels of (A) PDCD1, (C) CD274 between the low-risk and high-risk groups stratified by the risk scores of the GI-lncRNA-based model; Kaplan–Meier curve analysis of OS in (B) the low-risk group with low PDCD1 versus high-risk group with high PDCD1, and (D) low-risk group with low CD274 versus high-risk group with high CD274; (E) Comparison of difference in tumor mutational burden (TMB) between the high- and low-risk groups; (F) Kaplan–Meier analysis of OS in the low-risk group with low TMB versus high-risk group with high TMB.
Discussion
The prognostic value of GI-lncRNAs in patients with LUAC remains largely unexplored. The key novel findings of this study were as follows. Nineteen GI-lncRNAs associated with OS were identified in patients with LUAC. Nine GI-lncRNAs were significantly correlated with OS, which were used to generate the newly developed GI-lncRNA-based prognostic model. The GI-lncRNA model performed well with an AUC of 0.771, which was greater than the AUCs of the TP53 mutation status and three reported lncRNA-based models in predicting prognosis of patients with LUAC. The GI-lncRNA model showed a prognostic value independent of key clinical characteristics. The GI-lncRNA model was strongly correlated with PD-L1 and TMB, suggesting its value in predicting the response to immunotherapy in patients with LUAC. Together, these findings support the newly established GI-lncRNA model as a potentially better prognostic approach to predicting the prognosis and response to immunotherapy in LUAC patients.
Despite great advances in the diagnosis and treatment of lung cancer, OS is still considerably low and lung cancer remains the most common cause of cancer-related deaths globally (Lee et al., 2018). Traditional clinical parameters such as tumor-node-metastasis stage, tumor size, regional lymph node metastasis, and distant metastasis have long been used for prognosis prediction in lung cancer. However, these conventional prognostic factors are not always convincible, as patients with the same clinical and pathological characteristics may encounter completely different clinical outcomes (Liu et al., 2019). Therefore, a prognostic approach with high accuracy may assist clinicians in making optimal treatment decisions to improve the patient survival. The past several decades have witnessed great progress in understanding the biological mechanisms of tumorigenesis, development, and progression. One of the major breakthroughs was the involvement of GI in tumorigenesis and the therapeutic response (Negrini et al., 2010). Lung cancer has the second highest somatic mutation burden, just second to melanoma, indicating the potential prognostic and diagnostic value of GI (Kandoth et al., 2013). An increasing number of studies have suggested that the aberrant expression of lncRNAs plays a key role in GI (Liu, 2016; Nair et al., 2020). Therefore, identifying GI-lncRNAs in the whole genome and exploring their prognostic value are of great significance.
In this study, we identified 138 GI-lncRNAs in LUAC using the computational frame as previously reported by Bao et al. (2019). GO analysis showed that the biological processes of the co-expressed mRNAs of 138 lncRNAs were enriched in pathways of cell cycle regulation and DNA damage response, such as cell cycle checkpoint and double-strand break repair. Cell cycle checkpoints are the dynamically surveillance mechanism of major events in cell cycle that ensure the order, integrity and fidelity of cell replication, and dysregulation of cell cycle checkpoints that often accompany GI (Anand et al., 2020). For example, p53, a major cell cycle checkpoint of G2/M phase, plays an important role in DNA replication by halting the cell cycle at G2 phase and allowing the repair mechanism to restore genomic stability (Williams and Schumacher, 2016). Double-strand break repair is another important mechanism for maintaining genomic stability (Terasawa et al., 2014). In light of KEGG analysis, the co-expressed mRNAs were also enriched in molecular pathways of cell cycle regulation and DNA damage repair (e.g., cell cycle, transcriptional dysregulation in cancer, p53 signaling pathway, and nucleotide excision repair) in association with maintenance of genomic stability. These results further confirmed that these 138 lncRNAs were closely linked to GI. In addition, the co-expressed mRNAs were enriched in the MAPK signaling pathway, which is involved in cancer invasion and metastasis (Wei et al., 2017; Wen et al., 2019). Therefore, we postulated that the abnormal expression of these differentially expressed lncRNAs may disrupt the lncRNA-mRNA regulatory network, thereby leading to GI and promoting the invasive and metastatic capacity of tumor cells. Then we identified GI-lncRNAs related to OS and constructed a GI-lncRNA model consisting of nine lncRNAs to predict clinical outcomes in the training cohort. Among the nine GI-lncRNAs used for construction of the GI-lncRNA model, the biological functions of FAM83A-AS1 and LNC01116 have been well illustrated in LUAC. LncRNA FAM83A-AS1 promotes LUAC cell proliferation, migration, invasion and the epithelial–mesenchymal transition (EMT) by competitively combining miR-150-5p with MMP14, and Wang et al. reported that FAM83A-AS1 promotes the progression of LUAC by enhancing its pre-mRNA FAM83A via the ERK signaling pathway (Xiao et al., 2019). Zeng et al. reported that the overexpression of LINC01116 contributes to tumor proliferation and metastasis of LUAC cells, and Wang et al. reported that LINC01116 contributes to cisplatin resistance via the EMT process (Wang et al., 2020; Zeng et al., 2020). Although the biological functions of AF131215.5 and RHOXF1-AS1 have not been reported in cancer, our results are in accordance with previous findings showing that these two lncRNAs are associated with the OS of LUAC (Hou and Yao, 2021). The remaining five lncRNAs have not been reported in lung cancer or other human cancers. Their prognostic value and biological function need to be further investigated in future research. The GI-lncRNA model stratified LUAC patients into high- and low-risk groups with significantly different OS, and ROC curve analysis showed the high sensitivity and specificity of the GI-lncRNA model, which were further confirmed in the testing cohort and total cohort. Stratification analysis showed that the GI-lncRNA model was applicable for all clinical subgroups, and multivariate Cox regression analysis revealed that the GI-lncRNA model was an independent prognostic factor for OS in the training, testing, and total cohorts. These results indicate that the GI-lncRNA model may be a promising non-invasive biomarker for OS prediction in LUAC.
In the functional enrichment analysis, we found that the co-expressed mRNAs were enriched in the PD-L1 expression and PD-1 checkpoint pathway, and further analyses revealed higher PD-L1 expression with a significantly worse OS in the high-risk group, indicating that GI-lncRNAs may also have regulatory effects on PD-L1 expression. PD-L1 contributes to the immune evasion of tumor cells by binding to PD-1 and negatively modulating T-cell receptor signaling (Blank et al., 2005; Pardoll, 2012). In recent years, unprecedented achievements have been made in immune checkpoint inhibitors (ICIs) targeting PD-1 or its ligand PD-L1 such as nivolumab, atezolizumabm and pembrolizumab (Doroshow et al., 2021). In fact, PD-L1 expression is the only test approved by the U.S. Food and Drug Administration (FDA) for ICI first-line treatment decision-making in lung cancer (Borghaei et al., 2015). However, PD-L1 expression to predict immunotherapy has limitations (e.g., variability and intra-tumor heterogeneity) (McLaughlin et al., 2016). Moreover, patients with low or no expression of PD‐L1 may also have favorable responses to ICIs (Frigola et al., 2021). TMB is another promising biomarker for ICI response, which was approved by FDA for the treatment decision of pembrolizumab in 2020 (Subbiah et al., 2020). Though TMB and PD-L1 expression are unrelated, greater benefit was found in patients with high TMB and high PD-L1 expression treated with anti-PD-1 and anti-PD-L1 agents in lung cancer, indicating synergistic association of the two independent biomarkers (Carbone et al., 2017; Peters et al., 2017; Hellmann et al., 2018). Based on the hypothesis that accumulating genetic alterations resulting from GI may lead to a higher TMB, we further analyzed the pattern of TMB between different risk groups, and found a significantly higher TMB in the high-risk group. According to the above results, we infer that the GI-lncRNAs may have potential for selecting LUAC patients who will benefit more from immunotherapy.
This study had several limitations. First, the GI-lncRNA model was developed and validated with retrospective data from TCGA database. Although we tried to further validate it in the Gene Expression Omnibus (GEO) database, some of the lncRNAs in the GI-lncRNA model could not be found in the GEO datasets due to a limited number of lncRNAs conserved in the Affymetrix platform. Therefore, a prospective study in the real world will be required to verify its clinical application value. Second, the findings in this study could not explain how these lncRNAs affect GI and malignant behaviors of tumor cells, for which future in-depth mechanistic investigation is warranted to illustrate their biological function and underlying mechanism. Third, it is worth noting that the correlation between the GI-lncRNA model and TMB/PD-L1 expression was not experimentally verified in this study.
In conclusion, a novel prognostic model based on a panel of GI-lncRNAs was established in this study. Notably, the new prognostic model exhibits high prognostic accuracy and overall better performance compared with the reported lncRNA-based prognostic models. In addition, the GI-lncRNA model may help predict which LUAC patient may have a better response to immunotherapy. Therefore, this newly constructed prognostic model may assist oncologists/surgeons in navigating optimal treatment plans, thereby enhancing survival and eventually improving the care of LUAC patients.
Data Availability Statement
Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/repository.
Author Contributions
Conception and design: XW and GT. Data collection: XP, BH, and PZ. Data analysis: GT, QC, ZZ and SS. Manuscript writing: GT and WP. All authors contributed to manuscript revision, read, and approved the submitted version.
Funding
This work was supported by the Hunan Provincial Key Area R&D Programmes (No. 2019SK2253); and Science and Technology Innovation Program of Hunan Province (No. 2020SK53424).
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.
The reviewer (JL) declared a shared parent affiliation with the authors to the handling editor at the time of the review.
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.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.720013/full#supplementary-material
References
Anand, S. K., Sharma, A., Singh, N., and Kakkar, P. (2020). Entrenching Role of Cell Cycle Checkpoints and Autophagy for Maintenance of Genomic Integrity. DNA Repair 86, 102748. doi:10.1016/j.dnarep.2019.102748
Bao, S., Zhao, H., Yuan, J., Fan, D., Zhang, Z., Su, J., et al. (2019). Computational Identification of Mutator-Derived lncRNA Signatures of Genome Instability for Improving the Clinical Outcome of Cancers: a Case Study in Breast Cancer. Brief Bioinform. 21 (5), 1742–1755. doi:10.1093/bib/bbz118
Blank, C., Gajewski, T. F., and Mackensen, A. (2005). Interaction of PD-L1 on Tumor Cells with PD-1 on Tumor-specific T Cells as a Mechanism of Immune Evasion: Implications for Tumor Immunotherapy. Cancer Immunol. Immunother. 54 (4), 307–314. doi:10.1007/s00262-004-0593-x
Borghaei, H., Paz-Ares, L., Horn, L., Spigel, D. R., Steins, M., Ready, N. E., et al. (2015). Nivolumab versus Docetaxel in Advanced Nonsquamous Non-small-Cell Lung Cancer. N. Engl. J. Med. 373 (17), 1627–1639. doi:10.1056/NEJMoa1507643
Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., and Jemal, A. (2018). Global Cancer Statistics 2018: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: A Cancer J. Clinic. 68 (6), 394–424. doi:10.3322/caac.21492
Broustas, C. G., and Lieberman, H. B. (2014). DNA Damage Response Genes and the Development of Cancer Metastasis. Radiat. Res. 181 (2), 111–130. doi:10.1667/RR13515.1
Burgess, J. T., Rose, M., Boucher, D., Plowman, J., Molloy, C., Fisher, M., et al. (2020). The Therapeutic Potential of DNA Damage Repair Pathways and Genomic Stability in Lung Cancer. Front. Oncol. 10, 1256. doi:10.3389/fonc.2020.01256
Carbone, D. P., Reck, M., Paz-Ares, L., Creelan, B., Horn, L., Steins, M., et al. (2017). First-Line Nivolumab in Stage IV or Recurrent Non-small-Cell Lung Cancer. N. Engl. J. Med. 376 (25), 2415–2426. doi:10.1056/NEJMoa1613493
Chen, B., Dragomir, M. P., Fabris, L., Bayraktar, R., Knutsen, E., Liu, X., et al. (2020). The Long Noncoding RNA CCAT2 Induces Chromosomal Instability through BOP1-AURKB Signaling. Gastroenterology 159 (6), 2146–2162. doi:10.1053/j.gastro.2020.08.018
Denisenko, T. V., Budkevich, I. N., and Zhivotovsky, B. (2018). Cell Death-Based Treatment of Lung Adenocarcinoma. Cell Death Dis. 9 (2), 117. doi:10.1038/s41419-017-0063-y
Doroshow, D. B., Bhalla, S., Beasley, M. B., Sholl, L. M., Kerr, K. M., Gnjatic, S., et al. (2021). PD-L1 as a Biomarker of Response to Immune-Checkpoint Inhibitors. Nat. Rev. Clin. Oncol. 18, 345–362. doi:10.1038/s41571-021-00473-5
Elguindy, M. M., Kopp, F., Goodarzi, M., Rehfeld, F., Thomas, A., Chang, T.-C., et al. (2019). PUMILIO, but Not RBMX, Binding Is Required for Regulation of Genomic Stability by Noncoding RNA NORAD. Elife 8, e48625. doi:10.7554/eLife.48625
Freudenstein, D., Litchfield, C., Caramia, F., Wright, G., Solomon, B. J., Ball, D., et al. (2020). TP53 Status, Patient Sex, and the Immune Response as Determinants of Lung Cancer Patient Survival. Cancers 12 (6), 1535. doi:10.3390/cancers12061535
Frigola, J., Navarro, A., Carbonell, C., Callejo, A., Iranzo, P., Cedrés, S., et al. (2021). Molecular Profiling of Long‐term Responders to Immune Checkpoint Inhibitors in Advanced Non‐small Cell Lung Cancer. Mol. Oncol. 15 (4), 887–900. doi:10.1002/1878-0261.12891
Hellmann, M. D., Nathanson, T., Rizvi, H., Creelan, B. C., Sanchez-Vega, F., Ahuja, A., et al. (2018). Genomic Features of Response to Combination Immunotherapy in Patients with Advanced Non-small-Cell Lung Cancer. Cancer Cell 33 (5), 843–852. doi:10.1016/j.ccell.2018.03.018
Hou, J., and Yao, C. (2021). Potential Prognostic Biomarkers of Lung Adenocarcinoma Based on Bioinformatic Analysis. Biomed. Res. Int. 2021, 1–14. doi:10.1155/2021/8859996
Hu, W. L., Jin, L., Xu, A., Wang, Y. F., Thorne, R. F., Zhang, X. D., et al. (2018). GUARDIN Is a P53-Responsive Long Non-coding RNA that Is Essential for Genomic Stability. Nat. Cel Biol. 20 (4), 492–502. doi:10.1038/s41556-018-0066-7
Jachimowicz, R. D., Beleggia, F., Isensee, J., Velpula, B. B., Goergens, J., Bustos, M. A., et al. (2019). UBQLN4 Represses Homologous Recombination and Is Overexpressed in Aggressive Tumors. Cell 176 (3), 505–519. doi:10.1016/j.cell.2018.11.024
Jin, D., Song, Y., Chen, Y., and Zhang, P. (2020). Identification of a Seven-lncRNA Immune Risk Signature and Construction of a Predictive Nomogram for Lung Adenocarcinoma. Biomed. Res. Int. 2020, 1–17. doi:10.1155/2020/7929132
Kandoth, C., McLellan, M. D., Vandin, F., Ye, K., Niu, B., Lu, C., et al. (2013). Mutational Landscape and Significance across 12 Major Cancer Types. Nature 502 (7471), 333–339. doi:10.1038/nature12634
Lee, S., Lee, S., Roh, H.-S., Song, S.-S., Ryoo, R., Pang, C., et al. (2018). Cytotoxic Constituents from the Sclerotia of Poria Cocos against Human Lung Adenocarcinoma Cells by Inducing Mitochondrial Apoptosis. Cells 7 (9), 116. doi:10.3390/cells7090116
Lengauer, C., Kinzler, K. W., and Vogelstein, B. (1998). Genetic Instabilities in Human Cancers. Nature 396 (6712), 643–649. doi:10.1038/25292
Li, R., Han, K., Xu, D., Chen, X., Lan, S., Liao, Y., et al. (2020). A Seven-Long Non-coding RNA Signature Improves Prognosis Prediction of Lung Adenocarcinoma: An Integrated Competing Endogenous RNA Network Analysis. Front. Genet. 11, 625977. doi:10.3389/fgene.2020.625977
Liu, H. (2016). Linking lncRNA to Genomic Stability. Sci. China Life Sci. 59 (3), 328–329. doi:10.1007/s11427-016-5009-6
Liu, Y., Wu, L., Ao, H., Zhao, M., Leng, X., Liu, M., et al. (2019). Prognostic Implications of Autophagy-Associated Gene Signatures in Non-small Cell Lung Cancer. Aging 11 (23), 11440–11462. doi:10.18632/aging.102544
Macheleidt, I. F., Dalvi, P. S., Lim, S. Y., Meemboor, S., Meder, L., Käsgen, O., et al. (2018). Preclinical Studies Reveal that LSD 1 Inhibition Results in Tumor Growth Arrest in Lung Adenocarcinoma Independently of Driver Mutations. Mol. Oncol. 12 (11), 1965–1979. doi:10.1002/1878-0261.12382
Mattick, J. S., and Rinn, J. L. (2015). Discovery and Annotation of Long Noncoding RNAs. Nat. Struct. Mol. Biol. 22 (1), 5–7. doi:10.1038/nsmb.2942
McLaughlin, J., Han, G., Schalper, K. A., Carvajal-Hausdorf, D., Pelekanou, V., Rehman, J., et al. (2016). Quantitative Assessment of the Heterogeneity of PD-L1 Expression in Non-small-Cell Lung Cancer. JAMA Oncol. 2 (1), 46–54. doi:10.1001/jamaoncol.2015.3638
Nair, L., Chung, H., and Basu, U. (2020). Regulation of Long Non-coding RNAs and Genome Dynamics by the RNA Surveillance Machinery. Nat. Rev. Mol. Cel Biol. 21 (3), 123–136. doi:10.1038/s41580-019-0209-0
Negrini, S., Gorgoulis, V. G., and Halazonetis, T. D. (2010). Genomic Instability - an Evolving Hallmark of Cancer. Nat. Rev. Mol. Cel Biol. 11 (3), 220–228. doi:10.1038/nrm2858
Pardoll, D. M. (2012). The Blockade of Immune Checkpoints in Cancer Immunotherapy. Nat. Rev. Cancer 12 (4), 252–264. doi:10.1038/nrc3239
Parry, E. M., Gable, D. L., Stanley, S. E., Khalil, S. E., Antonescu, V., Florea, L., et al. (2017). Germline Mutations in DNA Repair Genes in Lung Adenocarcinoma. J. Thorac. Oncol. 12 (11), 1673–1678. doi:10.1016/j.jtho.2017.08.011
Peters, S., Gettinger, S., Johnson, M. L., Jänne, P. A., Garassino, M. C., Christoph, D., et al. (2017). Phase II Trial of Atezolizumab as First-Line or Subsequent Therapy for Patients with Programmed Death-Ligand 1-Selected Advanced Non-small-Cell Lung Cancer (BIRCH). JCO 35 (24), 2781–2789. doi:10.1200/JCO.2016.71.9476
Siegel, R. L., Miller, K. D., and Jemal, A. (2019). Cancer Statistics, 2019. CA A. Cancer J. Clin. 69 (1), 7–34. doi:10.3322/caac.21551
Soca-Chafre, G., Montiel-Dávalos, A., Rosa-Velázquez, I. A. D. L., Caro-Sánchez, C. H. S., Peña-Nieves, A., and Arrieta, O. (2019). Multiple Molecular Targets Associated with Genomic Instability in Lung Cancer. Int. J. Genomics 2019, 1–8. doi:10.1155/2019/9584504
Subbiah, V., Solit, D. B., Chan, T. A., and Kurzrock, R. (2020). The FDA Approval of Pembrolizumab for Adult and Pediatric Patients with Tumor Mutational burden (TMB) ≥10: a Decision Centered on Empowering Patients and Their Physicians. Ann. Oncol. 31 (9), 1115–1118. doi:10.1016/j.annonc.2020.07.002
Tam, A. S., Sihota, T. S., Milbury, K. L., Zhang, A., Mathew, V., and Stirling, P. C. (2019). Selective Defects in Gene Expression Control Genome Instability in Yeast Splicing Mutants. MBoC 30 (2), 191–200. doi:10.1091/mbc.E18-07-0439
Terasawa, M., Shinohara, A., and Shinohara, M. (2014). Double‐strand Break Repair‐adox: Restoration of Suppressed Double‐strand Break Repair during Mitosis Induces Genomic Instability. Cancer Sci. 105 (12), 1519–1525. doi:10.1111/cas.12551
Wang, J., Gao, J., Chen, Q., Zou, W., Yang, F., Wei, C., et al. (2020). LncRNA LINC01116 Contributes to Cisplatin Resistance in Lung Adenocarcinoma. Ott 13, 9333–9347. doi:10.2147/OTT.S244879
Wei, C.-H., Wu, G., Cai, Q., Gao, X.-C., Tong, F., Zhou, R., et al. (2017). RETRACTED ARTICLE: MicroRNA-330-3p Promotes Cell Invasion and Metastasis in Non-small Cell Lung Cancer through GRIA3 by Activating MAPK/ERK Signaling Pathway. J. Hematol. Oncol. 10 (1), 125. doi:10.1186/s13045-017-0493-0
Wen, S., Hou, Y., Fu, L., Xi, L., Yang, D., Zhao, M., et al. (2019). Cancer-associated Fibroblast (CAF)-derived IL32 Promotes Breast Cancer Cell Invasion and Metastasis via Integrin β3-p38 MAPK Signalling. Cancer Lett. 442, 320–332. doi:10.1016/j.canlet.2018.10.015
Williams, A. B., and Schumacher, B. (2016). p53 in the DNA-Damage-Repair Process. Cold Spring Harb Perspect. Med. 6 (5), a026070. doi:10.1101/cshperspect.a026070
Xiao, G., Wang, P., Zheng, X., Liu, D., and Sun, X. (2019). FAM83A-AS1 Promotes Lung Adenocarcinoma Cell Migration and Invasion by Targeting miR-150-5p and Modifying MMP14. Cell Cycle 18 (21), 2972–2985. doi:10.1080/15384101.2019.1664225
Zeng, L., Lyu, X., Yuan, J., Wang, W., Zhao, N., Liu, B., et al. (2020). Long Non-coding RNA LINC01116 Is Overexpressed in Lung Adenocarcinoma and Promotes Tumor Proliferation and Metastasis. Am. J. Transl. Res. 12 (8), 4302–4313.
Zhou, M., Shao, W., Dai, H., and Zhu, X. (2020). A Robust Signature Based on Autophagy-Associated LncRNAs for Predicting Prognosis in Lung Adenocarcinoma. Biomed. Res. Int. 2020, 1–13. doi:10.1155/2020/3858373
Keywords: lung cancer, lung adenocarcinoma, long non-coding RNA, genomic instability, prognosis
Citation: Tu G, Peng W, Cai Q, Zhao Z, Peng X, He B, Zhang P, Shi S and Wang X (2021) A Novel Model Based on Genomic Instability-Associated Long Non-Coding RNAs for Predicting Prognosis and Response to Immunotherapy in Patients With Lung Adenocarcinoma. Front. Genet. 12:720013. doi: 10.3389/fgene.2021.720013
Received: 03 June 2021; Accepted: 04 October 2021;
Published: 29 October 2021.
Edited by:
Vince Kornél Grolmusz, National Institute of Oncology (NIO), HungaryReviewed by:
Tin Lap Lee, The Chinese University of Hong Kong, ChinaJuanni Li, Central South University, China
Copyright © 2021 Tu, Peng, Cai, Zhao, Peng, He, Zhang, Shi and Wang. 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: Xiang Wang, wangxiang@csu.edu.cn