Immune-Stromal Score Signature: Novel Prognostic Tool of the Tumor Microenvironment in Lung Adenocarcinoma

Background: Immune and stromal cells in the tumor microenvironment (TME) significantly contribute to the prognosis of lung adenocarcinoma; however, the TME-related immune prognostic signature is unknown. The aim of this study was to develop a novel immune prognostic model of the TME in lung adenocarcinoma. Methods: First, the immune and stromal scores among lung adenocarcinoma patients were determined using the ESTIMATE algorithm in accordance with The Cancer Genome Atlas (TCGA) database. Differentially expressed immune-related genes (IRGs) between high and low immune/stromal score groups were analyzed, and a univariate Cox regression analysis was performed to identify IRGs significantly correlated with overall survival (OS) among patients with lung adenocarcinoma. Furthermore, a least absolute shrinkage and selection operator (LASSO) regression analysis was performed to generate TME-related immune prognostic signatures. Gene set enrichment analysis was performed to analyze the mechanisms underlying these immune prognostic signatures. Finally, the functions of hub IRGs were further analyzed to delineate the potential prognostic mechanisms in comprehensive TCGA datasets. Results: In total, 702 intersecting differentially expressed IRGs (589 upregulated and 113 downregulated) were screened. Univariate Cox regression analysis revealed that 58 significant differentially expressed IRGs were correlated with patient prognosis in the training cohort, of which three IRGs (CLEC17A, INHA, and XIRP1) were identified through LASSO regression analysis. A robust prognostic model was generated on the basis of this three-IRG signature. Furthermore, functional enrichment analysis of the high-risk-score group was performed primarily on the basis of metabolic pathways, whereas analysis of the low-risk-score group was performed primarily on the basis of immunoregulation and immune cell activation. Finally, hub IRGs CLEC17A, INHA, and XIRP1 were considered novel prognostic biomarkers for lung adenocarcinoma. These hub genes had different mutation frequencies and forms in lung adenocarcinoma and participated in different signaling pathways. More importantly, these hub genes were significantly correlated with the infiltration of CD4+ T cells, CD8+ T cells, macrophages, B cells, and neutrophils. Conclusions: The robust novel TME-related immune prognostic signature effectively predicted the prognosis of patients with lung adenocarcinoma. Further studies are required to further elucidate the regulatory mechanisms of these hub IRGs in the TME and to develop new treatment strategies.


INTRODUCTION
Lung cancer is still the leading disease worldwide in terms of the threat to human life and health (1,2), and lung adenocarcinoma is the most common pathological subtype. Studies in the past decade have reported that tyrosine kinase inhibitors (TKIs) targeting epidermal growth factor receptor (EGFR), anaplastic lymphoma kinase (ALK), and ROS proto-oncogene 1 (ROS1) are potential therapeutic targets for lung adenocarcinoma, upon genotyping (3)(4)(5). Molecular-targeted therapy based on these sensitive targets has considerably enhanced overall survival (OS) among patients with lung adenocarcinoma; however, this therapy is not suitable for all patients with lung adenocarcinoma. Furthermore, drug resistance is common among patients receiving molecular-targeted therapy, and their prognosis is poor (6,7). Nonetheless, numerous studies have led to the advancement of immunotherapy for several cancers, including lung adenocarcinoma. Immunotherapy is different from targeted therapy; it has more durable clinical benefits. Furthermore, some antibodies used for immunotherapy have been successfully approved as first-and second-line treatments for advanced lung adenocarcinoma (8). In particular, the immune system reportedly plays an important role in the pathogenesis and prognosis of lung adenocarcinoma (9). Therefore, it is essential to understand the immune prognostic signature of lung adenocarcinoma.
Previous studies have investigated the prognostic role of immune-related genes (IRGs) in lung adenocarcinoma from the ImmPort database (10,11); however, this database contains published data on IRGs, thus potentially not accounting for all IRGs. Moreover, these studies have reported no correlation between prognostic factors and OS among certain subgroups of patients with lung adenocarcinoma, indicating that this association is largely unknown. One of the important reasons may be the complex prognostic behavior of tumors; furthermore, when considering the characteristics of IRGs directly associated with tumors, it is also important to focus on the tumor microenvironment (TME) (12). The TME is closely associated with tumorigenesis and patient prognosis (13,14). Moreover, accumulating evidence indicates that tumor-infiltrating immune cells and stromal cells (15,16), as the primary nontumor components of the TME, play a significant role in lung cancer prognosis. These findings highlight the importance of understanding the association between the TME and lung adenocarcinoma prognosis. The development of a prognostic Abbreviations: IRGs, immune-related genes; TME, tumor microenvironment; TCGA, The Cancer Genome Atlas data; DAVID, Database for Annotation, Visualization and Integrated Discovery. model of IRGs based on the TME might provide novel insights into the generation of a more accurate prognostic system. Accurate management and appropriate personalized therapies for lung adenocarcinoma are required in accordance with prognostic stratification. Moreover, an enhanced understanding of IRGs involved in the TME would help elucidate their regulatory mechanisms in the TME and develop new treatment strategies. With advancements in machine learning, the ESTIMATE algorithm has been used to investigate IRGs in the TME, based on the immune and stromal scores of the TME, and to generate a TME-related immune prognostic model (17). Moreover, this algorithm can effectively predict the prognosis of patients with various cancers. Accordingly, in this study, we determined the immune and stromal scores of tumors using the ESTIMATE algorithm and developed a novel TME-related immune prognostic model of lung adenocarcinoma.

Acquisition of TCGA Data
Normalization of RNA sequence data, in terms of level 3 fragments per kilobase of exon per million fragments mapped (FPKM) reads, was performed for 594 samples obtained from The Cancer Genome Atlas (TCGA) database, including 535 adenocarcinoma and 59 normal lung samples, before December 15, 2019. Thereafter, the Ensemble IDs were converted to gene symbols in accordance with human gene annotations. Furthermore, clinical data of lung adenocarcinoma patients were obtained and merged into a single matrix for subsequent analysis. Patients with an incomplete follow-up duration or recorded date of death of any cause were excluded. Finally, 494 lung adenocarcinoma patients with expression profiles and clinical data were included.
Immune Score and Stromal Score in the TME Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) is a tool for predicting and estimating infiltrating immune and stromal cells in tumor tissues based on gene expression profiles. Herein, the ESTIMATE algorithm was used to analyze the characteristics of specific gene expression in immune and stromal cells for each tumor sample to predict their immune and stromal scores. Thereafter, the immune and stromal scores were analyzed using the estimate package in R software.

Screening of Differentially Expressed IRGs
Based on the median immune and stromal scores, patients with lung adenocarcinoma were divided into two groups: high and low immune/stromal score groups. Significant differences in OS between the high and low immune/stromal score groups were analyzed. On the basis of the significant differences in patient prognosis, differentially expressed IRGs between the two groups were assessed using the Limma package in R software. Finally, intersecting differentially expressed IRGs in both groups were considered for further analysis. A log(fold change) of >2 and an adjusted p-value of <0.05 were considered cutoffs. Heat maps and Venn diagrams were generated using R.

GO and KEGG Pathway Enrichment Analyses of Differentially Expressed IRGs
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed to understand gene functional annotation and functional enrichment, respectively. Common differentially expressed IRGs of GO and KEGG annotation were performed using Database for Annotation, Visualization, and Integrated Discovery (DAVID), an online database (https://david.ncifcrf. gov/), to predict functional domains and their biological implications. Fisher's exact test was performed to analyze pathways, diseases, and functions. A p-value of <0.05 indicated the significance of GO terms and KEGG pathway enrichment in genes herein. Furthermore, the top 10 GO terms and KEGG pathway enrichment results were mapped using Hmisc and ggplot2 in R software.

Construction of Prognostic Model for Lung Adenocarcinoma in the Training Set
A total of 494 patients with lung adenocarcinoma in the TCGA dataset were randomly divided into a training set and testing set at a ratio of 7:3 (training cohort, 346 patients; testing cohort, 148 patients).
First, univariate Cox proportional hazard regression analysis was performed to screen out prognostic IRGs in the training cohort (iteration = 1,000), with p < 0.05 indicating statistical significance. Second, the key IRGs were further selected from among significant prognostic IRGs on univariate analysis, through least absolute shrinkage and selection operator (LASSO) regression, a powerful tool for developing refined prognostic models, fitting generalized linear models, selecting variables, and regularizing complexity, using R software. Key IRGs were subjected to multivariate Cox regression analysis. Finally, the risk score formula was developed in accordance with the key IRGs identified through LASSO analysis.

Evaluation of the Prognostic Model in the Training Set
After the expression value of each specific gene was included, the risk score formula for each patient was weighted by its estimated regression coefficient on LASSO regression analysis. On the basis of the best separation of risk score, patients were divided into high-risk-score and low-risk-score groups. Survival differences between the two groups were assessed by Kaplan-Meier survival curves using log-rank tests. ROC curves were used to assess the accuracy of model prediction. Furthermore, LASSO regression analysis was performed to examine the role of the risk score in predicting clinical outcomes. Furthermore, the association between risk score and clinical stage was analyzed.
To further determine whether the independent prognostic model could be used as an independent prognostic factor, univariate and multivariate Cox regression analyses were performed to analyze the predictive value of age, sex, stage, TNM stage, and predictive model. In the univariate analysis, the correlation between some independent variables and dependent variables was considered, and some correlations might be masked by the influence of confounding factors. To avoid omitting important predictors of prognosis, the threshold in univariate analysis was relaxed to p < 0.1, while the p-value in multivariate analysis was still 0.05.

GSEA Analysis of Differences in Pathway Enrichment in the Training Set
To investigate the differences in the putative mechanism between the high-risk-score and low-risk-score groups, gene set enrichment analysis (GSEA) was performed to comprehensively analyze the differences in function enrichment. GSEA is a computational method of determining whether an a priori defined set of genes is significantly different between two biological states. The number of permutations was set to 1,000, and the permutation type was set to phenotype. In this study, all genes in the training set were sequenced according to the degree of differential expression in the high-risk-score group and the low-risk-score group. GSEA was used to comprehensively analyze differences in gene pathway enrichment between the two groups.

Validation of the Prognostic Model in the Testing Cohort
The prognostic model was further validated for the testing cohort (n = 148).
Similarly, the risk score of each patient was weighted on the basis of the risk score. Thereafter, based on the best separation of the risk score, patients in the testing cohort were divided into high-risk-score and low-risk-score groups. A Kaplan-Meier survival curve and ROC curve analysis were performed in the testing however.

Functional Analysis of IRG Signatures in the Model
To further analyze the mutation characteristics and the putative functional mechanisms of these hub IRGs in lung adenocarcinoma, gene expression profiles of patients with lung adenocarcinoma were imported from the following datasets:  combined study of five datasets including 1,825 patients were included in this study. GSCALite is an online cancer genomic analysis tool that integrates cancer genomics data for 33 cancer types from the TCGA and normal tissue data from GTEx (http://bioinfo.life. hust.edu.cn/web/GSCALite/), enabling gene set pathway analysis in data analysis. In this study, GSCALite was used to analyze the pathway of hub genes.
The TIMER database is a comprehensive tool for analyzing the immune cell infiltrates in tumors. The abundances of six immune infiltrates (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells) were estimated using the TIMER algorithm. In this study, the TIMER database was further applied to analyze the correlation between hub genes and immune cells. The correlation of hub gene expression with immune infiltration level was visualized in lung adenocarcinoma using the Gene module.
The scatterplots were generated and displayed after hub genes and cancer type were submitted successfully, showing the purity-corrected partial Spearman's correlation and statistical significance.

Statistical Analysis
All statistical analyses were conducted with the R language (version 3.6.1). All statistical tests were bilateral, and p < 0.05 was statistically significant.

Study Design and Workflow Overview
A total of 594 samples were obtained, including 535 adenocarcinoma and 59 normal lung samples. Data for a total of 494 lung adenocarcinoma patients with clinical information were retrieved ( Table 1). The workflow for constructing and Immune Score and Stromal Score in the TME The immune and stromal scores were analyzed using the ESTIMATE algorithm. The immune and stromal scores of lung adenocarcinoma are shown in Supplementary Table 1.
On the basis of the median value of immune and stromal scores, patients with lung adenocarcinoma were divided into two groups: the high immune/stromal score group and the low score group. These results show that the high immune score group had better OS than the low immune score group for lung adenocarcinoma. In terms of stromal score, although patients with high stromal score had better prognosis than those with low stromal score, the difference was not statistically significant (Figures 2A-C).

Screening of Differentially Expressed IRGs
According to the criteria of a log(fold change) of >2 and an adjusted p-value of <0.05, our results showed that a total of 3,034 genes with significant differentially expressed IRGs were screened, including 2,521 upregulated IRGs and 513 downregulated IRGs. Among them, 1,394 differentially expressed IRGs (1,092 upregulated IRGs and 302 downregulated IRGs) were included in the immune score group (Supplementary Table 2 Table 4).

GO Terms and KEGG Pathway Enrichment Analysis of Differentially Expressed IRGs
GO terms are divided into three parts: biological processes, cellular components, and molecular functions. GO analysis    Table 2). The KEGG pathway enrichment results in the upregulated IRGs were mainly involved in cytokine-cytokine receptor interactions, chemokine signaling pathways, and cell adhesion molecules (CAMs). However, the KEGG pathway enrichment results of downregulated IRGs were mainly involved in arachidonic acid metabolism, metabolic pathways, and tyrosine metabolism (Figures 3D,F; Table 3).

Evaluation of the Prognostic Model in the Training Cohort
LASSO regression analysis was performed to construct and evaluate the prognostic model (Figures 4A,B). Patients were divided into high-and low-risk-score groups in accordance with the best separation of risk scores. The high-risk-score group had significantly worse OS than the low-risk-score group (p < 0.0001; Figure 4C). The area under the ROC curve for predicting the 1-, 3-, and 5-year survival of lung adenocarcinoma was 0.699, 0.631, and 0.669, respectively ( Figure 4D). Additionally, the risk curve indicated that the high-risk-score group had a higher mortality and worse prognosis than the low-risk-score group (cutoff value: 0.889; Figure 4E). Further analysis of the relationship between risk score and pathological stage revealed that patients with early-stage lung adenocarcinoma (stages 1 and 2) scored lower than those with advanced stage lung adenocarcinoma (p = 0.043; Figure 4F).

GSEA of the Mechanism Underlying the Prognostic Differences Between the Two Groups
In this study, the possible molecular mechanisms of the prognosis difference between the two groups of patients were analyzed by GSEA analysis. The results showed that the GO and pathway enrichment in the high-risk-score group was mainly involved in metabolism-related pathways (Figures 5A,C). However, GO and pathway enrichment in the low-risk-score group was primarily focused on immunoregulation and immune cell activation (Figures 5B,D). The detailed GSEA results are described in Table 5.

Validation of the Prediction Model in the Testing Cohort
The Kaplan-Meier results showed that the high-risk group had worse OS than the low-risk group (p < 0.0001) (Figure 6A). The area under the ROC curve for predicting the 1-, 3-, and 5-year survival of lung adenocarcinoma was 0.725, 0.712, and 0.660, respectively ( Figure 6B). Additionally, risk curve revealed that the high-risk-score group had a worse prognosis than the lowrisk-score group (Figure 6C). These results were consistent with the results of the training set.

The Mechanism of Action of Hub IRGs in the TCGA Database
To further analyze the potential function of hub IRGs, our results were verified using the TCGA and GTEx databases. First, the mutation characteristics of these hub IRGs were analyzed among patients with lung adenocarcinoma. The mutation rates of these hub IRGs in patients with lung adenocarcinoma were 0.8, 1, and 2.6% ( Figure 7A). Moreover, each hub gene had different mutation forms, including mutation, deletion, and amplification, in lung adenocarcinoma. For example, the mutation form of CLEC17A was mainly amplification, the mutation form of INHA was mainly amplification and missense mutation, while the mutation form of XIRP1 was mainly deep deletion and missense mutation ( Figure 7B).
In addition, analysis of pathways in the GSCALite database revealed that CLEC17A is primarily involved in the activation of the epithelial-mesenchymal transition (EMT) and RAS pathways and cell cycle inhibition. INHA is primarily involved in the activation of the mTOR pathway and inhibition of the apoptosis pathway. XIRP1 was mainly involved in the activation of the apoptosis and the EMT pathways and inhibition of the DNA damage and the PI3K pathways (Figures 7C,D).
Finally, the TIMER database was used to analyze the correlation between hub IRGs and immune cells. The results showed that these key genes were significantly correlated with the infiltration of CD4+ T cells, CD8+ T cells, macrophages, B cells, and neutrophils. Assuming that a correlation coefficient >0.3 was considered a strong correlation, further analysis showed that CLEC17A was positively correlated with the infiltration of B cells and CD4+ T cells. However, INHA was negatively correlated with the infiltration of CD8+ T and dendritic cells. However, there was no strong correlation between XIRP1 and the infiltration of immune cells (Figure 7E). Moreover, the relationship between copy number variation (CNV) of hub IRGs and immune cell infiltration was further analyzed. The results showed that there were significant differences between the CNV of these hub IRGs and immune cell infiltration. Armlevel deletion of the CLEC17A gene was closely related to the infiltration of B cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells. Arm-level gain of the INHA gene was closely related to the infiltration of CD4+ T cells. Arm-level deletion of the XIRP1 gene was closely related to the infiltration of CD4+ T cells and macrophages (Figure 7F).

DISCUSSION
Tumor-infiltrating immune cells in the TME significantly contribute to the prognosis of lung adenocarcinoma. Therefore, it is essential to develop a TME-related immune prognostic model for appropriate clinical management of lung adenocarcinoma. Accordingly, the aim of this study was to develop a novel TMErelated immune prognostic model for lung adenocarcinoma.
Although some studies have explored the prognostic value of TME-related IRGs in lung adenocarcinoma, certain key issues of these models remain to be resolved. Yang et al. (18) used the CIBERSORT algorithm to analyze a TME-related prognostic immunity-based model for lung adenocarcinoma. This model was developed to evaluate the relative levels of the 22 immune cell phenotypes, primarily including B cells, T cells, macrophages, dendritic cells, plasma cells, natural killer cells, and mast cells. Moreover, this algorithm is primarily used to evaluate immune cells; however, it cannot be used to evaluate stromal cells in the TME. Yue et al. (19) used the ESTIMATE algorithm to investigate the TME-related immune prognostic characteristics of lung adenocarcinoma. However, they directly enumerated immune and stromal cells from their expression profiles of all genes expressed in lung adenocarcinoma and normal tissues. Moreover, differentially expressed IRGs were analyzed using Wilcoxon correlation analysis between tumor and normal tissue. However, some limitations are associated with the analysis of multi-dimensional tumor gene expression profiles through the Wilcoxon rank-sum test. These findings indicate that these TME-related prognostic immunity-based models have not been adequately evaluated. These issues can be resolved primarily by improving the algorithm of IRGs in the TME and to identify more specific TME-related IRGs for lung adenocarcinoma. Therefore, it is essential to develop a new TME-related immune prognostic model for lung adenocarcinoma.
To address these aforementioned limitations, in the present study, a new method was developed to identify differentially expressed IRGs. First, TME-related differentially expressed IRGs were identified exclusively from tumor samples by evaluating tumor-infiltrating immune cells and stromal cells via the ESTIMATE algorithm; this probably effectively reflected the TME-related IRGs in tumor tissue. Second, differentially expressed IRGs were analyzed on the basis of significant differences in OS between the high-and lowimmune-score groups in terms of lung adenocarcinoma, rather than differences between lung adenocarcinoma and normal tissue using the Wilcoxon rank-sum test. The prognostic model, based on prognosis-associated differentially expressed genes, might more accurately predict the prognosis of lung adenocarcinoma. Finally, intersecting differentially expressed IRGs with significant prognostic characteristics in both immune and stromal scores were used for subsequent analysis. Both immune cells and stromal cells in each tumor sample were assessed, thus better reflecting the characteristics of the TME. Therefore, the TME-related immune prognostic model developed herein was different from those developed previously. In this study, we developed a more robust prognostic model of TME-related IRGs in lung adenocarcinoma.
Furthermore, this study shows that patients with lung adenocarcinoma and high immune scores had a better prognosis than those with low immune scores, which might be due to the involvement of upregulated IRGs in immune cell infiltration factors, such as cytokines and B cell immune pathways. Clinical studies have also shown that lung cancer patients with high immune infiltration of helper T cells have a better prognosis than those with low infiltration (20,21). These findings were consistent with our results. Previous studies have suggested that IL-2 is involved in antitumor T cell infiltration, increasing the efficacy of immunotherapy (22). IL-33 also promotes myeloidderived suppressor cells (MDSCs) and interferes with CD8+ T and natural killer (NK) cell infiltration (23). These studies have suggested that certain cytokines are involved in antitumor immune pathways, potentially elucidating the mechanisms associated with prognosis. Moreover, GSEA was performed to further investigate the potential mechanism underlying the differences in prognosis between the two groups. The present results indicate that the immunoregulation and immune cell activation pathways are potentially associated with a better prognosis. The underlying putative mechanism potentially involves the enrichment of B and T cell immune pathways. Furthermore, the infiltration of these immune cells is associated with an enhanced prognosis among patients with lung adenocarcinoma. Our results are concurrent with those of the aforementioned studies. Furthermore, this study described the functional prediction of potential hub IRGs. An enhanced understanding of these potential hub genes is essential to elucidate their mechanisms of action in the TME in lung adenocarcinoma. The present results suggest that although these genes were prognosis-related IRGs, they harbored different mutations involved in different pathways in lung adenocarcinoma, indicating their potential involvement in different immunoregulatory pathways in lung adenocarcinoma. Further analysis of the function of these genes revealed that these hub genes and their CNVs were different. Moreover, the association between these hub genes and the infiltration of B cells, CD4+ T cells, CD8+ T cells, neutrophils, and other immune cells was also different. These results indicate that different CNVs of these hub genes warrant further differentiation to better understand the association between hub genes and immune cell infiltration. Based on the aforementioned results, our results indicate that CLEC17A, INHA, and XIRP1 are potential novel biomarkers for the prognosis of lung adenocarcinoma.
CLEC17A is a human lectin found in lymph node B cells and is involved in a variety of biological processes, including cell adhesion, intercellular interactions, and pathogen recognition (24). Previous studies have shown that CLEC17A is related to the B cell receptor signaling pathway and plays an important role in the pathogenesis of chronic lymphocytic leukemia (25). The present results further indicate that CLEC17A is associated with immune cell infiltration in lung adenocarcinoma, concurrent with previous reports. XIRP1 is a striated muscle protein and belongs to the Xin actin-binding repeat-containing protein (XIRP) family. Previous studies have shown that the XIRP1 gene is related to hypertension and nervous system development (26,27). The function of the gene has not been reported in tumors. However, our study showed that it was not only related to the TME in lung adenocarcinoma but is also related to the prognosis of lung adenocarcinoma. INHA encodes a member of the transforming growth factor-beta (TGF-beta) superfamily of proteins. The function of the gene has not been reported in lung cancer. Our results showed that INHA was a marker of poor prognosis in lung adenocarcinoma. The possible mechanism was that INHA was involved in tumor angiogenesis, leading to tumor metastasis and poor prognosis (28). Further studies are required to elucidate the roles of these hub genes in the TME in the pathogenesis of lung adenocarcinoma.
This study also had some limitations. First, this study only mined data in the TCGA database and did not combine GEO database analysis. However, in our study, patient data were segregated into training and testing cohorts. In addition, the results were verified and analyzed in comprehensive TCGA and GSCALite datasets. Second, the function of the hub gene in our study was analyzed based on the TCGA database, and the validation function of the hub gene needs to be further confirmed by basic experiments. Constructing an immunerelated prognosis model was the focus of our research; hence, there was no basic experiment on hub prognostic genes.
Third, our study only analyzed the correlation of differentially expressed IRGs with immune cell infiltration, thus lacking the correlation analysis of the expression of PDL1 and tumor mutational burden.

CONCLUSIONS
The robust TME-related immune prognostic model developed herein effectively predicted the prognosis of patients with lung adenocarcinoma, thus potentially guiding personalized treatment of lung adenocarcinoma in accordance with prognostic stratification. Further studies are required to elucidate the regulatory mechanisms of these IRGs in the TME and develop new treatment strategies.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://portal.gdc.cancer.gov/.

AUTHOR CONTRIBUTIONS
XQ designed the study, analyzed the data, and drafted the paper. CQ critically revised it for important intellectual content. BQ and XK assisted in data acquisition and analysis. YH and WH revised the manuscript. All authors read and approved the final version of the manuscript.

FUNDING
The authors received no funding for this work.