Mingyuan Luan1, Xue Tian2, Dexiang Zhang3, Xiaoning Sun2, Minglu Jiang4, Yunbo Duan5, Changgang Sun6, Hongzong Si5,7. 1. Organ Transplantation Center, the Affiliated Hospital of Qingdao University, Qingdao, China. 2. School of Basic Medicine, Qingdao University Medical College, Qingdao, China. 3. Qingdao University Intelligent Campus and Information Construction Center, Qingdao University, Qingdao, China. 4. Binzhou Medical University, Binzhou, China. 5. Institute for Computational Science and Engineering, Laboratory of New Fibrous Materials and Modern Textile State Key Laboratory, Qingdao University, Qingdao, China. 6. Department of Cancer Center, Weifang Traditional Chinese Medicine Hospital, Weifang, China. 7. Department of Public Health, Qingdao University Medical College, Qingdao, China.
Abstract
BACKGROUND: Neutrophils play a crucial role in the development and progression of hepatocellular carcinoma (HCC); however, the mechanism underlying neutrophil recruitment is not fully understood. Therefore, we aimed to explore the potential genes or pathways related to neutrophil recruitment in the cancer microenvironment. METHODS: We downloaded TCGA HCC gene expression profiles, the abundance of 22 different immune cells in HCC patients, and patient survival information. We used Kaplan-Meier survival analysis to determine if neutrophils were related to survival. Next, we screened different expression genes (DEGs) between patients with high and low level of neutrophils. We then identified the transcription factor and its targets in the fence of DEGs. Then, we carried out enrichment analysis and gene set variation analysis (GSVA) for targets. Finally, we explored the potential mechanism of targets via calculating correlation scores. RESULTS: Our survival analysis results showed that neutrophils were significantly associated with patient survival. A total of 736 DEGs were screened. Next, we identified transcription factor larger E26 transformation-specific (ETS) homologous factor (EHF) and 702 targets of EHF from 736 DEGs. Among these targets, the level of FGD6 expression had the highest correlation with the level of EHF expression. Enrichment and GSVA analysis for FGD6 showed that the level of GO:0043547 had a positive regulatory effect on GTPase activity and the GO:0007010 cytoskeleton organization was significantly difference between the high and low neutrophils counts. By calculating the correlation between FGD6 and genes in GO:0043547 and GO:0007010, we identified RIC8B and SIPA1L3. CONCLUSIONS: These findings demonstrated that transcription factor EHF can influence recruitment of neutrophils by mediating the transcription of FGD6. Further investigations are needed to shed new light on EHF and its target FGD6. 2021 Translational Cancer Research. All rights reserved.
BACKGROUND: Neutrophils play a crucial role in the development and progression of hepatocellular carcinoma (HCC); however, the mechanism underlying neutrophil recruitment is not fully understood. Therefore, we aimed to explore the potential genes or pathways related to neutrophil recruitment in the cancer microenvironment. METHODS: We downloaded TCGA HCC gene expression profiles, the abundance of 22 different immune cells in HCC patients, and patient survival information. We used Kaplan-Meier survival analysis to determine if neutrophils were related to survival. Next, we screened different expression genes (DEGs) between patients with high and low level of neutrophils. We then identified the transcription factor and its targets in the fence of DEGs. Then, we carried out enrichment analysis and gene set variation analysis (GSVA) for targets. Finally, we explored the potential mechanism of targets via calculating correlation scores. RESULTS: Our survival analysis results showed that neutrophils were significantly associated with patient survival. A total of 736 DEGs were screened. Next, we identified transcription factor larger E26 transformation-specific (ETS) homologous factor (EHF) and 702 targets of EHF from 736 DEGs. Among these targets, the level of FGD6 expression had the highest correlation with the level of EHF expression. Enrichment and GSVA analysis for FGD6 showed that the level of GO:0043547 had a positive regulatory effect on GTPase activity and the GO:0007010 cytoskeleton organization was significantly difference between the high and low neutrophils counts. By calculating the correlation between FGD6 and genes in GO:0043547 and GO:0007010, we identified RIC8B and SIPA1L3. CONCLUSIONS: These findings demonstrated that transcription factor EHF can influence recruitment of neutrophils by mediating the transcription of FGD6. Further investigations are needed to shed new light on EHF and its target FGD6. 2021 Translational Cancer Research. All rights reserved.
Currently, hepatocellular carcinoma (HCC) is the most frequent liver malignancy worldwide, and the morbidity and mortality rates are high (1). Within 5 years after curative resection, metastasis and recurrence of HCC are up to 70% (2). The tumor immune microenvironment plays a crucial role in initiation, progression, metastasis, and recurrence of HCC. As major components in the tumor immune microenvironment, neutrophils have been increasingly recognized as important factors in cancer behavior and progression (3). Recent evidence has demonstrated that neutrophils contribute to the metastatic process of tumor cells via cell-to-cell interactions (4). Evidence suggests that the neutrophil-to-lymphocyte ratio can predict the prognosis of HCC (5). Therefore, a deep understanding of the mechanisms leading to the interaction between neutrophils and cancer cells will provide better strategies to orchestrate the immune system to eradicate cancers.Recently, a number of studies have focused on the roles of neutrophils in the tumor microenvironment. Neutrophils have been confirmed to induce the stem cell characteristics in HCC cells by secreting BMP2 and TGFB2 (6). Xiao et al. (7) reported that neurotensin increases the neutrophil count of the local HCC microenvironment by inducing the expression of IL-8. Haider et al. (8) concluded that the CXCL5-dependent attraction of neutrophils depends on TGF-β and Axl signaling. Despite substantial progress in research involving the roles of neutrophils in the tumor microenvironment, the precise molecular mechanisms underlying neutrophil recruitment are only partially understood.In this study we showed that the transcription factor, EHF, and its target, FGD6, had significantly different levels of expression in patients with high- and low-level neutrophil counts. The complete nomenclature of transcription factor EHF is larger E26 transformation-specific (ETS) homologous factor. Increasing evidence over recent years has revealed that EHF has a critical function in tumorigenesis. It has been demonstrated that EHF activates the TGF-β signaling pathway in colorectal carcinoma cells (9). Increased expression of EHF leads to poor survival in gastric cancer by activating HER family signaling (10). Liu et al. (11) reported that the level of EHF expression has a close relationship with the efficacy of anti-PD-1 therapy in pancreatic ductal adenocarcinoma. Our outcome showed that the expression level of FGD6 had significant correlation with the expression level of RIC8B and SIPA1L3. Previous evidence demonstrates that FGD6 can interact with CDC42 which contributes to cell adhesion, cell polarity, and membrane recycling during bone degradation (12). Taken together, our data suggest that the transcription factor, EHF, influenced recruitment of neutrophils by mediating the transcription of FGD6. EHF-FGD6-CDC42 plays an important role in neutrophil recruitment.We present the following article in accordance with the MDAR reporting checklist (available at http://dx.doi.org/10.21037/tcr-20-2714).
Methods
HCC data sets downloading and preprocessing
RNA-seq data and clinical data of HCC patient samples were extracted from the National Cancer Institute Genomic Data Commons (https://gdc.cancer.gov/about-data/publications/pancanatlas). We performed log2 normalization for RNA-seq data. The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).
Assess and analyze the infiltration of immune cells in HCC
Infiltrations of immune cells in HCC were assessed using the CIBERSORT algorithm (13) with default parameters. Next, we carried out Kaplan-Meier survival analysis to screen the immune cells that were significantly correlated with patient survival. Patients were divided into high and low neutrophil groups based on the median neutrophil count. Kaplan-Meier survival analysis was implemented for the neutrophil fraction using the survival package in R (14). The P values for the Kaplan-Meier survival analysis were adjusted using the false discovery rate (FDR) method. We screened the different expression genes (DEGs) in the two groups utilizing the limma package in R (15) [filter criteria: fold change (FC) ≥1.5; adjusted P value <0.01]. The P values were further adjusted using the Benjamini-Hochberg (BH) multiple testing correction method.
Construct prognosis prediction model using least absolute shrinkage and selection operator (LASSO) regression
In this step, we constructed a prognostic predictive module using LASSO regression. We screened survival-related genes using univariate Cox proportional risk regression model in survival R package (14). P<0.001 was selected as the threshold. We screened the most powerful prognostic markers from the survival-related genes using LASSO regression in the glmnet R package (16). The workflow of the prognostic predictive module was as follows: (I) patients were randomly divided into training (70% patients) and test sets (30% patients); (II) the training set was trained on the prognostic model and the performance of the model on the test set was assessed; (III) and the concordance index (c index) of the prognostic model in the training and test sets were calculated utilizing the Hmisc R package (17). We used a univariate Cox proportional risk regression model for age, gender, stage, stage T, stage N, stage M, cigarette smoking, and subtypes. The survival-related clinical features were selected for multivariate Cox proportional risk regression for risk scores and abundance of immune cells.
Identifying transcription factor and the targets
The detailed information of human transcription factor was downloaded from Ensembl (18). We screened all transcription factors from DEGs based on information about transcription factors downloaded from Ensembl. Next, we selected the transcription factor with the maximum logFC for further analysis. The ChIP-seq outcome of transcription factor was downloaded from Cistrome Data Browser (19). Yan et al. (20) contributed to the ChIP-seq outcome. We selected the genes which were at the intersection of the ChIP-seq outcome and DEGs as the target of transcription factor for further analysis. We uploaded genes we selected to DAVID (21) and assessed the gene ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Next, we calculated the correlation between the transcription factor and the differentially-expressed targets. The target that had the greatest correlation was selected for further analysis.
Function analysis of transcription factor targets
We selected the target that had the greatest correlation with the transcription factor. Then, we submitted the target to DAVID and obtained the GO functions and KEGG pathways. The GSVA package in R (22) was used to calculate the GO function and KEGG pathway scores in all patients. The limma package in R to screen GO functions and KEGG pathways were shown to be different between the high and low neutrophil group (filter criteria: log2FC ≥0.1; adjusted P value <0.01). To further evaluate the possible mechanisms underlying the target with the greatest correlation, we calculated the correlation between the target with the greatest correlation and the screened genes with different GO functions and KEGG pathways. The functions and pathways of significantly correlated genes were annotated by DAVID.
Statistical analysis
Wilcoxon test was utilized to compare continuous variables. Correlations were calculated with Spearman correlation. All the tests were two-sided, and P<0.05 was regarded as indicating significance, unless otherwise stated. Survival analysis was performed using the Kaplan-Meier method, and the survival of the clusters was compared using the log-rank test. The FDR correction was utilized in multiple tests to decrease false positive rates. All of the analyses were performed with R software (version 3.4.2, http://www.R-project.org).
Results
HCC data sets downloading and pre-processing
RNA-seq data and clinical data of 323 HCC patients were extracted from the National Cancer Institute Genomic Data Commons (https://gdc.cancer.gov/about-data/publications/pancanatlas). We performed log2 normalization for RNA-seq data.
Assessment and analysis of immune cell infiltration in HCC
Utilizing the CIBERSORT algorithm, we assessed infiltration of 22 immune cells in HCC (). Next, we divided patients into high and low immune cell groups based on the median of each immune cell count. We carried out Kaplan-Meier survival analysis for 22 immune cells. The P values of Kaplan-Meier survival analysis were adjusted using the FDR method (). Our outcome showed that neutrophil infiltration was significantly associated with patient survival (). Next, we screened DEGs between the high (n=147) and low neutrophil groups (n=206) using the limma R package. A total of 736 DEGs were screened between the high and low neutrophil groups (). Five hundred sixty-three DEGs had significantly high expression, whereas 173 DEGs were downregulated in the high neutrophil group.
Figure 1
Infiltration of immune cells in HCC. (A) The abundance of 22 immune cells assessed by CIBERSORT. (B) The Kaplan-Meier survival analysis of neutrophils. The infiltration of neutrophils is significantly correlated with HCC patients’ survival. HCC, hepatocellular carcinoma.
Table 1
The outcome of Kaplan-Meier survival analysis
Immune cells
FDR
T cells (CD4 memory activated)
0.988706859
T cells (CD4 naïve)
0.952450144
NK cells (resting)
0.952450144
T cells (follicular helper)
0.94049182
Eosinophils
0.915125939
Dendritic cells (resting)
0.915125939
T cells (CD8)
0.915125939
B cells (memory)
0.748294157
Mast cells (activated)
0.748294157
Macrophages (M1)
0.748294157
Mast cells (resting)
0.664405339
Plasma cells
0.576581743
Dendritic cells (activated)
0.576581743
T cells regulatory (Tregs)
0.576581743
T cells (gamma delta)
0.576581743
Monocytes
0.525440109
NK cells (activated)
0.525440109
T cells (CD4 memory resting)
0.12998342
B cells (naïve)
0.12998342
Macrophages (M2)
0.042122642
Macrophages (M0)
0.042122642
Neutrophils
0.042122642
P values were adjusted using FDR method. FDR, false discovery rate.
Figure 2
Total 736 DEGs were screened between high and low neutrophils group. Red points represent up-regulated DEGs while blue points represent down-regulated DEGs. DEGs, different expression genes.
Infiltration of immune cells in HCC. (A) The abundance of 22 immune cells assessed by CIBERSORT. (B) The Kaplan-Meier survival analysis of neutrophils. The infiltration of neutrophils is significantly correlated with HCC patients’ survival. HCC, hepatocellular carcinoma.P values were adjusted using FDR method. FDR, false discovery rate.Total 736 DEGs were screened between high and low neutrophils group. Red points represent up-regulated DEGs while blue points represent down-regulated DEGs. DEGs, different expression genes.
Construction of a prognosis predictive model using LASSO regression
A total of 323 patients were selected to fit the prognosis predictive model. We screened 1,616 survival-related genes using a univariate Cox proportional risk regression model with a P<0.001 threshold. Next, we performed LASSO regression to screen survival-related markers. A total of three survival-related markers were screened. The risk score for each patient was calculated by combining the expression of the prognosis-related markers with the corresponding LASSO coefficients. The formula of the risk score is listed as follows:To assess the performance of our model, we plotted cross-validated time-dependent receiver operating characteristic (ROC) curves. The 12-month area under the curve (AUC) was 0.823 (). The 36- and 60-month AUCs were both >0.79. The AUCs indicated that our model had good performance in prognosis prediction both in the short- and long-term. We calculated the optimal cut-off value according to the optimum sensitivity and specificity of the 5-year survival ROC curve. Patients with a risk score ≥−0.027 were assigned to the high-risk group, and the remaining patients comprised the low-risk group. The survival of patients in the low-risk group was significantly better than patients in the high-risk group (). The number of patients in high risk group was less than the number of patients in low risk group (). The high-risk group had more deaths than the low-risk group (). As presented in , RAP1GAP and SLC39A1 in patients in the high-risk group exhibited high expression levels while ZNF23 showed the opposite tendency. Finally, we carried out multivariate Cox regression analysis and the outcome revealed that our prognosis predictive model independently predicted prognosis ().
Figure 3
Performance of prognosis prediction model. (A) Kaplan-Meier survival analysis of the high- and low-risk groups. Patients in the low-risk group showed significant better prognosis than patients in high-risk group. (B) The time-dependent ROC curves of prognosis prediction model. AUCs at 1, 3, and 5 years were more than 0.79 which indicated that our model had good performance in predicting prognosis. AUC, area under the curve; ROC, receiver operating characteristic; AUC, area under the curve.
Figure 4
Prognosis prediction model of HCC patients. (A) The distribution of risk scores. Our outcome showed that −0.027 was the optimal cut-off value. Patients were divided into a high-risk (blue) and low-risk groups (yellow) based on the optimal cut-off value. (B) Patient survival time and status. The black dotted line represents the optimum cut-off dividing patients into low-risk (left) and high-risk (right) groups. The yellow dots represent death and the blue dot represents alive status. The number of deaths in the high-risk group was significantly greater compared with that in the low-risk group. (C) Heatmap of the expression level of the three prognosis markers in the prognosis model. Red represents high expression level while blue represents low expression level. ZNF23 exhibited lower expression level in the high-risk group while RAP1GAP and SLC39A1 showed opposite trend. HCC, hepatocellular carcinoma.
Table 2
Prognostic significance of prognosis risk score in HCC
Variables
Univariate
Multivariate
HR (95% CI)
P value
HR (95% CI)
P value
Age (years)
1.01 (1.00–1.03)
0.1270
Not included
Gender
Male
Ref.
Ref.
Female
1.48 (1.00–2.19)
0.0496
1.38 (0.92–2.06)
0.1167
Stage
Stage 1
Ref.
Ref.
Stage 2
1.27 (0.77–2.09)
0.3540
0.88 (0.53–1.46)
0.6170
Stage 3
2.55 (1.64–3.97)
2.98×10−5
1.78 (1.13–2.81)
0.0124
Stage 4
4.49 (1.38–14.61)
0.0124
3.45 (1.04–11.43)
0.0426
Risk score
4.38 (3.25–5.90)
3.17×10−22
4.24 (3.11–5.77)
5.46×10−20
HRs and P values of risk score and the covariates in the univariate and multivariate Cox proportional hazards model. Risk score can predict prognosis independently. HCC, hepatocellular carcinoma; HR, hazard ratio; CI, confidence interval.
Performance of prognosis prediction model. (A) Kaplan-Meier survival analysis of the high- and low-risk groups. Patients in the low-risk group showed significant better prognosis than patients in high-risk group. (B) The time-dependent ROC curves of prognosis prediction model. AUCs at 1, 3, and 5 years were more than 0.79 which indicated that our model had good performance in predicting prognosis. AUC, area under the curve; ROC, receiver operating characteristic; AUC, area under the curve.Prognosis prediction model of HCC patients. (A) The distribution of risk scores. Our outcome showed that −0.027 was the optimal cut-off value. Patients were divided into a high-risk (blue) and low-risk groups (yellow) based on the optimal cut-off value. (B) Patient survival time and status. The black dotted line represents the optimum cut-off dividing patients into low-risk (left) and high-risk (right) groups. The yellow dots represent death and the blue dot represents alive status. The number of deaths in the high-risk group was significantly greater compared with that in the low-risk group. (C) Heatmap of the expression level of the three prognosis markers in the prognosis model. Red represents high expression level while blue represents low expression level. ZNF23 exhibited lower expression level in the high-risk group while RAP1GAP and SLC39A1 showed opposite trend. HCC, hepatocellular carcinoma.HRs and P values of risk score and the covariates in the univariate and multivariate Cox proportional hazards model. Risk score can predict prognosis independently. HCC, hepatocellular carcinoma; HR, hazard ratio; CI, confidence interval.
Identifying the transcription factor and the targets
The information of human transcription factor was downloaded from Ensembl. Through comparison, we identified 18 transcription factors in DEGs (). Transcription factor EHF with maximum logFC (logFC =1.1717) was selected for further analysis (). We downloaded the ChIP-seq outcome of EHF from Cistrome Data Browser. Based on the ChIP-seq outcome, we screened 702 targets of transcription factor EHF in DEGs. Among these targets, the expression level of FGD6 was most correlated with transcription factor EHF (). Next, we submitted the targets of EHF to DAVID and carried out pathways and GO function analysis. The outcome revealed that the targets of EHF mainly enriched in GO:0016021 integral component of membrane, GO:0006954 inflammatory response, GO:0005887 integral component of plasma membrane, GO:0005886 plasma membrane and GO:0005615 extracellular space ().
Figure 5
Differentially expressed transcription factor between high level of neutrophils and low level of neutrophils. (A) Eighteen differentially expressed transcription factors were identified. The size of the point is correlated with absolute value of logFC. The red color represents the less adjusted P value while the green color is opposite. (B) High neutrophils patients had significant higher expression level of EHF than low neutrophils patients. (C) The expression level of FGD6 was most correlated with transcription factor EHF. Scatter plots of log2-transformed gene expression values of FGD6 and EHF are shown. Linear regression line is plotted with the gray shaded region showing the 95% confidence interval. Pearson’s correlation coefficient r and P values are given at the top. The levels of neutrophils are labeled. Median values of the expression of FGD6 and EHF are indicated with dashed gray lines. We carried out log-rank statistics to identify the optimal cut-off for transforming the continuous variable of gene expression into categorical high- and low-expression groups in a survfit model. The test score at each candidate cut-off across the log-transformed gene expression values was plotted. The highest test score (indicated with a blue arrow) was applied for best separating patients into 4 different risk groups (using solid blue lines; named groups I to IV). FC, fold change; EHF, larger E26 transformation-specific homologous factor.
Figure 6
GO functions and KEGG pathways which the targets of EHF significantly enriched in. The outcome revealed that the targets of EHF mainly enriched in GO:0016021 integral component of membrane, GO:0006954 inflammatory response, GO:0005887 integral component of plasma membrane, GO:0005886 plasma membrane and GO:0005615 extracellular space. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; EHF, larger E26 transformation-specific homologous factor.
Differentially expressed transcription factor between high level of neutrophils and low level of neutrophils. (A) Eighteen differentially expressed transcription factors were identified. The size of the point is correlated with absolute value of logFC. The red color represents the less adjusted P value while the green color is opposite. (B) High neutrophils patients had significant higher expression level of EHF than low neutrophils patients. (C) The expression level of FGD6 was most correlated with transcription factor EHF. Scatter plots of log2-transformed gene expression values of FGD6 and EHF are shown. Linear regression line is plotted with the gray shaded region showing the 95% confidence interval. Pearson’s correlation coefficient r and P values are given at the top. The levels of neutrophils are labeled. Median values of the expression of FGD6 and EHF are indicated with dashed gray lines. We carried out log-rank statistics to identify the optimal cut-off for transforming the continuous variable of gene expression into categorical high- and low-expression groups in a survfit model. The test score at each candidate cut-off across the log-transformed gene expression values was plotted. The highest test score (indicated with a blue arrow) was applied for best separating patients into 4 different risk groups (using solid blue lines; named groups I to IV). FC, fold change; EHF, larger E26 transformation-specific homologous factor.GO functions and KEGG pathways which the targets of EHF significantly enriched in. The outcome revealed that the targets of EHF mainly enriched in GO:0016021 integral component of membrane, GO:0006954 inflammatory response, GO:0005887 integral component of plasma membrane, GO:0005886 plasma membrane and GO:0005615 extracellular space. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; EHF, larger E26 transformation-specific homologous factor.We submitted FGD6 to DAVID and obtained 16 GO function terms and zero KEGG pathways. GSVA and limma analysis demonstrated that the high and low neutrophil groups had significant differences in 2 GO terms [GO:0043547 (positive regulation of GTPase activity) and GO:0007010 (cytoskeleton organization); ]. To further evaluate the possible mechanisms underlying FGD6 in GO:0043547 and GO:0007010, we calculated the correlation between FGD6 and genes in the two GO terms. The outcome showed that the RIC8B and SIPA1L3 genes had the greatest correlation (). The KEGG pathways analysis utilizing DAVID showed that SIPA1L3 was involved in the hsa04015: Rap1 signaling pathway.
Figure 7
GO functions of FGD6. FGD6 enriched in 16 GO functions. Limma analysis showed the high neutrophils group and low neutrophils group had significant differences in two GO terms: GO:0043547 positive regulation of GTPase activity and GO:0007010 cytoskeleton organization. GO, gene ontology.
Figure 8
The correlation scatter plots of FGD6 with RIC8B and SIPA1L3. (A) The correlation scatter plots of FGD6-RIC8B. (B) The correlation scatter plots of FGD6-SIPA1L3. Scatter plots of log2-transformed gene expression values of FGD6, RIC8B and SIPA1L3 are shown. Linear regression lines are plotted with the gray shaded region showing the 95% confidence interval. Pearson’s correlation coefficient r and P values are given at the top. The levels of neutrophils are labeled. Median values of the expression of FGD6, RIC8B and SIPA1L3 are indicated with dashed gray lines. We carried out log-rank statistics to identify the optimal cut-off for transforming the continuous variable of gene expression into categorical high- and low-expression groups in a survfit model. The test score at each candidate cut-off across the log-transformed gene expression values was plotted. The highest test score (indicated with a blue arrow) was applied for best separating patients into 4 different risk groups (using solid blue lines; named groups I to IV).
GO functions of FGD6. FGD6 enriched in 16 GO functions. Limma analysis showed the high neutrophils group and low neutrophils group had significant differences in two GO terms: GO:0043547 positive regulation of GTPase activity and GO:0007010 cytoskeleton organization. GO, gene ontology.The correlation scatter plots of FGD6 with RIC8B and SIPA1L3. (A) The correlation scatter plots of FGD6-RIC8B. (B) The correlation scatter plots of FGD6-SIPA1L3. Scatter plots of log2-transformed gene expression values of FGD6, RIC8B and SIPA1L3 are shown. Linear regression lines are plotted with the gray shaded region showing the 95% confidence interval. Pearson’s correlation coefficient r and P values are given at the top. The levels of neutrophils are labeled. Median values of the expression of FGD6, RIC8B and SIPA1L3 are indicated with dashed gray lines. We carried out log-rank statistics to identify the optimal cut-off for transforming the continuous variable of gene expression into categorical high- and low-expression groups in a survfit model. The test score at each candidate cut-off across the log-transformed gene expression values was plotted. The highest test score (indicated with a blue arrow) was applied for best separating patients into 4 different risk groups (using solid blue lines; named groups I to IV).
Discussion
Our findings demonstrated that the neutrophil count was significantly correlated with patient survival. The transcription factor, EHF, and its target, FGD6, had significantly different levels of expression between patients with high and low neutrophil counts. GO:0043547 exerts positive regulation of GTPase activity and GO:0007010 facilitates cytoskeleton organization; FGD6 participates in both functions and there were significant differences in patients with high and low neutrophil counts. Correlation analysis showed that FGD6 had the greatest correlation with RIC8B and SIPA1L3 in the positive regulation of GTPase activity and cytoskeleton organization. SIPA1L3 is a member of the hsa04015: Rap1 signaling pathway.EHF is a member of the highly diverse ETS superfamily. Accumulation evidence over recent years has demonstrated that EHF has a great effect on tumorigenesis and the tumor microarray environment. Studies have revealed that EHF inhibits tumor invasion, metastasis, and the mesenchymal phenotype (23-25). Liu et al. (11) reported that deficiency of EHF induces T reg cell and myeloid derived suppressor cell (MDSC) accumulation and causes a decrease in CD8+ T cell infiltration. Wang et al. (9) found that EHF promotes colorectal carcinoma progression via activation pf TGF-β1 transcription and canonical TGF-β signaling. Liu et al. (11) revealed that overexpression of EHF can significantly lead to better response to anti-PD-1 therapy. High expression of EHF promotes oral squamous cell carcinoma metastasis, cancer stemness, and drug resistance (26). Taken together, EHF has emerged as a promising target in cancer research; however, the contributions of EHF to neutrophils in HCC are not fully understood. We found that the level of EHF expression was different between patients with high and low neutrophil counts. Our finding suggested that EHF plays an important role in neutrophil recruitment.To understand the mechanism underlying EHF, we screened its downstream target, FGD6, based on the outcome of ChIP-seq. FGD6 is a member of the faciogenital dysplasia family. FGD6 includes six Rho guanine nucleotide exchange factors (GEFs) that activate Rho proteins (12). FGD6, which is highly homologous to FGD1, has two pleckstrin homology (PH) domains and one Dbl homology (DH) domain. Previous studies have shown that mutations in the FGD1 gene are linked to the inherited disease, faciogenital dysplasia, and Aarskog-Scott syndrome (27,28). Recent evidence demonstrated that FGD6 plays a crucial role in cell polarity and membrane recycling by regulating the assembly of different actin-based protein networks and activating CDC42 at different locations (12). Our data showed that FGD6 is involved in GO:0043547 positive regulation of GTPase activity and GO:0007010 cytoskeleton organization. In addition, the level of FGD6 expression was highly correlated with the level of RIC8B and SIPA1L3 expression in GO:0043547 positive regulation of GTPase activity and GO:0007010 cytoskeleton organization. This finding indicated that FGD6, RIC8B, and SIPA1L3 may work together in the positive regulation of GTPase activity and cytoskeleton organization.SIPA1L3, which has N-terminal RapGAP (a Rap-GTPase activating domain), PDZ, and C-terminal leucine zipper domains, is a member of the signal-induced proliferation-associated (SIPA) proteins. A previous study revealed that loss of SIPA1L3 interferes with cell polarity and cytoskeletal organization (29). Previous evidence showed that RIC8B is GEFs for the α subunits of heterotrimeric guanine nucleotide-binding proteins (30). In addition, it has been suggested that RIC8B contributes to the steady state of G protein α subunits (31). Some studies have showed that overexpression of RIC8B enhances activation of adenylyl cyclase (AC) via Gs and Golf (32,33). We discovered that the level of FGD6 expression had a high correlation with the level of RIC8B and SIPA1L3 expression. Therefore, it is essential to explore the relationship and potential mechanisms underlying FGD6, RIC8B, and SIPA1L3.Our data showed that SIPA1L3 participates in the has04015: Rap1 signaling pathway. Accumulating evidence demonstrates that the Rap1 signaling pathway contributes a major role in the endothelial barrier regulation, vasculogenesis, and angiogenesis (34,35). A recent study involving the Rap1 signaling pathway has focused on integrin-mediated cell adhesion (36). Based on the diagram of the has04015: Rap1 signaling pathway on the KEGG website (https://www.genome.jp/kegg-bin/show_pathway?map=hsa04015&show_description=show) (), we found that the G protein-coupled receptors (GPCRs) signaling pathway activates the Rap1 signaling pathway via EPAC. EPAC, which includes GEFs, can exchange bound GDP of Rap1 for free GTP. In contrast, SIPA1L3, which contains the RapGAP domain, can activate the Rap-GTPase inhibiting Rap1 signaling pathway. Active Rap1 can activate CDC42 via Src and FRG. CDC42 will activate the process of cell adhesion, migration, and polarity. Collectively, these processes can be characterized in three sections: (I) the process of GPCRs regulating Rap1; (II) the process of regulation between Rap1-GDP and Rap1-GTP; and (III) the regulation of CDC42 in Rap1 downstream. We showed that FGD6, RIC8B, and SIPA1L3 are crucial regulatory factors in these three sections. RIC8B enhances activation of AC via Gs. This process will trigger EPAC to promote Rap1-GDP conversion into Rap1-GTP, which will activate the Rap1 signaling pathway. SIPA1L3, which has an N-terminal RapGAP domain, promotes the transformation of Rap1-GTP to Rap1-GDP. This process inhibits the Rap1 signaling pathway. FGD6, which is regulated by the transcription factor, EHF, can activate CDC42. This process will positively regulate the Rap1 signaling pathway and contribute to cell adhesion, migration, and polarity.
Figure 9
The diagram of has04015: Rap1 signaling pathway. The pertinent proteins have been marked by red color on the plot. RIC8B can enhance activation of adenylyl cyclase (AC) via Gs. This process will trigger EPAC to promote Rap1-GDP converting into Rap1-GTP. SIPA1L3, which has an N-terminal Rap-GAP domain (Rap-GTPase activating domain), can promote the transforming of Rap1-GTP to Rap1-GDP. This process will inhibit Rap1 signaling pathway. FGD6 which includes six Rho GEFs activating Rho proteins can activate CDC42.
The diagram of has04015: Rap1 signaling pathway. The pertinent proteins have been marked by red color on the plot. RIC8B can enhance activation of adenylyl cyclase (AC) via Gs. This process will trigger EPAC to promote Rap1-GDP converting into Rap1-GTP. SIPA1L3, which has an N-terminal Rap-GAP domain (Rap-GTPase activating domain), can promote the transforming of Rap1-GTP to Rap1-GDP. This process will inhibit Rap1 signaling pathway. FGD6 which includes six Rho GEFs activating Rho proteins can activate CDC42.We inferred that the transcription factor, EHF, mediates GO:0043547 positive regulation of GTPase activity, GO:0007010 cytoskeleton organization, and the has04015: Rap1 signaling pathway by regulating the transcription of FGD6. These processes will mediate local adhesion. Increasing local adhesion may contribute to the recruitment and persistence of neutrophils in tumor lesions. Also, tumor cells might be protected and transported by the surrounding neutrophils following cell-cell interactions (3). As the key step in tumor metastasis, increasing adhesion will raise the ability of metastasis and colonization of tumor cells. This feature can contribute to the attachment of tumor cells to matrix collagen fibers of vascular endothelial cells and secondary lesions (37). These reasons suggest that patients with a high neutrophil count had worse survival. Increased EHF upregulates FGD6. Increased FGD6 and RIC8B or decreased SIPA1L3 will contribute to activation of the Rap1 signaling pathway. Activation of the Rap1 signaling pathway will increase cell adhesion, which is advantageous to neutrophil recruitment and metastasis, and colonization of tumor cells.
Conclusions
These findings demonstrated that transcription factor EHF can influence recruitment of neutrophils by mediating the transcription of FGD6. Further investigations are needed to shed new light on EHF and its target FGD6 which may contribute to immunotherapy in HCC.
Authors: Jian Yan; Martin Enge; Thomas Whitington; Kashyap Dave; Jianping Liu; Inderpreet Sur; Bernhard Schmierer; Arttu Jolma; Teemu Kivioja; Minna Taipale; Jussi Taipale Journal: Cell Date: 2013-08-15 Impact factor: 41.582
Authors: Meital Gabay; Mary E Pinter; Forrest A Wright; PuiYee Chan; Andrew J Murphy; David M Valenzuela; George D Yancopoulos; Gregory G Tall Journal: Sci Signal Date: 2011-11-22 Impact factor: 8.192
Authors: Matthew E Ritchie; Belinda Phipson; Di Wu; Yifang Hu; Charity W Law; Wei Shi; Gordon K Smyth Journal: Nucleic Acids Res Date: 2015-01-20 Impact factor: 16.971
Authors: Domenico Albino; Gianluca Civenni; Cecilia Dallavalle; Martina Roos; Hartmut Jahns; Laura Curti; Simona Rossi; Sandra Pinton; Gioacchino D'Ambrosio; Fausto Sessa; Jonathan Hall; Carlo V Catapano; Giuseppina M Carbone Journal: Cancer Res Date: 2016-05-02 Impact factor: 12.701
Authors: Anne-Marie Benoliel; Nicolas Pirro; Valérie Marin; Bernard Consentino; Anne Pierres; Joana Vitte; Pierre Bongrand; Igor Sielezneff; Bernard Sastre Journal: Anticancer Res Date: 2003 Nov-Dec Impact factor: 2.480
Authors: Daniel R Zerbino; Nathan Johnson; Thomas Juetteman; Dan Sheppard; Steven P Wilder; Ilias Lavidas; Michael Nuhn; Emily Perry; Quentin Raffaillac-Desfosses; Daniel Sobral; Damian Keefe; Stefan Gräf; Ikhlak Ahmed; Rhoda Kinsella; Bethan Pritchard; Simon Brent; Ridwan Amode; Anne Parker; Steven Trevanion; Ewan Birney; Ian Dunham; Paul Flicek Journal: Database (Oxford) Date: 2016-02-17 Impact factor: 3.451
Authors: Christine Haider; Julia Hnat; Roland Wagner; Heidemarie Huber; Gerald Timelthaler; Markus Grubinger; Cédric Coulouarn; Wolfgang Schreiner; Karin Schlangen; Wolfgang Sieghart; Markus Peck-Radosavljevic; Wolfgang Mikulits Journal: Hepatology Date: 2018-12-20 Impact factor: 17.425
Authors: Cho-Hao Lin; Jimmy Chun-Tien Kuo; Ding Li; Aaron B Koenig; Alexander Pan; Pearlly Yan; Xue-Feng Bai; Robert J Lee; Kalpana Ghoshal Journal: Front Cell Dev Biol Date: 2022-03-24