Skip to main content

ORIGINAL RESEARCH article

Front. Immunol., 04 April 2023
Sec. Cancer Immunity and Immunotherapy

Construction of a novel anoikis-related prognostic model and analysis of its correlation with infiltration of immune cells in neuroblastoma

Ji Chen&#x;Ji Chen1†Mengjiao Sun&#x;Mengjiao Sun2†Chuqin Chen&#x;Chuqin Chen2†Meiyun Kang&#x;Meiyun Kang2†Bo QianBo Qian3Jing SunJing Sun2Xiaopeng MaXiaopeng Ma2Jianfeng ZhouJianfeng Zhou1Lei Huang*Lei Huang1*Bin Jiang*Bin Jiang1*Yongjun Fang*Yongjun Fang2*
  • 1Department of General Surgery, Children’s Hospital of Nanjing Medical University, Nanjing, China
  • 2Department of Hematology and Oncology, Children’s Hospital of Nanjing Medical University, Nanjing, China
  • 3Department of Cardiothoracic Surgery, Children’s Hospital of Nanjing Medical University, Nanjing, China

Background: Anoikis resistance (AR) plays an important role in the process of metastasis, which is an important factor affecting the risk stage of neuroblastoma (NB). This study aims to construct an anoikis-related prognostic model and analyze the characteristics of hub genes, important pathways and tumor microenvironment of anoikis-related subtypes of NB, so as to provide help for the clinical diagnosis, treatment and research of NB.

Methods: We combined transcriptome data of GSE49710 and E-MTAB-8248, screened anoikis-related genes (Args) closely related to the prognosis of NB by univariate cox regression analysis, and divided the samples into anoikis-related subtypes by consistent cluster analysis. WGCNA was used to screen hub genes, GSVA and GSEA were used to analyze the differentially enriched pathways between anoikis-related subtypes. We analyzed the infiltration levels of immune cells between different groups by SsGSEA and CIBERSORT. Lasso and multivariate regression analyses were used to construct a prognostic model. Finally, we analyzed drug sensitivity through the GDSC database.

Results: 721 cases and 283 Args were included in this study. All samples were grouped into two subtypes with different prognoses. The analyses of WGCNA, GSVA and GSEA suggested the existence of differentially expressed hub genes and important pathways in the two subtypes. We further constructed an anoikis-related prognostic model, in which 15 Args participated. This model had more advantages in evaluating the prognoses of NB than other commonly used clinical indicators. The infiltration levels of 9 immune cells were significantly different between different risk groups, and 13 Args involved in the model construction were correlated with the infiltration levels of immune cells. There was a relationship between the infiltration levels of 6 immune cells and riskscores. Finally, we screened 15 drugs with more obvious effects on NB in high-risk group.

Conclusion: There are two anoikis-related subtypes with different prognoses in the population of NB. The anoikis-related prognostic model constructed in this study can accurately predict the prognoses of children with NB, and has a good guiding significance for clinical diagnosis, treatment and research of NB.

Introduction

NB is a kind of malignant solid tumor originating from neural crest stem cells, which is highly heterogeneous and most prone to occur outside the central nervous system in children (1). Surgery combined with chemotherapy can achieve a good prognosis for low-risk (2) and intermediate-risk (3) NB, while the prognoses of patients with high-risk NB are still poor even with a combination of surgery (4, 5), chemotherapy (68), radiotherapy (9), bone marrow transplantation (10) and other treatments. As an important therapeutic method for malignant tumors, immunotherapy has been widely used in the treatment of relapsed and refractory NB. With the application of immunotherapy methods such as GD2 monoclonal antibody (1113) and CAR-T (1416), the prognoses of children with relapsed and refractory NB have been significantly improved. However, how to improve the prognosis of high-risk NB with metastasis is still a problem that needs to be solved. Previous studies have shown that the prognosis of NB with bone metastasis or bone marrow metastasis is generally poor (17). Therefore, exploring the mechanism of NB metastasis will help us better understand the biological characteristics of NB, and provide the basis for molecular targeted therapy of high-risk NB.

Metastasis is a very complex biological process. In general, cells attach to the extracellular matrix (ECM) through signal transduction and carry out cell division and differentiation in specific areas (18, 19). Apoptosis occurs when cells leave the ECM, a process known as anoikis (20). Subsequently, tumor cells can invade the vasculature and colonize in distant organs in the absence of ECM and intercellular signaling. This process is known as AR. Anoikis was first found in epithelial cells and endothelial cells, and then AR was found to be an important pathway for tumor metastasis (2124). A large number of studies have shown that the tumor microenvironment (TME) is closely related to the metastasis of malignant tumors (25, 26), and immune cells play important roles in the TME (2730). Therefore, we hypothesize that there might be some relationships between AR and the changes of TME.

At present, the specific mechanism of anoikis in NB has not been elucidated (31, 32), and the relationship between AR and immune cell infiltration in NB is not clear. In this study, two anoikis-related subtypes of NB were screened, and the mechanisms of anoikis in NB were investigated by analyzing hub genes, important pathways and immune cell infiltration of anoikis-related subtypes based on transcriptome data from GEO and ArrayExpress databases. Subsequently, we constructed a reliable anoikis-related prognostic model of NB, and confirmed that this model had a good guiding significance in the evaluation of NB prognosis and the prediction of drug sensitivity. Figure 1 illustrates the workflow of this research.

FIGURE 1
www.frontiersin.org

Figure 1 Flow chart of this study.

Methods

Data acquisition

Transcriptome data from GSE49710 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE49710) and E-MTAB-8248 cohort (https://www.ebi.ac.uk/arrayexpress/experiments/E-MTAB-8248/) in GEO and ArrayExpress databases was included in this study. The two cohorts included 498 and 223 patients respectively, and the main clinical information included age, sex, COG risk, MYCN status, INSS status, survival time, and survival status. The transcriptome data of 39 NB cell lines in GSE89413 cohort (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE89413) was also obtained. After data correction and probe annotation, we finally got the transcriptome matrix data with row names of gene and column names of patient number/cell number. In order to obtain more accurate anoikis-related subtypes and build an anoikis-related prognostic model, we merged the transcriptome matrix data of GSE49710 and E-MTAB-8248, standardized them through the “normalizeBetweenArray ()” function, and finally got the merged matrix data.

Args related to the prognosis of NB

Args were obtained by searching key word “anoikis” in Genecards (https://www.genecards.org/, score>0.4) and Harmonizome (https://maayanlab.cloud/Harmonizome/) databases. Through univariate cox regression analysis of the combined cohort, we screened Args closely related to survival status of children with NB by setting the standard of “HR≠1 and p<0.05”.

Consensus clustering analysis for Args

We conducted consensus clustering analysis in 721 children by analyzing the expression levels of Args closely related to the survival status of NB. We finally obtained the best k-value for clustering by using the “CONSENSUSClusterPlus ()” function and setting the parameters reps=50, pItem=0.8, pFeature=1, clusterAlg=“km”, distance=“euclidean”, seed=123, and used the “prcomp ()” function to perform PCA analysis on all samples. Kaplan-Meier (KM) survival analysis of the anoikis-related subtypes was conducted by the “survdiff ()” function, and Args (|LogFC|>0.5 and p<0.05) differentially expressed in anoikis-related subtypes were screened through the “Lima” package. “Ggboxplot () “ and “ Pheeatmap () “ function were used to draw the boxplot and heatmap of differently expressed Args between different subtypes, and “ Pheeatmap () “ function was also used to observe the distribution of clinical indicators such as age, sex, COG risk, MYCN status and INSS stage between different subtypes.

GSVA and GSEA enrichment analyses

We obtained the hallmark gene sets (c5.go.v7.4.symbols.gmt and c2.cp.Kegg.V7.4.Symbols.gmt) from MSigDB database, conducted GSVA and GSEA enrichment analyses between different subtypes through “GSVA” R package and “GSEA ()” function, and visualize the data through “pheatmap ()” and “gseaplot ()” function respectively.

Analysis of TME

Single-sample gene set enrichment analysis (ssGSEA) was used to analyze the differences of immune cell infiltration between different anoikis-related subtypes. We firstly scored the infiltration levels of 28 kinds of immune cells in each NB sample by using the “GSVA ()” function, then analyzed the differences of immune cell infiltration in different subtypes, and finally visualized the data by using the “ggboxplot ()” function. We used “CIBERSORT ()” function to score immune cell infiltration level for each sample in different risk groups, analyzed the differences in immune cell infiltration level between different risk groups divided according to the prognostic model, and visualized the data by “vioplot ()” function. Finally, we used the “corrplot ()” function to visualize the correlation between immune cells. The “Cor.test()” function was used to analyze the correlation between immune cells and riskscores.

Weighted correlation network analysis

By setting the standard of Height=2000, we used the “hclust ()” function to cluster all the samples and draw a clustering tree after filtering out the outliers. Based on the criterion of approximate scale-free network (R2>0.85), we used the “pick soft threshold ()” function to select an appropriate soft threshold to make the network more in line with the power law distribution, thus reducing the error and making the results more characteristic of biological data. We set 50 as the minimum number of genes in each module and 0.25 as the threshold of cutting height, used the dynamic cutting tree method to identify modules, combined the modules with high similarity and finally obtained multiple important modules. In order to get key modules, we evaluated the correlation between genes and samples by calculating gene significance (GS), defined module significance as the average GS of all genes in the module, calculated the correlation between the modules and anoikis, and selected the module with the highest correlation as the key module. By calculating the GS and module membership (MM) of each gene and setting the interception criteria of GS>0.5 and |MM|>0.80, we finally got the core genes of the key module. By overlapping the prognosis-related Args between two subtypes and the core genes screened by WGCNA, we finally screened the hub genes closely related to anoikis in NB.

Construction and validation of anoikis-related prognostic model

All samples were randomly divided into training group and test group with equal sample size. Lasso regression analysis was performed on all Args related to the survival state in the training group through the “Glmnet ()” function. The minimum regularization parameter lambda (λ) and genes highly associated with anoikis-related subtypes were obtained by cross-validation (alpha=1). Multivariate cox regression analysis was then used to screen the core genes and the weight of each core gene through the “coxph ()” function. Riskscore of each sample was obtained according to the expression level and corresponding weight of core genes. The total sample group, training group and test group were divided into high-risk group and low-risk group according to the median riskscore value respectively. KM survival analysis was performed for the three groups by “survdiff ()” function. “pROC ()” function was used to calculate area under the curve (AUC) in three groups and evaluated the predictive capacity of this model. We used the “coxp ()” function to conduct multivariate cox regression analysis on the clinical information that might affect the prognoses of children with NB, and drew the nomogram through the “replot ()” function. The “calibrate ()” function and “plot ()” function were used to analyze the parameters of calibration curve and plot calibration curve respectively. The “survfit ()” function and “ggsurvplot ()” function were used to analyze the parameters of the cumulative hazard curve and plot the cumulative hazard curve respectively. The parameters of the decision curve were analyzed by “dca ()” function, and the decision curve was plotted by “ggplot ()” function.

Drug sensitivity analysis

We obtained the training set gene expression matrix file (GDSC2_Expr.rds) and drug processing information file (GDSC2_Res.rds) from Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/). The “calcPhenotype ()” function in the “oncoPredict” R package was adopted, and the criteria of removeLowVaryingGenes = 0.2 and minNumSamples = 10 were set for the prediction of drug susceptibility.

Data analysis

All data analyses in this study were completed by R 4.2.1 version, and p < 0.05 was considered statistically significant. The t test and one-way ANOVA were used to analyze the parametric data. The Wilcoxon test and Kruskal-Wallis test were used to analyze the nonparametric data of two and multiple independent samples respectively. Related R packages described above were downloaded from Bioconductor packages or R packages.

Results

Data acquisition and processing

In order to obtain relatively complete clinical data and transcriptome data with a large sample size, we selected GSE49710 and E-MTAB-8248 cohorts from GEO and ArrayExpress databases. The two cohorts included 498 and 223 clinical samples respectively. Both groups contained complete transcriptome array data and Log2 conversion was performed. The GSE49710 cohort contained relatively complete information of age (> 18 months n=193, ≤18 months n=305), sex (boy n=287, girl n=211), COG risk (High risk n=176, Low risk n=322), MYCN status (amplified n=92, non amplified n=401), INSS stage (Stage I n =121, Stage II n=78, Stage III n=63, Stage IV n=183, Stage IVs n=53), survival state (alive n=393, dead n=105) and survival time. E-MTAB-8248 had complete information of age (> 18 months n=119, ≤18 months n=104), MYCN status (amplified n=46, non amplified n=177), INSS stage (Stage I n =29, Stage II n=39, Stage III n=36, Stage IV n=89, Stage IVs n=30), alive state (alive n=181, dead n=42) and survival time. By combining and standardizing transcriptome data in two cohorts, we obtained transcriptome microarray matrix data with 721 samples. The transcriptome data of 39 NB cell lines in the GSE89413 cohort were obtained and Log2 conversion was also performed for further study.

Consensus clustering analysis for Args

We searched the key word “anoikis” in Genecards and Harmonizome databases, and a total of 640 Args were obtained. Univariate cox regression analysis was conducted on the expression level of 640 genes and the survival state of the children in this study. A total of 283 Args closely related to the survival state were obtained according to the standard of HR≠1 and p<0.05. Ten Args with the greatest impact on the prognosis were: PARP1 (HR=8.51), ELAVL1 (HR=6.49), XRCC5 (HR=5.62), CSNK2A1 (HR=5.60), DAP3 (HR=4.96), MYO5A (HR=0.20), CDC42 (HR=0.24), ARHGEF7 (HR=0.25), PDPK1 (HR=0.27), ATF2 (HR=0.27). Ten Args with the most statistical difference were: NTRK1 (p=1.58e-35), ETV4 (p=2.80e-31), PARP1 (p=4.19e-30), CDC42 (p=7.44e-30), SKP2 (p=3.00e-29), CHEK2 (p=6.68e-29), CDKN3 (p=7.78e-29), HK2 (p=8.59e-27), ARHGEF7 (p=1.62e-26), KIF18A (p=3.48e-26). We showed 28 genes with the lowest p-value through the forest map (Figure 2A). According to the expression levels of these 283 genes, 721 samples were clustered into different anoikis-related subtypes by consistent cluster analysis. Through the clustering heatmap and consistent cumulative distribution function (CDF) graph, we found that the CDF slope was the smallest and the clustering effect was the best when the number of clusters k was equal to 2 (Figures 2B, C). Through principal component analysis (PCA), we found that the clustering method could well distinguish the samples into two subtypes (Figure 2D) and KM survival analysis suggested a worse prognosis for ARCluster B (Figure 2E). These results indicated that the two anoikis-related subtypes had significant differences in the prognosis, and aroused strong interests of us in the differential expression of Args and the differential distribution of clinical information between ARCluster B and ARCluster A. We further analyzed the differences of Args between the two subtypes and found 119 differentially expressed Args. The 10 Args with the largest value of |LogFC| were NTRK1 (LogFC=3.66), LMO3 (LogFC=2.47), TWIST1 (LogFC=2.46), BIRC3 (LogFC=2.02), HK2 (LogFC=1.94), CCR7(LogFC=-1.88), ITGA8(LogFC=-1.86), KIF18A(LogFC=1.82), E2F1(LogFC=1.78) and UBE2C(LogFC=1.77). 32 genes with the lowest p-value were shown by boxplot (Figure 2F). It could also be found that the expression levels of Args were different between ARCluster A and ARCluster B according to the heatmap (Figure 2G). The results also showed that ARCluster B with a poor prognosis contained more samples with age > 18 months(p<0.001), INSS stage III-IV(p<0.001), MYCN amplified(p<0.001) and high COG risk(p<0.001) (33, 34). However, ARCluster A with a better prognosis mainly included samples with age ≤18 months(p<0.001), INSS stage I/II/IVs(p<0.001), MYCN non amplified(p<0.001) and low COG risk(p<0.001)(Figure 2G). The above results fully indicated that children with NB could be clustered into ARCluster A with a good prognosis or ARCluster B with a poor prognosis through univariate cox regression analysis and consistent cluster analysis, which was consistent with clinical indicators such as age, COG risk, MYCN status and INSS stage (33, 34). These results made us interested in the study of the important pathways and key genes between the two anoikis-related subtypes.

FIGURE 2
www.frontiersin.org

Figure 2 Consensus clustering analysis for Args (A). Args closely related to the survival state were obtained according to the standard of HR≠1 and P<0.05. 28 genes with the lowest p-value through the forest map. (B, C). Clustering heatmap and consistent cumulative distribution function (CDF) graph of 721 samples (k=2). (D). Distinguish the samples into two subtypes through principal component analysis (PCA). (E) KM survival analysis of two anoikis-related subtypes. (F) Boxplot of 32 Args between the two subtypes with the lowest p-value. (G) Heatmap of clinical information and differentially expressed Args between ARCluster A and ARCluster B. ***, p < 0.001.

GSVA and GSEA analyses and immune cell infiltration in ARClusters

In order to explore whether there were differences in the pathways between ARCluster B and ARCluster A, we performed GSVA and GSEA enrichment analyses for the two subtypes. It was found that there were obvious differences in the pathways between ARCluster B and ARCluster A through GSVA enrichment analysis. According to “c5.go.v7.4.symbols.gmt”, the five pathways with the most obvious differences were as follows: GOCC_METHYLOSOME, GOBP_ACTIN_FILAMENT_ORGANIZATION, GOBP_DNA_STRAND_ELONGATION, GOBP_REGULATION_OF_ACTIN_FILAMENT_ORGANIZATION, GOCC_SPLICEOSOMAL_TRI_SNRNP_COMPLEX (Figure 3A). We found that the expression level of GOBP_REGULATION_OF_CELL_MATRIX_ADHESION pathway decreased significantly in ARCluster B. According to “c2.cp.kegg.v7.4.symbols.gmt”, the most obvious five pathways differentially expressed between two subtypes were: KEGG_BASE_EXCISION_REPAIR, KEGG_SPLICEOSOME, KEGG_NUCLEOTIDE_EXCISION_REPAIR, KEGG_DNA_REPLICATION, KEGG_RNA_DEGRADATION (Figure 3B). We also found that the expression level of “KEGG_CELL_ADHESION_MOLECULES_CAMS” pathway related to cell adhesion was significantly reduced in ARCluster B. GSEA enrichment analysis was further performed. According to “c5.go.v7.4.symbols.gmt”, the five pathways with the most obvious expression differences were GOBP_ADAPTIVE_IMMUNE_RESPONSE, GOBP_IMMUNE_RESPONSE_REGULATING_CELL_SURFACE_RECEPTOR_SIGNALING_PATHWAY, GOBP_POSITIVE_REGULATION_OF_LEUKOCYTE_CELL_CELL_ADHESION, GOBP_LEUKOCYTE_MEDIATED_IMMUNITY, GOBP_LYMPHOCYTE_MEDIATED_IMMUNITY (Figure 3C). According to “c2.cp.kegg.v7.4.symbols.gmt”, the five pathways contained KEGG_HEMATOPOIETIC_CELL_LINEAGE, KEGG_GRAFT_VERSUS_HOST_DISEASE, KEGG_CELL_ADHESION_MOLECULES_CAMS, KEGG_CELL_CYCLE, KEGG_CYTOKINE_CYTOKINE_RECEPTOR_INTERACTION (Figure 3D). Through GSVA and GSEA enrichment analyses, it was found that there were A large number of differentially expressed pathways between ARCluster A and ARCluster B. Moreover, the results showed that the expression levels of the pathways related to cell adhesion were all reduced in ARCluster B, which was consistent with the phenomenon that malignant tumors with AR were more likely to break down the barriers of ECM and form metastasis. We analyzed the differences of the infiltration levels of immune cells between ARCluster A and ARCluster B through ssGSEA. Except for the increased level of Activated. CD4 T cell and Type2 T helper cell, the infiltration level of other immune cells in ARCluster B decreased significantly (Figure 3E). These results fully showed that there were obvious differences in the infiltration of immune cells between different anoikis-related subtypes of NB, and the changes of TME might play an important role in breaking down the barriers of ECM and forming AR.

FIGURE 3
www.frontiersin.org

Figure 3 GSVA and GSEA analyses and immune cell infiltration in different ARClusters. (A, B). Enrichment pathways between ARCluster B and ARCluster A through GSVA enrichment analysis. (C, D). Enrichment pathways between ARCluster B and ARCluster A through GSEA enrichment analysis. (E). Difference of the infiltration levels of immune cells between ARCluster A and ARCluster B through SsGSEA. **, p < 0.01; ***, p < 0.001.

Hub genes in anoikis-related subtypes

We further screened the hub genes that might play important roles in anoikis-related subtypes through WGCNA. The remaining samples were clustered after removing outliers (Figure 4A). When the threshold value was 7, R2 = 0.92, the gene expression was the most consistent with the scale-free network (35) (Figure 4B). Through the construction of scale-free network, we obtained a total of 18 modules closely related to anoikis and further combined the 18 modules into 14 modules according to the standard of Height=0.25 (Figures 4C, D). Among the 14 modules, MElightcyan, MEcyan, MEpink, MEyellow, and MEgrey were positively related to anoikis, while MEturquoise, MEblue, and MEmagenta were negatively related to anoikis. Among all modules, MEturquoise module had the largest number of genes and was the most closely related to anoikis (correlation index=-0.77) (Figure 4E). According to the screening criteria of GS>0,5 and |MM|>0.8, 93 core genes were screened from the MEturquoise module (Figure 4F). We further crossed the prognosis-related Args with the 93 core genes and finally obtained 10 anoikis-related hub genes (Figure 4G): BIRC5 (GS=0.62, MM=-0.81), BUB1 (GS=0.60, MM=-0.80), CDKN3 (GS=0.64, MM=-0.82), CENPF (GS=0.61, MM=-0.82), CHEK2 (GS=0.66, MM=-0.85), E2F1 (GS=0.64, MM=-0.83), MAD2L1 (GS=0.64, MM=-0.81), NTRK1 (GS=-0.67, MM=0.82), PLK1 (GS=0.65, MM=-0.85) and UBE2C (GS=0.60, MM=-0.81).The above results indicated that there were multiple modules affecting anoikis between two anoikis-related subtypes of NB, and the MEturquoise module had the most obvious inhibitory effect on the anoikis of NB. The screening of the hub genes in this module would play an important role in our subsequent research on the mechanism of AR in NB.

FIGURE 4
www.frontiersin.org

Figure 4 Hub genes in anoikis-related subtypes based on WGCNA. (A) Clustering dendrogram of 721 samples. (B) Construction of scale-free network (Soft threshold value=7, R2 = 0.92). (C, D). 18 anoikis-related modules of NB were screened and further combined into 14 modules (Red line in C: MEDissThres = 0.25). (E) Relationships between anoikis and 14 modules. Turquoise module was most closely related to anoikis in NB. (F) Scatterplot analysis of Turquoise module. Key genes were filtered in the upper right region of GS > 0.5 and MM > 0.8. (G) The Venn plot of 10 hub genes in Turquoise module and 283 Args.

Construction of an accurate anoikis-related model

We had found that there were two anoikis-related subtypes with different prognoses in NB through the above cluster analysis. Therefore, we hoped to construct a prognostic model that could classify NB into anoikis-related subtypes with different prognoses through their transcriptome data. Firstly, we randomly divided 721 cases of NB into the training group and the test group with the same sample size, and constructed lasso regression analysis based on the expression levels of 283 Args and the prognoses of these children. We got a total of 26 key genes through further cross-validation (alpha=1). The multivariate cox regression analysis was carried out between the expression of these 26 key genes and the prognoses of these patients. Finally, we obtained 15 key Args for building the model and listed the location of each gene (Figures 5A, B): STK11, CCND1, PIK3CG, DAPK1, FADD, CPT1A, TWIST1, BAG1, RAD9A, MALAT1, NTRK3, MMP3, HTRA2, KIF18A, DAP3. The weight of each gene in this model was stated in Table 1. By calculating the riskscore of each sample, we divided the samples of training group and test group into high-risk group and low-risk group according to the median value of each group. The results suggested that the expression levels of 15 Args participating in the construction of this model were different between high and low-risk groups (Figure 5C), and the riskscores in ARCluster B were significantly higher than those in ARCluster A (Figure 5D). Almost all of the samples in ARCluster B were classified into the high-risk group, and most of the samples in ARCluster A were classified into the low-risk group. Almost all the dead patients came from the high-risk group, while the majority of the surviving children belonged to the low-risk group (Figure 5E). We found that Args differentially expressed in ARCluster A and ARCluster B groups also had consistent differences between high-risk and low-risk groups, and the clinical indicators related to poor prognosis including ≥ 18 months, MYCN amplified, high COG risk, INSS stage III-IV were mainly enriched in the anoikis-related high-risk group (Figure 5F), which was consistent with our previous conclusions. According to the commonly used clinical indicators: Age, MYCN status, COG risk, INSS stage and the risk level calculated by anoikis-related model, we drew the nomogram to comprehensively evaluate the prognoses of children with NB. Different clinical indicators were given different points in the nomogram. The probabilities of the 1-, 3-, and 5-year survival were predicted by the accumulated scores. It was found that the scores of the two opposite states of other indicators were very close except for the risk level. The nomogram showed that the risk level was an independent diagnostic factor and synergistically predicted the survival probability of NB patients. (Figure 5G). KM survival analysis indicated that the prognoses of patients with high riskscores were poor (p<0.001 in three groups) whether in the total sample group, training group or test group (Figure 5H). AUC values at 1-year, 3-years, 5-years in the total sample group (0.947, 0.900, 0.917), training group (0.947, 0.933, 0.958), and test group (0.942, 0.873, 0.888) indicated that the model we built could well predict the prognoses of children with NB (Figure 5I). By drawing the decision curve (Figure 5J) and calibration curve (Figure 5K), we could also find that the anoikis-related model we built could well predict the 1,3,5-year survival rate of children with NB. The cumulative hazard curve showed that the overall survival (OS) risk of children in high-risk group increased gradually with the prolongation of survival time (Figure 5L). The above results suggested that this anoikis-related model could accurately reflect the prognoses of children with NB, which laid the foundation for the popularization of this model and the researches on the potential mechanisms and immune microenvironment of NB in high-risk group.

FIGURE 5
www.frontiersin.org

Figure 5 Construction of an accurate anoikis-related model. (A) 26 key genes screened through further cross-validation (alpha=1). (B) The locations of 15 genes in chromosomes obtained by multivariate cox regression analysis. (C) Heatmap of the distribution of the 15 genes involved in model construction between the high riskscore and low riskscore groups. (D) Boxplot of riskscore distribution of samples in different ARClusters. (E) Alluvial diagram of subtypes and living status. (F) Heatmap of clinical information and differentially expressed Args between two risk subgroups. (G) Nomogram of the commonly used clinical indicators and the risk level calculated by anoikis-related model. (H) OS analysis of two risk subgroups in the total sample group, training group and test group. (I) The time-dependent ROC curves for OS at 1-, 3-, and 5-years in the total sample group, training group and test group. (J) Net benefit (y-axis) as calculated are plotted against the threshold probabilities of patients having 5-year survival on the x-axis. (K) Calibration curve of OS at 1-, 3-, and 5-years for the anoikis-related model in the total sample group. (L). Cumulative hazard curve represented the probability of survival risk with the prolongation of survival time. ***, p < 0.001.

TABLE 1
www.frontiersin.org

Table 1 15 key genes involved in anoikis-related prognostic model construction.

Validation of the anoikis-related model

We found that the model constructed in this study could effectively divide the samples of GSE49710 and E-MTAB-8248 cohorts into high-risk group and low-risk group (Figures 6A, B), and the prognoses of children in high-risk group were worse than that in low-risk group (Figures 6C, D). Args differentially expressed between ARCluster A and ARCluster B were also differentially expressed between high and low-risk groups of the GSE49710 and E-MTAB-8248 cohorts according to the heatmap (Figures 6A, B). Clinical indicators closely related to poor prognoses of children: ≥18 months, MYCN amplified, high COG risk, INSS stage III-IV were also enriched in high-risk group of both cohorts (Figures 6A, B). The results of independent prognostic analysis suggested that sex, MYCN status, COG risk, INSS stage and anoikis-related riskscore could be used as independent prognostic indicators in the GSE49710 cohort (Figure 6E). Age, MYCN status and anoikis-related riskscore were independent prognostic indicators in the E-MTAB-8248 cohort (Figure 6F). 39 NB cell lines in GSE89413 cohort were also scored according to the constructed model. RD, NB-SD, COG-N-519, CHP-212, COG-N-561, IMR-05, SK-N-FI, COG-N-440, LA-N-6, RD, NB-SD, COG-N-519, CHP-212, COG-N-561, IMR-05, SK-N-Fi, COG-N-440, LA-N-6 and Nb-1 were the top ten NB cell lines with the highest riskscore, while COG-N-557nb, SK-N-SH, NB-1643, IMR-32, COG-N-471nb, NB-69, NMB, SMS-KAN, SK-N-BE (2) and SH-SY5Y were the top ten cell lines with the lowest riskscore (Figure 6G). The above results fully indicated that the anoikis-related model constructed by us could generally predict the prognoses of children with NB, and the scores obtained by this model could be used as an independent indicator to better predict the prognosis. The above results also suggested that there were generally two anoikis-related subtypes in the population and cell lines of NB, among which the subtype with high riskscore and poor prognosis was likely to be closely related to the AR. The classification of NB cell lines would help us study the mechanism of AR in NB.

FIGURE 6
www.frontiersin.org

Figure 6 Validation of the anoikis-related model. (A, B). Heatmap of clinical information and differentially expressed Args between two risk subgroups in GSE49710 and E-MTAB-8248. (C, D): KM survival analysis of two risk subgroups in GSE49710 and E-MTAB-8248. (E, F): Forest plots of multivariable cox regression analyses of the clinical features and riskscores calculated by the anoikis-related model in GSE49710 and E-MTAB-8248. (G): Riskscores of 39 NB cell lines in GSE89413 cohort according to the constructed model. *, p < 0.05; **, p < 0.01; ***, p < 0.001.

Correlation between immune cell status and different riskscores

Previous studies had shown that TME was closely related to the metastasis of tumor cells, and AR was an important mechanism for metastasis. Therefore, understanding the TME characteristics of anoikis-related subtypes of NB was very important to immunotherapy for anoikis-related subtype of NB with a poor prognosis. CIBERSORT algorithm was used to analyze the immune cell infiltration level of NB samples from different anoikis-related risk groups. The results showed that B cells memory, Plasma cells, T cells follicular helper, NK cells resting, and Neutrophils were highly expressed, while T cells CD4 memory resting, T cells gamma delta, Macrophages M2 and Mast cells activated were weekly expressed in anoikis-related subtype with high riskscore (Figure 7A). We also found that there were correlations between immune cells, especially positive relationship between B cells naive and T cells follicular helper (correlation index=0.37). But the negative relationship between B cells naive and Macrophages M0 was more obvious (correlation index=-0.51) (Figure 7B). The estimate algorithm was used to further analyze the differences of TME between different risk groups. The results showed that the estimate score of high-risk group was lower, and the content of stromal cells and immune cells in tumor tissue was lower than that of low-risk group (Figure 7C). By analyzing the correlations between 15 Args participating in the model construction and immune cells, we found that except STK11 and CCND1, the other 13 Args were all related to the infiltration levels of immune cells. In particular, the expression level of PIK3CG was positively correlated with the infiltration levels of T cells regulatory (Tregs), T cells follicular helper, T cells CD8, T cells CD4 naive, T cells CD4 memory activated, Plasma cells, Monocells, Macrophages M1, B cells naive, and negatively correlated with the expression levels of NK cells activated, Mast cells activated, Macrophages M2, Macrophages M0 (Figure 7D). By further analyzing the correlation between the infiltration levels of 23 kinds of immune cells and the riskscores of patients, we found that a total of 6 kinds of immune cells were related to the riskscore. B cells memory (R=0.13, p=0.0079), Plasma cells (R=0.23, p=1.8E-6) and T cells follicular helper (R=0.12, p=0.0013) were positively correlated with the riskscore, while T cells CD4 memory resting (R=-0.27, p=1.2E-8), Macrophages M2 (R=-0.12, p=0.0012) and Mast cells activated (R=-0.12, p=0.0016) were negatively correlated with the riskscore (Figure 7E). The above results fully demonstrated that TME and immune cell infiltration level were different between anoikis-related subtypes. The key genes involved in the construction of anoikis-related model were also related to the infiltration levels of immune cells. The differences of the TME and immune cell infiltration level might be important factors leading to the AR of NB, and the key genes for constructing the model might play important roles in the above mechanisms.

FIGURE 7
www.frontiersin.org

Figure 7 Correlation between immune cell status and different riskscores. (A) Analyzing the immune cell infiltration levels of NB samples from different anoikis-related risk groups by CIBERSORT algorithm. (B) Correlations between immune cells. (C) Differences of TME between different risk groups. (D) Analyzing the correlations between 15 Args participating in the model construction and immune cells. (E) Relationships between 6 kinds of immune cells and the riskscores of patients. *, p < 0.05; **, p < 0.01; ***, p < 0.001.

Chemotherapeutic response analysis

According to the results calculated based on the GDSC database, we predicted that the anoikis-related high-risk group and low-risk group had different drug sensitivities to 158 kinds of chemotherapeutic drugs. The low-risk group had higher sensitivity to most of the predicted drugs, while the high-risk group had higher sensitivity to only17 kinds of drugs. The 8 kinds of chemotherapeutic drugs with the highest drug sensitivity in high-risk group were: JQ1, GSK269962A, Doramapimod, BMS-754807, AZD8055, KU-55933, AZD6482 and Axitinib (Figure 8). The above results indicated that these selected chemotherapeutic drugs might have better therapeutic effects on NB in high-risk group, and might have potential targeted therapeutic effects on the AR of NB.

FIGURE 8
www.frontiersin.org

Figure 8 Chemotherapeutic Response Analysis. (A-H). The eight predicted chemotherapeutic drugs with the highest drug sensitivity in the high-risk group.

Discussion

NB is the most common solid extra-cranial malignancy in infants (36). At present, there are many treatments for NB, including surgery (4, 5), chemotherapy (7, 8), bone marrow transplantation (10), immunotherapy (1116) and so on. However, the five-year survival rate is still low even with multiple therapies for high-risk NB, especially those with metastasis. Therefore, exploring the mechanisms of metastasis has a high clinical significance for the treatment of NB. In recent years, as an important mechanism of metastasis, AR has become a hot topic of scientific research and has been paid more and more attentions by scholars. Anoikis can prevent normal cells and tumor cells from leaving the ECM, suppress their distant colonization and growth, while malignant tumor cells colonizing in organs other than the primary site can get rid of the regulation of anoikis. In recent years, the AR has been found in many solid tumors, including breast cancer (37, 38), prostate cancer (39, 40), gastric cancer (41, 42), liver cancer (43, 44), renal cancer (45, 46) and so on. Some scholars have also found this phenomenon in NB (31, 32), but there is a lack of research on the relevant mechanism. In this study, transcriptome microarray data of children with NB were used to construct an anoikis-related prognostic model for the first time. Studies had shown that the model constructed by us could well predict the prognoses of children and provide an auxiliary role for the study of mechanisms related to the AR.

This study followed the principle of large sample size and complete clinical information when we selected the cohorts, and finally included 721 cases for consistency cluster analysis. The Args obtained were also strictly screened, and the genes less related to anoikis were removed, which ensured that the constructed model had high accuracy and high correlation with anoikis. The results of consistent cluster analysis showed that there were two anoikis-related subtypes with completely opposite prognosis in the population of NB. Through the analysis of the distribution of clinical characteristics of different subtypes, we found that the clinical information indicating a poor prognosis was enriched in the poor prognostic subtype, and the distribution of Args in the two subtypes was also significantly different. By dividing samples into these two subtypes, we could more effectively and accurately analyze the changes of core genes, key pathways and immune microenvironment related to the anoikis of NB. But a small number of children with poor prognostic indicators were enrolled in the good prognostic subtype. We thought that the expressions of Args could distinguish between more and less aggressive MYCN or stage III-IV tumors, or maybe a certain subgroup of aggressive tumors could be predicted based on Args.

Through the enrichment analysis of GSVA and GSEA for the two subtypes, we found that there were significant differences in the expression levels of pathways between the two subtypes, and the expression level was significantly reduced especially in the pathway of cell adhesion, which was consistent with the phenomenon of AR. The analysis of these differentially expressed pathways will help us better understand the important roles of AR in the metastasis of NB. By analyzing the immune cell infiltration of the two subtypes, we speculated that the occurrence of anoikis might be closely related to the TME of NB, especially the changes of immune cell infiltration. We also found that there were multiple modules that played important roles in anoikis in the subtype with a poor prognosis through WGCNA analysis. Ten hub Args that played important roles in the mechanism of anoikis in NB were obtained. Anoikis is actually a kind of responses of cells to their loss of contact by different signal pathways, and cell cycle regulation is also an important way for cells to generate anoikis or AR (47, 48). Previous studies have mentioned that BIRC5, BUB1, CDKN3, CENPF, PLK1, CHEK2, Ntrk1, E2F1 and UBE2C are closely related to cell cycle. However, the functions of these genes are not single and these genes may play other functions such as promting anoikis or AR in additon to regulating cell cycle. For example, Survivin (BIRC5) is required for enhancing AR in ovarian cancer cells and hepatocellular carcinoma cells (24, 49). CHEK2 is a mediator of anoikis of intestinal epithelial cells and is associated with the progression of papillary thyroid cancer by promoting anoikis (50, 51). The protein encoded by E2F1 which plays a crucial role in the control of cell cycle is a member of the E2F family of transcription factors. However, current studies have also found that E2F1 is closely related to anoikis (52), and CDKN3, as the upstream of E2F1, regulates the expression of E2F1 (53). The PI3K/Akt signal pathway is one of the most important downstream targets of TrkA translated by Ntrk1 (54, 55), and mTOR which is thought to be directly involved in anoikis is the downstream target of PI3K/Akt pathway (56, 57). Studies have shown that PLK1 is closely related to cell mitosis, but recent studies have also found that PLK1 protects esophageal cancer cells from anoikis through regulating β-Catenin protein levels by inhibiting their degradation (58). BUB1 can form polymer with PLK1 and promote the function of PLK1 (59). Recent studies have also found that UBE2C/ZEB1/2 signal axis plays an important role in the metastasis and AR in cervical cancer (60). Silencing of CENPF could simultaneously improve sensitivity to anoikis-induced apoptosis in human PC3 cells (61). All these results indicated that the hub genes involved in this study were indeed closely related to anoikis.The study of these genes will help us better understand the role of hub genes in the mechanisms of anoikis in NB.

A large number of previous studies have shown that TME not only plays an important role in the metastasis of NB (62, 63), but also is closely related to the drug resistance (6466), especially immunotherapy drugs (65, 66). Therefore, the study of TME in the anoikis-related subtypes will help us better understand the mechanism of AR in NB and guide the screening of targeted therapeutic drugs. Through the analysis of TME, we found that the infiltration levels of some immune cells in the two groups were significantly different. The infiltration levels of some immune cells were also correlated with the riskscore of each sample, and 13 of the 15 genes involved in the model construction were correlated with the infiltration levels of some immune cells. These findings with statistical differences will provide guidance for us to further analyze the correlation between TME and AR of NB.

At present, we still face many challenges for the clinical treatments of NB, especially for the refractory and relapsed high-risk NB (67, 68), which is also the focus of scientific research. In this study, we found that there was a subtype of NB with a poor prognosis, which was closely related to anoikis. If there are effective chemotherapeutic drugs or targeted therapeutic drugs, the prognosis will be greatly improved for these children. Through drug sensitivity prediction, we found that most of the drugs with statistical significance were sensitive to the subtype of NB with good prognosis, while only a small number of drugs were sensitive to the poor prognosis subtype of NB. At present, there is no reseach on directly applying any of the eight drugs obtained to the study of anoikis. However, the targets of some drugs may be related to anoikis. For example, the PI3K/Akt/mTOR signaling pathway promotes the development of AR (57). Current studies have found that AZD6482 is a selective inhibitor of PI3K (69, 70). Recent studies have shown that AZD8055, as a potent dual mTORC1-mTORC2 inhibitor (71), can inhibit NB proliferation in vitro and in vivo by targeting mTOR (72). Overexpression of mTOR in non-transformed wt fibroblasts promoted AR, and mTOR knockdown in MEFs deficient in all retinoblastoma-family members restored anoikis (73). GSK269962A is a selective inhibitor of ROCK (74), and multiple studies have shown that ROCK is closely related to AR (75, 76). Doramapimod is a selective inhibitor of P58 MAPK (77) and ERK/MAPK signaling was found to operate downstream of ErbB2 and EGFR to protect cells from anoikis (78). Some studies have shown that IGF-1R is associated with AR (79, 80), and BMS-754807 is an important selective inhibitor of IGF-1R (81). Recent study has also found that foretinib, as an orally available multikinase inhibitor of c-Met and VEGFR-2, blocks proliferation, induces anoikis and impairs ovarian cancer metastasis (82). Axitinib, as a selective inhibitor of VEGFR, may inhibit AR (83). In the later stage, further experiments in vivo and in vitro on these drugs may help us better understand the relevant mechanisms of AR in NB and guide the selection of chemotherapeutic drugs and molecular targeted drugs for the anoikis-related NB subtype with poor prognosis.

Conclusion

In conclusion, this study clustered NB patients into two anoikis-related subtypes with different prognoses based on consistency clustering analysis, and analyzed the differences in clinical information, hub genes, signal pathways and immune microenvironment between the two subtypes. Subsequently, we constructed an anoikis-related model that could evaluate the prognoses of children through 15 Args, further verified the accuracy, popularization and application of this model, analyzed the differences of TME between two risk groups, and screened 17 drugs with high drug sensitivity in high-risk group. The model constructed in this study has a good guiding significance for the risk staging and prognosis evaluation of NB, a high auxiliary role in the study of the mechanism related to the AR, and a good reference significance for the selection of chemotherapeutic and targeted therapeutic drugs for NB.

Data availability statement

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

Author contributions

LH, BJ and YF: conceptualization, methodology, writing-reviewing and editing. JC, MS, CC and MK: investigation, data curation, writing-original draft preparation. JC, CC, BQ, JS, XM and JZ: visualization, validation, supervision, software. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the General Project of Nanjing Medical University (NMUB20210060), National Natural Science Foundation of China (81903383), Natural Science Foundation of Jiangsu Province (BK20211009), Scientific Research Projects of Jiangsu Health Commission (ZDB2020018), China Postdoctoral Science Foundation funded project (2021M701764), Special Fund for Health Science and Technology Development in Nanjing (JQX19008), Nanjing Medical Science and Technology Development Project (YKK21149), Young Talent Support Project of Children’s Hospital of Nanjing Medical University (TJGC2020016, TJGC2020007, TJGC2020014).

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1135617/full#supplementary-material

References

1. Maris JM, Hogarty MD, Bagatell R, Cohn SL. Neuroblastoma. Lancet (2007) 369(9579):2106–20. doi: 10.1016/s0140-6736(07)60983-0

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Strother DR, London WB, Schmidt ML, Brodeur GM, Shimada H, Thorner P, et al. Outcome after surgery alone or with restricted use of chemotherapy for patients with low-risk neuroblastoma: Results of children's oncology group study P9641. J Clin Oncol (2012) 30(15):1842–8. doi: 10.1200/jco.2011.37.9990

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Twist CJ, Schmidt ML, Naranjo A, London WB, Tenney SC, Marachelian A, et al. Maintaining outstanding outcomes using response- and biology-based therapy for intermediate-risk neuroblastoma: A report from the children's oncology group study ANBL0531. J Clin Oncol (2019) 37(34):3243–55. doi: 10.1200/jco.19.00919

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Vollmer K, Gfroerer S, Theilen TM, Bochennek K, Klingebiel T, Rolle U, et al. Radical surgery improves survival in patients with stage 4 neuroblastoma. World J Surg (2018) 42(6):1877–84. doi: 10.1007/s00268-017-4340-9

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Holmes K, Pötschger U, Pearson ADJ, Sarnacki S, Cecchetto G, Gomez-Chacon J, et al. Influence of surgical excision on the survival of patients with stage 4 high-risk neuroblastoma: A report from the HR-NBL1/SIOPEN study. J Clin Oncol (2020) 38(25):2902–15. doi: 10.1200/jco.19.03117

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Matthay KK. Chemotherapy for neuroblastoma: does it hit the target? Lancet Oncol (2008) 9(3):195–6. doi: 10.1016/s1470-2045(08)70046-9

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Furman WL, McCarville B, Shulkin BL, Davidoff A, Krasin M, Hsu CW, et al. Improved outcome in children with newly diagnosed high-risk neuroblastoma treated with chemoimmunotherapy: Updated results of a phase II study using hu14.18K322A. J Clin Oncol (2022) 40(4):335–44. doi: 10.1200/jco.21.01375

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Garaventa A, Poetschger U, Valteau-Couanet D, Luksch R, Castel V, Elliott M, et al. Randomized trial of two induction therapy regimens for high-risk neuroblastoma: HR-NBL1.5 international society of pediatric oncology European neuroblastoma group study. J Clin Oncol (2021) 39(23):2552–63. doi: 10.1200/jco.20.03144

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Liu KX, Naranjo A, Zhang FF, DuBois SG, Braunstein SE, Voss SD, et al. Prospective evaluation of radiation dose escalation in patients with high-risk neuroblastoma and gross residual disease after surgery: A report from the children's oncology group ANBL0532 study. J Clin Oncol (2020) 38(24):2741–52. doi: 10.1200/jco.19.03316

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Park JR, Kreissman SG, London WB, Naranjo A, Cohn SL, Hogarty MD, et al. Effect of tandem autologous stem cell transplant vs single transplant on event-free survival in patients with high-risk neuroblastoma: A randomized clinical trial. Jama (2019) 322(8):746–55. doi: 10.1001/jama.2019.11642

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Ladenstein R, Pötschger U, Valteau-Couanet D, Luksch R, Castel V, Yaniv I, et al. Interleukin 2 with anti-GD2 antibody ch14.18/CHO (dinutuximab beta) in patients with high-risk neuroblastoma (HR-NBL1/SIOPEN): a multicentre, randomised, phase 3 trial. Lancet Oncol (2018) 19(12):1617–29. doi: 10.1016/s1470-2045(18)30578-3

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Mabe NW, Huang M, Dalton GN, Alexe G, Schaefer DA, Geraghty AC, et al. Transition to a mesenchymal state in neuroblastoma confers resistance to anti-GD2 antibody via reduced expression of ST8SIA1. Nat Cancer (2022) 3(8):976–93. doi: 10.1038/s43018-022-00405-x

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Theruvath J, Menard M, Smith BAH, Linde MH, Coles GL, Dalton GN, et al. Anti-GD2 synergizes with CD47 blockade to mediate tumor eradication. Nat Med (2022) 28(2):333–44. doi: 10.1038/s41591-021-01625-x

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Heitzeneder S, Bosse KR, Zhu Z, Zhelev D, Majzner RG, Radosevich MT, et al. GPC2-CAR T cells tuned for low antigen density mediate potent activity against neuroblastoma without toxicity. Cancer Cell (2022) 40(1):53–69.e9. doi: 10.1016/j.ccell.2021.12.005

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Tumino N, Weber G, Besi F, Del Bufalo F, Bertaina V, Paci P, et al. Polymorphonuclear myeloid-derived suppressor cells impair the anti-tumor efficacy of GD2.CAR T-cells in patients with neuroblastoma. J Hematol Oncol (2021) 14(1):191. doi: 10.1186/s13045-021-01193-0

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Heczey A, Courtney AN, Montalbano A, Robinson S, Liu K, Li M, et al. Anti-GD2 CAR-NKT cells in patients with relapsed or refractory neuroblastoma: an interim analysis. Nat Med (2020) 26(11):1686–90. doi: 10.1038/s41591-020-1074-2

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Irwin MS, Naranjo A, Zhang FF, Cohn SL, London WB, Gastier-Foster JM, et al. Revised neuroblastoma risk classification system: A report from the children's oncology group. J Clin Oncol (2021) 39(29):3229–41. doi: 10.1200/jco.21.00278

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Gerarduzzi C, Hartmann U, Leask A, Drobetsky E. The matrix revolution: Matricellular proteins and restructuring of the cancer microenvironment. Cancer Res (2020) 80(13):2705–17. doi: 10.1158/0008-5472.Can-18-2098

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Janiszewska M, Primi MC, Izard T. Cell adhesion in cancer: Beyond the migration of single cells. J Biol Chem (2020) 295(8):2495–505. doi: 10.1074/jbc.REV119.007759

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Han HJ, Sung JY, Kim SH, Yun UJ, Kim H, Jang EJ, et al. Fibronectin regulates anoikis resistance via cell aggregate formation. Cancer Lett (2021) 508:59–72. doi: 10.1016/j.canlet.2021.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Adeshakin FO, Adeshakin AO, Afolabi LO, Yan D, Zhang G, Wan X. Mechanisms for modulating anoikis resistance in cancer and the relevance of metabolic reprogramming. Front Oncol (2021) 11:626577. doi: 10.3389/fonc.2021.626577

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Jin L, Chun J, Pan C, Kumar A, Zhang G, Ha Y, et al. The PLAG1-GDH1 axis promotes anoikis resistance and tumor metastasis through CamKK2-AMPK signaling in LKB1-deficient lung cancer. Mol Cell (2018) 69(1):87–99.e7. doi: 10.1016/j.molcel.2017.11.025

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Sharma R, Gogoi G, Saikia S, Sharma A, Kalita DJ, Sarma A, et al. BMP4 enhances anoikis resistance and chemoresistance of breast cancer cells through canonical BMP signaling. J Cell Commun Signal (2022) 16(2):191–205. doi: 10.1007/s12079-021-00649-9

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Mak CS, Yung MM, Hui LM, Leung LL, Liang R, Chen K, et al. MicroRNA-141 enhances anoikis resistance in metastatic progression of ovarian cancer through targeting KLF12/Sp1/survivin axis. Mol Cancer (2017) 16(1):11. doi: 10.1186/s12943-017-0582-2

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Quail DF, Joyce JA. Microenvironmental regulation of tumor progression and metastasis. Nat Med (2013) 19(11):1423–37. doi: 10.1038/nm.3394

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Deepak KGK, Vempati R, Nagaraju GP, Dasari VR, Nagini S, Rao DN, et al. Tumor microenvironment: Challenges and opportunities in targeting metastasis of triple negative breast cancer. Pharmacol Res (2020) 153:104683. doi: 10.1016/j.phrs.2020.104683

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Mao X, Xu J, Wang W, Liang C, Hua J, Liu J, et al. Crosstalk between cancer-associated fibroblasts and immune cells in the tumor microenvironment: New findings and future perspectives. Mol Cancer (2021) 20(1):131. doi: 10.1186/s12943-021-01428-1

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Chen D, Zhang X, Li Z, Zhu B. Metabolic regulatory crosstalk between tumor microenvironment and tumor-associated macrophages. Theranostics (2021) 11(3):1016–30. doi: 10.7150/thno.51777

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Terrén I, Orrantia A, Vitallé J, Zenarruzabeitia O, Borrego F. NK cell metabolism and tumor microenvironment. Front Immunol (2019) 10:2278. doi: 10.3389/fimmu.2019.02278

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Downs-Canner SM, Meier J, Vincent BG, Serody JS. B cell function in the tumor microenvironment. Annu Rev Immunol (2022) 40:169–93. doi: 10.1146/annurev-immunol-101220-015603

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Leblond P, Dewitte A, Le Tinier F, Bal-Mahieu C, Baroncini M, Sarrazin T, et al. Cilengitide targets pediatric glioma and neuroblastoma cells through cell detachment and anoikis induction. Anticancer Drugs (2013) 24(8):818–25. doi: 10.1097/CAD.0b013e328362edc5

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Bozzo C, Sabbatini M, Tiberio R, Piffanelli V, Santoro C, Cannas M. Activation of caspase-8 triggers anoikis in human neuroblastoma cells. Neurosci Res (2006) 56(2):145–53. doi: 10.1016/j.neures.2006.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Shohet J, Foster J. Neuroblastoma. Bmj (2017) 357:j1863. doi: 10.1136/bmj.j1863

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Pinto NR, Applebaum MA, Volchenboum SL, Matthay KK, London WB, Ambros PF, et al. Advances in risk classification and treatment strategies for neuroblastoma. J Clin Oncol (2015) 33(27):3008–17. doi: 10.1200/jco.2014.59.4648

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Langfelder P, Horvath S. WGCNA: an r package for weighted correlation network analysis. BMC Bioinf (2008) 9:559. doi: 10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

36. Ward E, DeSantis C, Robbins A, Kohler B, Jemal A. Childhood and adolescent cancer statistics, 2014. CA Cancer J Clin (2014) 64(2):83–103. doi: 10.3322/caac.21219

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Luey BC, May FE. Insulin-like growth factors are essential to prevent anoikis in oestrogen-responsive breast cancer cells: Importance of the type I IGF receptor and PI3-kinase/Akt pathway. Mol Cancer (2016) 15:8. doi: 10.1186/s12943-015-0482-2

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Yu SJ, Hu JY, Kuang XY, Luo JM, Hou YF, Di GH, et al. MicroRNA-200a promotes anoikis resistance and metastasis by targeting YAP1 in human breast cancer. Clin Cancer Res (2013) 19(6):1389–99. doi: 10.1158/1078-0432.Ccr-12-1959

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Hasanuzzaman M, Kutner R, Agha-Mohammadi S, Reiser J, Sehgal I. A doxycycline-inducible urokinase receptor (uPAR) upregulates uPAR activities including resistance to anoikis in human prostate cancer cell lines. Mol Cancer (2007) 6:34. doi: 10.1186/1476-4598-6-34

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Sakamoto S, Kyprianou N. Targeting anoikis resistance in prostate cancer metastasis. Mol Aspects Med (2010) 31(2):205–14. doi: 10.1016/j.mam.2010.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Zhang T, Wang B, Su F, Gu B, Xiang L, Gao L, et al. TCF7L2 promotes anoikis resistance and metastasis of gastric cancer by transcriptionally activating PLAUR. Int J Biol Sci (2022) 18(11):4560–77. doi: 10.7150/ijbs.69933

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Du S, Yang Z, Lu X, Yousuf S, Zhao M, Li W, et al. Anoikis resistant gastric cancer cells promote angiogenesis and peritoneal metastasis through C/EBPβ-mediated PDGFB autocrine and paracrine signaling. Oncogene (2021) 40(38):5764–79. doi: 10.1038/s41388-021-01988-y

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Song J, Liu Y, Liu F, Zhang L, Li G, Yuan C, et al. The 14-3-3σ protein promotes HCC anoikis resistance by inhibiting EGFR degradation and thereby activating the EGFR-dependent ERK1/2 signaling pathway. Theranostics (2021) 11(3):996–1015. doi: 10.7150/thno.51646

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Mo CF, Li J, Yang SX, Guo HJ, Liu Y, Luo XY, et al. IQGAP1 promotes anoikis resistance and metastasis through Rac1-dependent ROS accumulation and activation of Src/FAK signalling in hepatocellular carcinoma. Br J Cancer (2020) 123(7):1154–63. doi: 10.1038/s41416-020-0970-z

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Wang L, Li C, Wang J, Yang G, Lv Y, Fu B, et al. Transformable ECM deprivation system effectively suppresses renal cell carcinoma by reversing anoikis resistance and increasing chemotherapy sensitivity. Adv Mater (2022) 34(43):e2203518. doi: 10.1002/adma.202203518

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Chen Z, Liu X, Zhu Z, Chen J, Wang C, Chen X, et al. A novel anoikis-related prognostic signature associated with prognosis and immune infiltration landscape in clear cell renal cell carcinoma. Front Genet (2022) 13:1039465. doi: 10.3389/fgene.2022.1039465

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Kim HR, Lin HM, Biliran H, Raz A. Cell cycle arrest and inhibition of anoikis by galectin-3 in human breast epithelial cells. Cancer Res (1999) 59(16):4148–54.

PubMed Abstract | Google Scholar

48. Macabenta F, Sun HT, Stathopoulos A. BMP-gated cell-cycle progression drives anoikis during mesenchymal collective migration. Dev Cell (2022) 57(14):1683–93.e3. doi: 10.1016/j.devcel.2022.05.017

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Hu H, Li Z, Chen J, Wang D, Ma J, Wang W, et al. P16 reactivation induces anoikis and exhibits antitumour potency by downregulating akt/survivin signalling in hepatocellular carcinoma cells. Gut (2011) 60(5):710–21. doi: 10.1136/gut.2010.220020

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Zhao W, Chen S, Hou X, Chen G, Zhao Y. CHK2 promotes anoikis and is associated with the progression of papillary thyroid cancer. Cell Physiol Biochem (2018) 45(4):1590–602. doi: 10.1159/000487724

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Yoo BH, Berezkin A, Wang Y, Zagryazhskaya A, Rosen KV. Tumor suppressor protein kinase Chk2 is a mediator of anoikis of intestinal epithelial cells. Int J Cancer (2012) 131(2):357–66. doi: 10.1002/ijc.26368

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Zhang YF, Zhang AR, Zhang BC, Rao ZG, Gao JF, Lv MH, et al. MiR-26a regulates cell cycle and anoikis of human esophageal adenocarcinoma cells through Rb1-E2F1 signaling pathway. Mol Biol Rep (2013) 40(2):1711–20. doi: 10.1007/s11033-012-2222-7

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Cen J, Liang Y, Huang Y, Pan Y, Shu G, Zheng Z, et al. Circular RNA circSDHC serves as a sponge for miR-127-3p to promote the proliferation and metastasis of renal cell carcinoma via the CDKN3/E2F1 axis. Mol Cancer (2021) 20(1):19. doi: 10.1186/s12943-021-01314-w

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Yu X, Qi Y, Zhao T, Fang J, Liu X, Xu T, et al. NGF increases FGF2 expression and promotes endothelial cell migration and tube formation through PI3K/Akt and ERK/MAPK pathways in human chondrocytes. Osteoarthritis Cartilage (2019) 27(3):526–34. doi: 10.1016/j.joca.2018.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Song QQ, Lin LP, Chen YL, Qian JC, Wei K, Su JW, et al. Characterization of LTr1 derived from cruciferous vegetables as a novel anti-glioma agent via inhibiting TrkA/PI3K/AKT pathway. Acta Pharmacol Sin (2022). doi: 10.1038/s41401-022-01033-y

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Ng TL, Leprivier G, Robertson MD, Chow C, Martin MJ, Laderoute KR, et al. The AMPK stress response pathway mediates anoikis resistance through inhibition of mTOR and suppression of protein synthesis. Cell Death Differ (2012) 19(3):501–10. doi: 10.1038/cdd.2011.119

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Ma B, Sen T, Asnaghi L, Valapala M, Yang F, Hose S, et al. βA3/A1-crystallin controls anoikis-mediated cell death in astrocytes by modulating PI3K/AKT/mTOR and ERK survival pathways through the PKD/Bit1-signaling axis. Cell Death Dis (2011) 2(10):e217. doi: 10.1038/cddis.2011.100

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Lin DC, Zhang Y, Pan QJ, Yang H, Shi ZZ, Xie ZH, et al. PLK1 is transcriptionally activated by NF-κB during cell detachment and enhances anoikis resistance through inhibiting β-catenin degradation in esophageal squamous cell carcinoma. Clin Cancer Res (2011) 17(13):4285–95. doi: 10.1158/1078-0432.Ccr-10-3236

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Jia L, Li B, Yu H. The Bub1-Plk1 kinase complex promotes spindle checkpoint signalling through Cdc20 phosphorylation. Nat Commun (2016) 7:10818. doi: 10.1038/ncomms10818

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Chen M, Liu LX. MiR-525-5p repressed metastasis and anoikis resistance in cervical cancer via blocking UBE2C/ZEB1/2 signal axis. Dig Dis Sci (2020) 65(8):2442–51. doi: 10.1007/s10620-019-05916-9

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Shahid M, Lee MY, Piplani H, Andres AM, Zhou B, Yeon A, et al. (CENPF), a microtubule binding protein, modulates cancer metabolism by regulating pyruvate kinase M2 phosphorylation signaling. Cell Cycle (2018) 17(24):2802–18. doi: 10.1080/15384101.2018.1557496

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Quinn CH, Beierle AM, Beierle EA. Artificial tumor microenvironments in neuroblastoma. Cancers (Basel) (2021) 13(7). doi: 10.3390/cancers13071629

CrossRef Full Text | Google Scholar

63. Horwacik I. The extracellular matrix and neuroblastoma cell communication-a complex interplay and its therapeutic implications. Cells (2022) 11(19). doi: 10.3390/cells11193172

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Liu H, Zhao H, Sun Y. Tumor microenvironment and cellular senescence: Understanding therapeutic resistance and harnessing strategies. Semin Cancer Biol (2022) 86(Pt 3):769–81. doi: 10.1016/j.semcancer.2021.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Dai E, Zhu Z, Wahed S, Qu Z, Storkus WJ, Guo ZS. Epigenetic modulation of antitumor immunity for improved cancer immunotherapy. Mol Cancer (2021) 20(1):171. doi: 10.1186/s12943-021-01464-x

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Wu HW, Sheard MA, Malvar J, Fernandez GE, DeClerck YA, Blavier L, et al. Anti-CD105 antibody eliminates tumor microenvironment cells and enhances anti-GD2 antibody immunotherapy of neuroblastoma with activated natural killer cells. Clin Cancer Res (2019) 25(15):4761–74. doi: 10.1158/1078-0432.Ccr-18-3358

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Mody R, Naranjo A, Van Ryn C, Yu AL, London WB, Shulkin BL, et al. Irinotecan-temozolomide with temsirolimus or dinutuximab in children with refractory or relapsed neuroblastoma (COG ANBL1221): An open-label, randomised, phase 2 trial. Lancet Oncol (2017) 18(7):946–57. doi: 10.1016/s1470-2045(17)30355-8

PubMed Abstract | CrossRef Full Text | Google Scholar

68. DuBois SG, Granger MM, Groshen S, Tsao-Wei D, Ji L, Shamirian A, et al. Randomized phase II trial of MIBG versus MIBG, vincristine, and irinotecan versus MIBG and vorinostat for patients with relapsed or refractory neuroblastoma: A report from NANT consortium. J Clin Oncol (2021) 39(31):3506–14. doi: 10.1200/jco.21.00703

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Tian J, Zhu Y, Rao M, Cai Y, Lu Z, Zou D, et al. N(6)-methyladenosine mRNA methylation of PIK3CB regulates AKT signalling to promote PTEN-deficient pancreatic cancer progression. Gut (2020) 69(12):2180–92. doi: 10.1136/gutjnl-2019-320179

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Zhou H, Meng M, Wang Z, Zhang H, Yang L, Li C, et al. The role of m5C-related lncRNAs in predicting overall prognosis and regulating the lower grade glioma microenvironment. Front Oncol (2022) 12:814742. doi: 10.3389/fonc.2022.814742

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Jordan NJ, Dutkowski CM, Barrow D, Mottram HJ, Hutcheson IR, Nicholson RI, et al. Impact of dual mTORC1/2 mTOR kinase inhibitor AZD8055 on acquired endocrine resistance in breast cancer in vitro. Breast Cancer Res (2014) 16(1):R12. doi: 10.1186/bcr3604

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Xu DQ, Toyoda H, Yuan XJ, Qi L, Chelakkot VS, Morimoto M, et al. Anti-tumor effect of AZD8055 against neuroblastoma cells in vitro and in vivo. Exp Cell Res (2018) 365(2):177–84. doi: 10.1016/j.yexcr.2018.02.032

PubMed Abstract | CrossRef Full Text | Google Scholar

73. El-Naggar S, Liu Y, Dean DC. Mutation of the Rb1 pathway leads to overexpression of mTor, constitutive phosphorylation of akt on serine 473, resistance to anoikis, and a block in c-raf activation. Mol Cell Biol (2009) 29(21):5710–7. doi: 10.1128/mcb.00197-09

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Zhao X, Li C, Chiu MC, Qiao R, Jiang S, Wang P, et al. Rock1 is a novel host dependency factor of human enterovirus A71: Implication as a drug target. J Med Virol (2022) 94(11):5415–24. doi: 10.1002/jmv.27975

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Miñambres R, Guasch RM, Perez-Aragó A, Guerri C. The RhoA/ROCK-I/MLC pathway is involved in the ethanol-induced apoptosis by anoikis in astrocytes. J Cell Sci (2006) 119(Pt 2):271–82. doi: 10.1242/jcs.02723

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Schackmann RC, van Amersfoort M, Haarhuis JH, Vlug EJ, Halim VA, Roodhart JM, et al. Cytosolic p120-catenin regulates growth of metastatic lobular carcinoma through Rock1-mediated anoikis resistance. J Clin Invest (2011) 121(8):3176–88. doi: 10.1172/jci41695

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Browne AJ, Göbel A, Thiele S, Hofbauer LC, Rauner M, Rachner TD. p38 MAPK regulates the wnt inhibitor dickkopf-1 in osteotropic prostate cancer cells. Cell Death Dis (2016) 7(2):e2119. doi: 10.1038/cddis.2016.32

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Buchheit CL, Angarola BL, Steiner A, Weigel KJ, Schafer ZT. Anoikis evasion in inflammatory breast cancer cells is mediated by bim-EL sequestration. Cell Death Differ (2015) 22(8):1275–86. doi: 10.1038/cdd.2014.209

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Tang W, Feng X, Zhang S, Ren Z, Liu Y, Yang B, et al. Caveolin-1 confers resistance of hepatoma cells to anoikis by activating IGF-1 pathway. Cell Physiol Biochem (2015) 36(3):1223–36. doi: 10.1159/000430292

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Li H, Baldwin BR, Zahnow CA. LIP expression is regulated by IGF-1R signaling and participates in suppression of anoikis. Mol Cancer (2011) 10:100. doi: 10.1186/1476-4598-10-100

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Hou X, Huang F, Macedo LF, Harrington SC, Reeves KA, Greer A, et al. Dual IGF-1R/InsR inhibitor BMS-754807 synergizes with hormonal agents in treatment of estrogen-dependent breast cancer. Cancer Res (2011) 71(24):7597–607. doi: 10.1158/0008-5472.Can-11-1080

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Zillhardt M, Park SM, Romero IL, Sawada K, Montag A, Krausz T, et al. Foretinib (GSK1363089), an orally available multikinase inhibitor of c-met and VEGFR-2, blocks proliferation, induces anoikis, and impairs ovarian cancer metastasis. Clin Cancer Res (2011) 17(12):4042–51. doi: 10.1158/1078-0432.Ccr-10-3387

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Patson B, BC R, Olszanski AJ. Pharmacokinetic evaluation of axitinib. Expert Opin Drug Metab Toxicol (2012) 8(2):259–70. doi: 10.1517/17425255.2012.652947

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: neuroblastoma, anoikis, tumor microenvironment, immune cell, WGCNA, drug sensitivity

Citation: Chen J, Sun M, Chen C, Kang M, Qian B, Sun J, Ma X, Zhou J, Huang L, Jiang B and Fang Y (2023) Construction of a novel anoikis-related prognostic model and analysis of its correlation with infiltration of immune cells in neuroblastoma. Front. Immunol. 14:1135617. doi: 10.3389/fimmu.2023.1135617

Received: 01 January 2023; Accepted: 23 March 2023;
Published: 04 April 2023.

Edited by:

Hai Fang, Shanghai Jiao Tong University, China

Reviewed by:

Timofey Dmitrievich Lebedev, Engelhardt Institute of Molecular Biology (RAS), Russia
Hongyan Guo, Louisiana State University Health Science Center Shreveport, United States
Jingjing Song, Second Affiliated Hospital and Yuying Children’s Hospital of Wenzhou Medical University, China

Copyright © 2023 Chen, Sun, Chen, Kang, Qian, Sun, Ma, Zhou, Huang, Jiang and Fang. 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: Lei Huang, surgeonhuang@126.com; Bin Jiang, jiangbin363@163.com; Yongjun Fang, fyj322@189.cn

These authors have contributed equally to this work and share the first authorship

Disclaimer: 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.