Wen-Jie Wang1†
Han Wang2†
Ting-Yan Hua3†
Wei Song4
Jie Zhu5
Jing-Jing Wang6
Yue-Qing Huang7*
Zhi-Liang Ding8*- 1Department of Radio-Oncology, The Affiliated Suzhou Hospital of Nanjing Medical University, Suzhou, China
- 2Department of Oncology, Jining Cancer Hospital, Jining, China
- 3Department of Gastroenterology, The Affiliated Suzhou Hospital of Nanjing Medical University, Suzhou, China
- 4Department of Gastrointestinal Surgery II, Renmin Hospital of Wuhan University, Wuhan, China
- 5Department of Oncology, Changzhou Traditional Chinese Medical Hospital, Changzhou, China
- 6Department of Oncology, Taizhou Hospital of Traditional Chinese Medicine, Taizhou, China
- 7Department of General Medicine, The Affiliated Suzhou Hospital of Nanjing Medical University, Suzhou, China
- 8Department of Neurosurgery, The Affiliated Suzhou Hospital of Nanjing Medical University, Suzhou, China
Hepatocellular carcinoma (HCC) is one of the most prevalent neoplasms worldwide, particularly in China. Immune-related genes (IRGs) and immune infiltrating lymphocytes play specific roles in tumor growth. Considering how important immunotherapy has become for HCC treatment in the past decade, our objective was to establish a prognostic model by screening survival-related IRGs in patients with HCC. Using edgeR, we identified differentially expressed IRGs (DEIRGs), DEmiRNAs, and DElncRNAs. Functional enrichment analysis of DEIRGs was performed to investigate the biological functions of IRGs via gene ontology annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses. Protein-protein interaction and competing endogenous RNA networks were established using Cytoscape. Survival-associated IRGs were selected via univariate COX regression analysis, a The Cancer Genome Atlas (TCGA) prognostic model and GSE76427 validation model were developed using multivariate COX regression analysis test by AIC (Akaike Information Criterion). We identified 116 DEIRGs in patients with HCC; the “cytokine-cytokine receptor interaction” pathway was found to be the most enriched pathway. Via the prognostic model helped us classify patients into high- and low-risk score groups based on overall survival (OS); high risk score was associated with worse OS, and a positive correlation was observed between the prognostic model and immune cell infiltration. To summarize, we established a prognostic model using survival-related IRGs that provides sufficient information for prognosis prediction and immunotherapy of patients with HCC.
Introduction
Primary liver cancer is one of the most common and lethal cancers in China, with a 5-year survival rate of <5%. It can be classified into hepatocellular carcinoma (HCC) and cholangiocarcinoma (Chen et al., 2016; Bray et al., 2018). Recent studies have reported that the incidence and mortality of liver cancer have declined in China, but considering the large population, a large number of new cases are still reported each year (Chen et al., 2016).
HCC is usually associated with chronic hepatitis and cirrhosis attributable to hepatitis B or hepatitis C virus (HBV and HCV, respectively) infections (Marrero et al., 2018). HCC caused by chronic hepatitis is associated with complex mechanisms; for example, hepatitis virus infection can disrupt the immune system, leading to immune cell infiltration and cytokine secretion (Wang et al., 2013).
Patients with HCC show high mortality rates, which can be attributed to immune evasion, drug tolerance, and distant metastasis (Nagao et al., 2000; Chi et al., 2016). In the past decade, immunotherapy has emerged as a powerful modality to treat various conditions, considering that our understanding of how the immune system functions has substantially increased. Immunological checkpoint inhibitors such as pembrolizumab and nivolumab appear promising for treating patients with advanced liver cancer (El-Khoueiry et al., 2017; Zhu et al., 2018). Several studies have reported that immune cell infiltration is associated with HCC prognosis (Zhou et al., 2016; Zhou et al., 2019). Some studies have moreover described the correlation between immune-related genes (IRGs) and tumor prognosis, but few are based on HCC (Li B. et al., 2017; Lin et al., 2019).
In the present study, our objective was to establish a prognostic model by screening survival-related IRGs. We analyzed the correlation between the model and immune infiltrating cells with the aim of providing sufficient information for prognosis prediction and immunotherapy of patients with HCC.
Material and Methods
Subjects
A total of 10 pairs of HCC tissues and the corresponding adjacent normal tissues were obtained from HCC patients undergoing surgery at the Affiliated Suzhou Hospital of Nanjing Medical University (China). None of the patients had received chemotherapy or radiotherapy prior to surgery. Before surgery, written informed consent was collected from HCC patients. This work was approved by the Medical Ethics Committees of the Affiliated Suzhou Hospital of Nanjing Medical University.
Data Acquisition and Processing
RNA-sequencing and clinical data of patients with HCC were obtained from The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/) and Gene Expression Omnibus (GEO) data portal (https://www.ncbi.nlm.nih.gov/geo/). A total of 374 HCC and 50 adjacent non-cancer liver tissues were included in TCGA. GSE76427 was included 115 HCC patients. The TCGA hepatocellular carcinoma patients' clinical feathers were list in Table 1. The GSE76427 hepatocellular carcinoma patients' clinical feathers were list in Table S1. A list of IRGs was downloaded from the Immunology Database and Analysis Portal system (Bhattacharya et al., 2014).
Table 1 The clinical features of TCGA hepatocellular carcinoma patients.
Identification of Differentially Expressed Genes (DEGs), Differentially Expressed IRGs (DEIRGs), and Survival-Associated IRGs
DEGs were identified using the edgeR package in the R statistical environment (http://bioconductor.org/packages/edgeR/) (R Development Core Team, Vienna, Austria) and then further analyzed. |Log2 fold change (FC)| > 2.0 and false discovery rate (FDR) adjusted to P < 0.05 were set as the thresholds (Robinson et al., 2010). A list of immune-related genes was downloaded from the Immunology Database and Analysis Portal (ImmPort, https://www.immport.org/shared/genelists), which shares basic immunological data for cancer research. The DEIRGs were obtained by intersect the previously acquired DGEs list with the immune-related genes list. In addition, we generated volcano and heat maps of DEGs and DEIRGs using the gplots and heatmap packages in the R platform. Survival-associated DEIRGs were selected via univariate Cox regression analysis, which was performed using the survival package in the R platform.
Functional Enrichment Analysis
To understand the biological mechanisms underlying IRGs in the prognostic model, gene ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed using the DAVID (https://david.ncifcrf.gov/) online tool and ClusterProfiler, which is an R package for functional classification and enrichment of gene clusters using hypergeometric distribution (Dennis et al., 2003). The results of GO annotation and KEGG pathway analyses were visualized using the GOplot package in R. We constructed a visualized network using Cytoscape 3.6.1 (National Resource for Network Biology, USA). GO annotation and KEGG pathway analyses were based on the threshold of P < 0.05.
Constructing a Protein-Protein Interaction (PPI) Network
To clarify the potential relationships of DEmRNAs, a protein-protein interaction (PPI) network was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) 10.5 (https://string-db.org/cgi/input.pl) and visualized by Cytoscape. We then used CytoHubba to identify hub genes.
Construction of the Competing Endogenous RNA (ceRNA) Network
To further analyze the potential targets of key genes, we established a competing endogenous RNA (ceRNA) network. First, through edgeR package in the R statistical environment, we identified differentially expression miRNA (DEmiRNA) and differentially expression LncRNA (DELncRNA). |Log2 fold change (FC)| > 2.0 and FDR adjusted to P< 0.05 were set as the thresholds. Using miRcode (http://www.mircode.org/), we predict the DEmiRNA related DELncRNA in the network. Then through miRDB, miRTarBase, and TargetScan databases, we identified potential target genes for DEmiRNA. Intersecting these target genes with DEIRGs, we obtained the mRNA in the network. Based on miRNA-LncRNA and miRNA-mRNA, we established a ceRNA network, which was visualized using Cytoscape.
Development of the IRG Prognostic Model (IRGPM) and Validation Model
Survival-associated IRGs were selected via univariate Cox regression analysis using the survival package in the R platform. Patients with HCC were divided into high- and low-risk groups using the median risk score as the cut-off via multivariate Cox regression analysis test by AIC (Akaike Information Criterion). To verify the feasibility of the prognostic model, we also divided GSE76427 patients into two groups according to the median risk score. The survival of the two groups was analyzed by KM curve. The risk score calculation formula for all patients is as follows
In the formula, k represents the number of mRNA, Ci represents the coefficient of mRNA in multivariate Cox regression analysis, and Vi represents the expression level of mRNA.
Relationship Between IRGPM and Immune Cell Infiltration
We used the Tumor Immune Estimation Resource (TIMER) web server to infer the abundance of tumor infiltrating immune cells (Li T. et al., 2017). TIMER re-analyzes gene expression data and the database includes 10,897 samples across 32 cancer types from TCGA to estimate the abundance of six subtypes of tumor-infiltrating immune cells, including CD4 T cells, CD8 T cells, B cells, macrophages, dendritic cells (DCs), and neutrophils. Thus, it can be effectively used to determine the relationship between immune cell infiltration and other parameters. We downloaded data pertaining to immune cell infiltration levels among patients with HCC and assessed the association between IRGPM and immune cell infiltration.
Quantitative Real-Time PCR (qRT-PCR)
Total RNA was obtained from 10 patients with HCC using TRIzol reagent (Invitrogen, Carlsbad, USA), and then reverse transcribed with the First Strand cDNA synthesis kit (New England Biolabs, Beijing, China). We performed amplifications with a SYBR Green PCR kit (Applied Biological Materials, Richmond, Canada) according to the manufacturer's instructions on Applied Biosystems 7500Real-Time PCR System (Applied Biosystems, Foster city, USA). The expression of RNA was normalized against GAPDH using the 2-ΔΔCt method. The PCR primers used are shown in Table 2. Three separate experiments were performed.
Table 2 Primer sequences for qRT-PCR.
Statistical Analysis
Survival analysis of data related to patients in the model was performed using the survival package in the R platform. Survival curves were generated using the Kaplan-Meier method, and the log-rank test was used to compare differences between the two groups. The area under the receiver operating characteristic (ROC) curve was calculated using the survival ROC package to validate the performance of the prognostic signature (Heagerty et al., 2000). HCC and normal tissue data are reported as means±SD and were analyzed by GraphPad Prism 8.0 software (GraphPad Software Inc., USA). The significant difference between groups was evaluated using Student's t-test for a single comparison or one-way analysis of variance (ANOVA) with Bonferroni post-hoc test for multiple comparisons. Analysis items with P< 0.05 were considered statistically significant.
Results
Identification DEGs
Using edgeR v3.53, we identified 2,068 DEGs among patients with HCC; 1,991 and 77 were up- and downregulated, respectively, with the thresholds of |log2 FC| > 2.0 and adjusted P < 0.05 (Figures 1A, B). Upon further comparison with the list of IRGs, we identified 116 DEIRGs; 96 and 20 were up- and downregulated, respectively, with the same thresholds (Figures 1C, D).
Figure 1 DEGs and DEIRGs in patients with HCC: Heatmap (A) and volcano plot (B) of DEGs in HCC and adjacent non-cancer liver tissues. Heatmap (C) and volcano plot (D) of DEIRGs in HCC and adjacent non-cancer liver tissues.
Construction of the Immune-Related Genes Prognostic Model (IRGPM) and Validation Model
We screened 22 survival-associated IRGs based on the overall survival (OS) of patients using univariate Cox regression analysis (Figure 2). Via multivariate Cox regression analyses of survival-associated IRGs text by AIC, we established a prognostic model that classified patients into high- and low-risk score groups based on OS using the median risk score as the cut-off (Figure 3).The survival risk score= (Expression level of BIRC5 × 0.02296) + (Expression level of CSPG5 × 0.33178) + (Expression level of IL-11 × 0.01577) + (Expression level of FABP6 × 0.07392) + (Expression level of FIGNL2 × 0.44366) + (Expression level of GAL × 0.18222) + (Expression level of IL17D × 0.08771) + (Expression level of MAPT × 0.27133) + (Expression level of SPP1 × 0.00015) + (Expression level of STC2 × 0.02978).
Figure 2 Prognostic values of survival-associated IRGs: Forest plot of survival-associated IRGs.
Figure 3 Construction of a TCGA prognostic model using survival-associated IRGs: Heatmap (A) of survival-associated IRGs in the prognostic model. (B) Rank of risk score and distribution of groups. (C) Survival status of patients in different groups.
By using GSE76427, we established a validation model that classified patients into high- and low-risk score groups based on OS using the median risk score as the cut-off (Figure S1).
Functional Enrichment Analysis of DEIRGs
We further investigated the biological functions of 116 IRGs via GO annotation and KEGG pathway analyses. The obtained results showed that the extracellular region was the most frequent GO biological process category (P < 0.05). The top frequent, enriched GO networks are shown in Figures 4A, B.
Figure 4 Functional enrichment analysis of DEIRGs: (A, B) GO biological process categories.
KEGG pathway analyses led to the identification of top significantly enriched pathways (Figure 5A); the “cytokine-cytokine receptor interaction” pathway was the most enriched one. Based on the relationship between IRGs and KEGG pathways, we constructed a network using Cytoscape (Figure 5B).
Figure 5 KEGG pathway analysis of DEIGRs: (A) Relationship between IRGs and KEGG pathways. (B) The significantly enriched KEGG pathway.
DElncRNAs and DEmiRNAs in HCC
Significant DEGs were identified in 374 HCC and 50 adjacent non-cancer liver tissues from TCGA database by using the edgeR package in the R statistical environment. A total of 172 DElncRNAs (126 upregulated and 46 downregulated) and 251 DEmiRNAs (187 upregulated and 64 downregulated) were identified with the thresholds of |log2 FC| > 2.0 and adjusted P < 0.05. The distribution of DElncRNAs and DEmiRNAs was based on two parameters, FDR and logFC, and is represented as a volcano plot (Figures 6A, C). A heatmap was constructed to depict the expression data of DElncRNAs and DEmiRNAs (Figures 6B, D).
Figure 6 DElncRNAs and DEmiRNAs in patients with HCC: Heatmap (A) and volcano plot (B) of DElncRNAs in HCC and adjacent non-cancer liver tissues. Heatmap (C) and volcano plot (D) of DEmiRNAs in HCC and adjacent non-cancer liver tissues.
Construction of the ceRNA Network
To further explore how lncRNAs mediate mRNA by combining miRNA in HCC, an lncRNA–miRNA–mRNA (ceRNA) network was constructed based on the above data. We used 172 DElncRNAs that were retrieved from the miRcode database and used the Perl program to identify 35 pairs of interacting lncRNAs and miRNAs. Moreover, miRTarBase, miRDB, and TargetScan databases were used to identify target mRNAs based on 251 miRNAs. The final target genes were selected, which were then included in all three datasets. Finally, six DEIRGmRNAs were included in the ceRNA network (Figure 7C). In total, 11 lncRNA nodes, six miRNA nodes, and six mRNA nodes were identified as having differential expression profiles in the ceRNA networks (Figures 7A, B).
Figure 7 The ceRNA network of mRNA, miRNA, and lncRNA: (A) Relationship between survival-associated IRG mRNA, miRNA, and lncRNA. (B) IL-11/miRNA-204/LncRNA axis. (C) Venn diagram of DEmiRNA target mRNA and DEIRGmRNAs.
Construction of a PPI Network and Module Analyses
Based on the obtained DEIRGmRNAs, we used STRING to construct a PPI network; visualization was performed using Cytoscape. The Molecular Complex Detection plug-in for Cytoscape was then used to screen the modules of hub genes from the PPI network (Figures 8A, B).
Figure 8 The PPI network of DEIRGs. (A) DEIRGs and (B) hub genes.
Clinical Outcome of Patients With HCC Using the IRGPM
By using KM curve, the prognostic of TCGA model (Figure 9A) and GSE76427 model (Figure S2) were analyzed.
Figure 9 Survival analysis by the prognostic model: (A) Kaplan–Meier analysis using high- and low-risk score groups of patients with HCC. (B) ROC curve analysis.
In TCGA model, the area under the ROC curve was 0.818, suggesting that the TCGA IRGPM could be used for monitoring survival (Figure 9B). Univariate analyses (Table 3) showed that high American Joint Committee on Cancer (AJCC) stage [hazard ratio (HR) = 1.866; 95% confidence interval (CI) = 1.488–2.341; P < 0.001), high T stage (HR = 1.826; 95% CI = 1.472–2.265; P < 0.001), and high risk score (HR = 2.238; 95% CI = 1.293-3.872; P < 0.001) were significant risk factors for poor prognoses. In the multivariate Cox regression analysis (Table 2), high risk score (HR = 2.071; 95% CI = 1.183–3.625; P < 0.001) was found to be independently associated with worse OS. The risk scores were significantly higher for patients with advanced T stage (Figure 10A), AJCC stage (Figure 10B), and grade (Figure 10C).
Table 3 Univariate and Multivariate Cox regression analysis of survival associated IRGS.
Figure 10 Relationships between risk score and clinical features. (A) Relationships between risk score and T stage, (B) between risk score and AJCC stage, and (C) between risk score and grade.
Correlation Analysis of the IRGPM and Immune Cell Infiltration
The risk factors identified via the model were positively correlated with macrophages (r = 0.442, P < 0.001) (Figure 11A), neutrophils (r = 0.300, P < 0.001) (Figure 11B), DCs (r = 0.286, P < 0.001) (Figure 11C), B cells (r = 0.124, P = 0.024) (Figure 11D), CD4 T cells (r = 0.146, P = 0.008) (Figure 11E), and CD8 T cells (r = 0.225, P < 0.001) (Figure 11F).
Figure 11 Relationships between the prognostic model and immune cell infiltration: (A) macrophages, (B) neutrophils, (C) dendritic cells, (D) B cells, (E) CD4 T cells, and (F) CD8 T cells.
Analysis and Validation of Gene Expression
To further validate the expression of relevant key genes, mi-RNA and lncRNA in the prognostic model and ceRNA network, we randomly selected five mRNAs, one miRNA, and one lncRNA to measure the expression level use HCC tissue and adjacent normal tissues. We found that the expression is consistent with the TCGA database (Figures 12A–G).
Figure 12 The expression of IL11, IL17D, SPP1, MAPT, BIRC5, miR-204, and SNHG12. (A) the expression of IL11, (B) the expression of IL17D, (C) the expression of SPP1, (D) the expression of MAPT, (E) the expression of BIRC5, (F) the expression of miR-204, and (G) the expression of SNHG12.
Discussion
Although the relationship between IRGs and tumor prognosis has been explored to a certain extent, that between IGRs and HCC still remains unknown. In the present study, we attempted to establish a prognostic model using IRGs in patients with HCC. Previous studies have indicated that IRGs are promising as prognostic indicators of early stage lung cancer. In a study conducted by (Li B. et al., 2017; Li T. et al., 2017) a prognostic immune signature of 25 gene pairs comprising 40 unique genes was created using the meta-training data set; the immune signature stratified patients with early-stage non-squamous non-small cell lung cancer into high- and low-risk groups in terms of OS, and remained as an independent prognostic factor in multivariate analyses (Li B. et al., 2017).
Through KEGG pathway analysis, we found that IRGs were mainly enriched in the “cytokine-cytokine receptor interaction” pathway. Cytokines secreted by various tumor cells have synergistic effects with immune cells in the tumor microenvironment (TME), resulting in increased tumor activity (de Looff et al., 2019). Macrophages and neutrophils can regulate the TRAIL/TRAIL-R system through cytokines to eliminate tumor cells (de Looff et al., 2019). Li et al. points out that the network between cytokines and tumor immunity has a bistable. This model has shown that tumors can use the bistable state to generate immunosuppression. If this eliminates this interaction, the immune system can return to immune boosting (Li and Levine, 2017). In addition, the authors have established various models to illustrate the interaction between tumor cells and macrophages in the tumor microenvironment, through which the stable state of tumor cells can be eliminated (Li et al., 2019). Two genes, i.e., IL-11 and IL-17, in the prognostic model were noted to be involved in this signaling pathway. Previous studies have indicated that cytokine signaling pathways are closely related to inflammation-mediated HCC (Hu et al., 2017). Inflammation plays a key role in liver cancer development (Capece et al., 2013). Chronic inflammation caused by HBV or HCV infection is a major factor in the development of HCC (Marrero et al., 2018). IL-17 is secreted by TH17 cells and is closely related to tumor growth (Miossec et al., 2009). Hu et al. indicated that IL-17 promotes the proliferation of liver cancer cells (HBV+) by activating the IL-6/STAT3 signaling pathway (Hu et al., 2017), and Tian et al. found that IL-17 expression and promoter methylation were closely related to the progression of chronic HBV infection, particularly in patients with HCC (Tian et al., 2019). Moreover, Wang et al. revealed that the IL-11 has a chief role in postsurgical HCC recurrence; they also found that the IL-11/STAT3 signaling pathway promotes the proliferation of tumor cells and that the inhibition of IL-11/STAT3 signaling could reduce the postsurgical recurrence of HCC (Wang et al., 2019). IL-11 regulates a variety of immune cells to participate in the immune response. For example, IL-11 can regulate macrophage function by inhibiting IL-1β (Schwertschlag et al., 1999). In addition, IL-11 can directly affect CD4T cells, promote the production of TH1 cells, and inhibit TH2 cells (Bozza et al., 2001; Curti et al., 2001).
To further identify the key genes participating in the prognostic model, we established PPI and ceRNA networks. In the PPI network, we screened five hub genes: EGF, SPP1, BMP4, MMP, and CXCL12. Three of these genes are associated with prognosis, and SPP1, also known as OPN, was one of the key genes for establishing the prognostic model. High expression levels of SPP1 can be reportedly detected in many tumors, including HCC (Gotoh et al., 2002; Reinholz et al., 2002; Urquidi et al., 2002). SPP1 has been reported to mediate macrophage polarization and induce immune escape in lung adenocarcinoma (Zhang et al., 2017). By activating the PI3K/Akt signaling pathway, OPN can cause HCC metastasis (Yu et al., 2018).
In case of the ceRNA network, we used DEmiRNAs and DElncRNAs to construct an lncRNA–miRNA–mRNA network. Salmena et al. proposed that all types of RNAs compete with each other for miRNAs, resulting in large-scale trans-regulatory crosstalk across the transcriptome as a whole (Salmena et al., 2011). In the ceRNA network, we identified that the expression of IL-11 was upregulated whereas that of miRNA and SHNG12 was upregulated. Several studies have reported that miRNA-204 expression is downregulated in tumor tissues than in normal tissues, and it is thus generally regarded to be a tumor suppressor gene (Yin et al., 2014; Liu et al., 2015; Xia et al., 2015; Hong et al., 2019). Ge et al. indicated that miRNA-204 can attenuate the proliferation of HCC cells by inhibiting the lncRNA HOTTIP (Ge et al., 2015). Lan et al. found that the expression levels of SNHG12 is higher in HCC compared with adjacent normal tissues, SNHG12 can promote the proliferation and metastasis of HCC through target miR-199 (Lan et al., 2017). Therefore, we hypothesized that IL-11 and SHNG12 competitively bind to miRNA204, resulting in the downregulation of miRNA204 expression with impaired anticancer effects, via further promoting the occurrence and development of HCC.
In the past decade, liver cancer treatment has shown marginal improvement, particularly for patients who are inoperable or at stage IV. Sorafenib has been reported to prolong the OS of patients with advanced HCC in the Asia-Pacific region (Cheng et al., 2009), and a study reported the median survival time for lenvatinib was 13.6 months, being non-inferior to sorafenib (12.3 months), in patients with advanced HCC (Kudo et al., 2018). In the past decade, the emergence of immunotherapy has revolutionized cancer treatment, and immunological checkpoint inhibitors appear promising for treating patients with liver cancer (El-Khoueiry et al., 2017; Zhu et al., 2018).
A recent study proposed that the functional role of innate lymphoid cells in antitumor immunity is complex (Gajewski et al., 2013). miRNAs in the tumor microenvironment have been reported to have a crucial role in the development of cancer and its progression (Kuninty et al., 2016). Interestingly, in this study, we found that our prognostic model was positively correlated with the infiltration of six subtypes of tumor-infiltrating immune cells.
Tumor-associated macrophages are the major components of tumor inflammatory infiltration and are key mediators between tumors and inflammation (Mantovani et al., 1992). Guo et al. found that M2 macrophages in the HCC microenvironment can secrete large amounts of IL-17 and inhibit oxaliplatin-induced apoptosis by activating chaperone-mediated autophagy (Guo B. et al., 2017). Tumor-associated neutrophils are also an important part of the tumor microenvironment and can promote various biological activities of tumors (Galdiero et al., 2013). Neutrophils can secrete proinflammatory and immunoregulatory factors, such as neutrophil elastase, hepatocyte growth factor, and β2-integrins, which have paracrine effects on tumor biology (Wislez et al., 2003; Houghton et al., 2010; Strell et al., 2010). Tumor-associated neutrophils/macrophages can reportedly promote the survival and proliferation of malignant tumors; promote the remodeling of extracellular matrix, blood vessels, lymphangiogenesis, and invasion; and also facilitate the transfer of adaptive immunity (Mantovani et al., 2008; Qian and Pollard, 2010; Mantovani et al., 2011).
Zhou et al. indicated that intratumoral infiltration of plasmacytoid DCs is a novel indicator of poor prognosis in patients with HCC. They also reported that an increase in intratumoral plasmacytoid DCs was associated with increased intratumoral infiltration of Foxp3+ regulatory T (Treg) cells and IL-17-producing cells (Zhou et al., 2019). Tumor-infiltrating B lymphocytes evidently play a dual role in some types of cancers: they combat cancer growth via antigen-restricted tumoricidal immune response or promote tumor progression via the induction of immune suppression (Guo S. et al., 2017). Another study reported that patients with HCC exhibited markedly higher levels of IL-10-expressing B cells as compared to healthy controls (Xue et al., 2016).
T cells have many subpopulations, including cytotoxic T cells (CTL), helper T (TH) cells, and Treg cells. Numerous previous studies have confirmed that CD4+CTLs or CD8+CTLs can inhibit the proliferation of HCC cells (Fu et al., 2013; Yao et al., 2017). Treg cells have a distinct ability to evade immune suppression to impair antitumor ability, helping tumor cells to escape immune surveillance (Du et al., 2012). Previous studies have indicated that high levels of CD4+CD25+Treg cells are positively correlated with the expression of HBV DNA in the serum of patients with chronic HBV infection, and antiviral therapy can reduce the frequency of Treg cells (Yan et al., 2014). Zhou et al. demonstrated that tumor-associated neutrophils can recruit Treg cells and macrophages, leading to progression of HCC and promoting resistance to sorafenib (Zhou et al., 2016). The complexity of the tumor microenvironment and immune system makes the response to treatment difficult to predict; thus, personalized treatment of tumors is imminent.
This study has some limitations that should be borne in mind in interpreting the findings. Firstly, the model established by us showed a positive correlation with the infiltration of T cells and DCs. We know that T cells and DCs have different subtypes. However, a progressive deficit in CD4+ CTLs induced by increased FoxP3+ Treg cells is reportedly correlated with poor survival and high recurrence rates in patients with HCC (Fu et al., 2013). We propose that subpopulations of cells should be individually studied in future studies. Finally, we identified some key genes and also explained their potential mechanisms; however, further studies on these are warranted.
In summary, we used survival-associated IRGs to establish a prognostic model. We identified some key genes via the PPI and ceRNA networks and that the expression of IL-11 was upregulated whereas that of miRNA and LINC00261 was downregulated; moreover, the model was positively correlated with the infiltration of several subtypes of tumor-infiltrating immune cells. The prognostic model was also related to clinical characteristics such as AJCC stage, T stage, and HCC differentiation. We believe that the model provides sufficient information for prognosis prediction and immunotherapy of patients with HCC.
Data Availability Statement
The datasets analysed for this study were obtained from TCGA data portal (https://portal.gdc.cancer.gov/) and GEO data portal (https://www.ncbi.nlm.nih.gov/geo/).
Ethics Statement
The studies involving human participants were reviewed and approved by Medical Ethics Committees of the Affiliated Suzhou Hospital of Nanjing Medical University. The patients/participants provided their written informed consent to participate in this study.
Author Contributions
Z-LD and Y-QH contributed toward the conception and design of this study. JZ and J-JW organized the database. WS performed statistical analyses. HW, W-JW, and T-YH wrote the first draft of the manuscript. T-YH and WS wrote sections of the manuscript. All authors contributed to manuscript revision and have read and approved the submitted version.
Funding
This study was supported by the Science and Education for Health Foundation of Suzhou for Youth (grant no. KJXW2018030 and KJXW2018032), the Science and Technology Project Foundation of Suzhou (grant no. SS201651), Suzhou Oncology Clinical Center (grant no. Szzx201506), Project of Clinical Expert Team Introduction in Suzhou (grant no. SZYJTD201807), Project of science and technology development plan in Suzhou (grant no. SYSD2018138) and the Jiangsu Province Medical key discipline (grant no. ZDXKC2016007).
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.
Acknowledgements
We thank Wei Shan for his assistance with statistical analyses.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00055/full#supplementary-material
Supplementary Figure S1 | Construction of a validation model of GSE76427. Heatmap (A) of survival-associated IRGs in the prognostic model. (B) Rank of risk score and distribution of groups. (C) Survival status of patients in different groups.
Supplementary Figure S2 | The prognosis of validation model.
References
Bhattacharya, S., Andorf, S., Gomes, L., Dunn, P., Schaefer, H., Pontius, J., et al. (2014). ImmPort: disseminating data to the public for the future of immunology. Immunol. Res. 58 (2–3), 234–239. doi: 10.1007/s12026-014-8516-1
Bozza, M., Bliss, J. L., Dorner, A. J., Trepicchio, W. L. (2001). Interleukin-11 modulates Th1/Th2 cytokine production from activated CD4+ T cells. J. Interferon Cytokine Res. 21 (1), 21–30. doi: 10.1089/107999001459123
Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 68 (6), 394–424. doi: 10.3322/caac.21492
Capece, D., Fischietti, M., Verzella, D., Gaggiano, A., Cicciarelli, G., Tessitore, A., et al. (2013). The inflammatory microenvironment in hepatocellular carcinoma: a pivotal role for tumor-associated macrophages. BioMed. Res. Int. 2013, 187204. doi: 10.1155/2013/187204
Chen, W., Zheng, R., Baade, P. D., Zhang, S., Zeng, H., Bray, F., et al. (2016). Cancer statistics in China, 2015. CA Cancer J. Clin. 66 (2), 115–132. doi: 10.3322/caac.21338
Cheng, A. L., Kang, Y. K., Chen, Z., Tsao, C. J., Qin, S., Kim, J. S., et al. (2009). Efficacy and safety of sorafenib in patients in the Asia-Pacific region with advanced hepatocellular carcinoma: a phase III randomised, double-blind, placebo-controlled trial. Lancet Oncol. 10 (1), 25–34. doi: 10.1016/S1470-2045(08)70285-7
Chi, H. C., Chen, S. L., Cheng, Y. H., Lin, T. K., Tsai, C. Y., Tsai, M. M., et al. (2016). Chemotherapy resistance and metastasis-promoting effects of thyroid hormone in hepatocarcinoma cells are mediated by suppression of FoxO1 and Bim pathway. Cell Death Dis. 7 (8), e2324. doi: 10.1038/cddis.2016.227
Curti, A., Ratta, M., Corinti, S., Girolomoni, G., Ricci, F., Tazzari, P., et al. (2001). Interleukin-11 induces Th2 polarization of human CD4(+) T cells. Blood 97 (9), 2758–2763. doi: 10.1182/blood.V97.9.2758
de Looff, M., de Jong, S., Kruyt, F. A. E. (2019). Multiple interactions between cancer cells and the tumor microenvironment modulate TRAIL signaling: implications for TRAIL receptor targeted therapy. Front. Immunol. 10, 1530. doi: 10.3389/fimmu.2019.01530
Dennis, G., Jr., Sherman, B. T., Hosack, D. A., Yang, J., Gao, W., Lane, H. C., et al. (2003). DAVID: Database for Annotation, Visualization, and Integrated Discovery. Genome Biol. 4 (5), P3. doi: 10.1186/gb-2003-4-5-p3
Du, Y., Chen, X., Huang, Z. M., Ye, X. H., Niu, Q. (2012). Increased frequency of Foxp3+ regulatory T cells in mice with hepatocellular carcinoma. Asian Pac. J. Cancer Prev. 13 (8), 3815–3819. doi: 10.7314/APJCP.2012.13.8.3815
El-Khoueiry, A. B., Sangro, B., Yau, T., Crocenzi, T. S., Kudo, M., Hsu, C., et al. (2017). Nivolumab in patients with advanced hepatocellular carcinoma (CheckMate 040): an open-label, non-comparative, phase 1/2 dose escalation and expansion trial. Lancet 389 (10088), 2492–2502. doi: 10.1016/S0140-6736(17)31046-2
Fu, J., Zhang, Z., Zhou, L., Qi, Z., Xing, S., Lv, J., et al. (2013). Impairment of CD4+ cytotoxic T cells predicts poor survival and high recurrence rates in patients with hepatocellular carcinoma. Hepatology 58 (1), 139–149. doi: 10.1002/hep.26054
Gajewski, T. F., Schreiber, H., Fu, Y. X. (2013). Innate and adaptive immune cells in the tumor microenvironment. Nat. Immunol. 14 (10), 1014–1022. doi: 10.1038/ni.2703
Galdiero, M. R., Garlanda, C., Jaillon, S., Marone, G., Mantovani, A. (2013). Tumor associated macrophages and neutrophils in tumor progression. J. Cell Physiol. 228 (7), 1404–1412. doi: 10.1002/jcp.24260
Ge, Y., Yan, X., Jin, Y., Yang, X., Yu, X., Zhou, L., et al. (2015). MiRNA-192 [corrected] and miRNA-204 directly suppress lncRNA HOTTIP and interrupt GLS1-mediated glutaminolysis in hepatocellular carcinoma. PloS Genet. 11 (12), e1005726. doi: 10.1371/journal.pgen.1005726
Gotoh, M., Sakamoto, M., Kanetaka, K., Chuuma, M., Hirohashi, S. (2002). Overexpression of osteopontin in hepatocellular carcinoma. Pathol. Int. 52 (1), 19–24. doi: 10.1046/j.1440-1827.2002.01316.x
Guo, B., Li, L., Guo, J., Liu, A., Wu, J., Wang, H., et al. (2017). M2 tumor-associated macrophages produce interleukin-17 to suppress oxaliplatin-induced apoptosis in hepatocellular carcinoma. Oncotarget 8 (27), 44465–44476. doi: 10.18632/oncotarget.17973
Guo, S., Contratto, M., Miller, G., Leichman, L., Wu, J. (2017). Immunotherapy in pancreatic cancer: unleash its potential through novel combinations. World J. Clin. Oncol. 8 (3), 230–240. doi: 10.5306/wjco.v8.i3.230
Heagerty, P. J., Lumley, T., Pepe, M. S. (2000). Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics 56 (2), 337–344. doi: 10.1111/j.0006-341X.2000.00337.x
Hong, B. S., Ryu, H. S., Kim, N., Kim, J., Lee, E., Moon, H., et al. (2019). Tumor suppressor miRNA-204-5p regulates growth, metastasis, and immune microenvironment remodeling in breast cancer. Cancer Res. 79 (7), 1520–1534. doi: 10.1158/0008-5472.CAN-18-0891
Houghton, A. M., Rzymkiewicz, D. M., Ji, H., Gregory, A. D., Egea, E. E., Metz, H. E., et al. (2010). Neutrophil elastase-mediated degradation of IRS-1 accelerates lung tumor growth. Nat. Med. 16 (2), 219–223. doi: 10.1038/nm.2084
Hu, Z., Luo, D., Wang, D., Ma, L., Zhao, Y., Li, L. (2017). IL-17 Activates the IL-6/STAT3 signal pathway in the proliferation of hepatitis B virus-related hepatocellular carcinoma. Cell Physiol. Biochem. 43 (6), 2379–2390. doi: 10.1159/000484390
Kudo, M., Finn, R. S., Qin, S., Han, K. H., Ikeda, K., Piscaglia, F., et al. (2018). Lenvatinib versus sorafenib in first-line treatment of patients with unresectable hepatocellular carcinoma: a randomised phase 3 non-inferiority trial. Lancet 391 (10126), 1163–1173. doi: 10.1016/S0140-6736(18)30207-1
Kuninty, P. R., Schnittert, J., Storm, G., Prakash, J. (2016). MicroRNA Targeting to Modulate Tumor Microenvironment. Front. Oncol. 6, 3. doi: 10.3389/fonc.2016.00003
Lan, T., Ma, W., Hong, Z., Wu, L., Chen, X., Yuan, Y. (2017). Long non-coding RNA small nucleolar RNA host gene 12 (SNHG12) promotes tumorigenesis and metastasis by targeting miR-199a/b-5p in hepatocellular carcinoma. J. Exp. Clin. Cancer Res. 36 (1), 11. doi: 10.1186/s13046-016-0486-9
Li, B., Cui, Y., Diehn, M., Li, R. (2017). Development and validation of an individualized immune prognostic signature in early-stage nonsquamous non-small cell lung cancer. JAMA Oncol. 3 (11), 1529–1537. doi: 10.1001/jamaoncol.2017.1609
Li, T., Fan, J., Wang, B., Traugh, N., Chen, Q., Liu, J. S., et al. (2017). TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 77 (21), e108–e110. doi: 10.1158/0008-5472.CAN-17-0307
Li, X., Jolly, M. K., George, J. T., Pienta, K. J., Levine, H. (2019). Computational modeling of the crosstalk between macrophage polarization and tumor cell plasticity in the tumor microenvironment. Front. Oncol. 9, 10. doi: 10.3389/fonc.2019.00010
Li, X., Levine, H. (2017). Bistability of the cytokine-immune cell network in a cancer microenvironment. IPO Sci. 3, 024002. doi: 10.1088/2057-1739/aa6c07
Lin, P., Guo, Y. N., Shi, L., Li, X. J., Yang, H., He, Y., et al. (2019). Development of a prognostic index based on an immunogenomic landscape analysis of papillary thyroid cancer. Aging (Albany N.Y.) 11 (2), 480–500. doi: 10.18632/aging.101754
Liu, L., Wang, J., Li, X., Ma, J., Shi, C., Zhu, H., et al. (2015). MiR-204-5p suppresses cell proliferation by inhibiting IGFBP5 in papillary thyroid carcinoma. Biochem. Biophys. Res. Commun. 457 (4), 621–626. doi: 10.1016/j.bbrc.2015.01.037
Mantovani, A., Bottazzi, B., Colotta, F., Sozzani, S., Ruco, L. (1992). The origin and function of tumor-associated macrophages. Immunol. Today 13 (7), 265–270. doi: 10.1016/0167-5699(92)90008-U
Mantovani, A., Allavena, P., Sica, A., Balkwill, F. (2008). Cancer-related inflammation. Nature 454 (7203), 436–444. doi: 10.1038/nature07205
Mantovani, A., Cassatella, M. A., Costantini, C., Jaillon, S. (2011). Neutrophils in the activation and regulation of innate and adaptive immunity. Nat. Rev. Immunol. 11 (8), 519–531. doi: 10.1038/nri3024
Marrero, J. A., Kulik, L. M., Sirlin, C. B., Zhu, A. X., Finn, R. S., Abecassis, M. M., et al. (2018). Diagnosis, staging, and management of hepatocellular carcinoma: 2018 practice guidance by the American Association for the study of liver diseases. Hepatology 68 (2), 723–750. doi: 10.1002/hep.29913
Miossec, P., Korn, T., Kuchroo, V. K. (2009). Interleukin-17 and type 17 helper T cells. N. Engl. J. Med. 361 (9), 888–898. doi: 10.1056/NEJMra0707449
Nagao, M., Nakajima, Y., Kanehiro, H., Hisanaga, M., Aomatsu, Y., Ko, S., et al. (2000). The impact of interferon gamma receptor expression on the mechanism of escape from host immune surveillance in hepatocellular carcinoma. Hepatology 32 (3), 491–500. doi: 10.1053/jhep.2000.16470
Qian, B. Z., Pollard, J. W. (2010). Macrophage diversity enhances tumor progression and metastasis. Cell 141 (1), 39–51. doi: 10.1016/j.cell.2010.03.014
Reinholz, M. M., Iturria, S. J., Ingle, J. N., Roche, P. C. (2002). Differential gene expression of TGF-beta family members and osteopontin in breast tumor tissue: analysis by real-time quantitative PCR. Breast Cancer Res. Treat 74 (3), 255–269. doi: 10.1023/A:1016339120506
Robinson, M. D., McCarthy, D. J., Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26 (1), 139–140. doi: 10.1093/bioinformatics/btp616
Salmena, L., Poliseno, L., Tay, Y., Kats, L., Pandolfi, P. P. (2011). A ceRNA hypothesis: the Rosetta stone of a hidden RNA language? Cell 146 (3), 353–358. doi: 10.1016/j.cell.2011.07.014
Schwertschlag, U. S., Trepicchio, W. L., Dykstra, K. H., Keith, J. C., Turner, K. J., Dorner, A. J. (1999). Hematopoietic, immunomodulatory and epithelial effects of interleukin-11. Leukemia 13 (9), 1307–1315. doi: 10.1038/sj.leu.2401514
Strell, C., Lang, K., Niggemann, B., Zaenker, K. S., Entschladen, F. (2010). Neutrophil granulocytes promote the migratory activity of MDA-MB-468 human breast carcinoma cells via ICAM-1. Exp. Cell Res. 316 (1), 138–148. doi: 10.1016/j.yexcr.2009.09.003
Tian, C. H., Dai, J., Zhang, W., Liu, Y., Yang, Y. (2019). Expression of IL-17 and its gene promoter methylation status are associated with the progression of chronic hepatitis B virus infection. Med. (Baltimore) 98 (23), e15924. doi: 10.1097/MD.0000000000015924
Urquidi, V., Sloan, D., Kawai, K., Agarwal, D., Woodman, A. C., Tarin, D., et al. (2002). Contrasting expression of thrombospondin-1 and osteopontin correlates with absence or presence of metastatic phenotype in an isogenic model of spontaneous human breast cancer metastasis. Clin. Cancer Res. 8 (1), 61–74.
Wang, B., Zhao, X. P., Fan, Y. C., Zhang, J. J., Zhao, J., Wang, K. (2013). IL-17A but not IL-22 suppresses the replication of hepatitis B virus mediated by over-expression of MxA and OAS mRNA in the HepG2.2.15 cell line. Antiviral Res. 97 (3), 285–292. doi: 10.1016/j.antiviral.2012.12.018
Wang, D., Zheng, X., Fu, B., Nian, Z., Qian, Y., Sun, R., et al. (2019). Hepatectomy promotes recurrence of liver cancer by enhancing IL-11-STAT3 signaling. EBioMedicine 46, 119–132. doi: 10.1016/j.ebiom.2019.07.058
Wislez, M., Rabbe, N., Marchal, J., Milleron, B., Crestani, B., Mayaud, C., et al. (2003). Hepatocyte growth factor production by neutrophils infiltrating bronchioloalveolar subtype pulmonary adenocarcinoma: role in tumor progression and death. Cancer Res. 63 (6), 1405–1412.
Xia, Z., Liu, F., Zhang, J., Liu, L. (2015). Decreased expression of MiRNA-204-5p contributes to glioma progression and promotes glioma cell growth, migration and invasion. PloS One 10 (7), e0132399. doi: 10.1371/journal.pone.0132399
Xue, H., Lin, F., Tan, H., Zhu, Z. Q., Zhang, Z. Y., Zhao, L. (2016). Overrepresentation of IL-10-expressing B cells suppresses cytotoxic CD4+ T cell activity in HBV-Induced hepatocellular carcinoma. PloS One 11 (5), e0154815. doi: 10.1371/journal.pone.0154815
Yan, Z., Zhou, J., Zhang, M., Fu, X., Wu, Y., Wang, Y. (2014). Telbivudine decreases proportion of peripheral blood CD4+CD25+CD127low T cells in parallel with inhibiting hepatitis B virus DNA. Mol. Med. Rep. 9 (5), 2024–2030. doi: 10.3892/mmr.2014.2042
Yao, W., He, J. C., Yang, Y., Wang, J. M., Qian, Y. W., Yang, T., et al. (2017). The prognostic value of tumor-infiltrating lymphocytes in hepatocellular carcinoma: a systematic review and meta-analysis. Sci. Rep. 7 (1), 7525. doi: 10.1038/s41598-017-08128-1
Yin, Y., Zhang, B., Wang, W., Fei, B., Quan, C., Zhang, J., et al. (2014). miR-204-5p inhibits proliferation and invasion and enhances chemotherapeutic sensitivity of colorectal cancer cells by downregulating RAB22A. Clin. Cancer Res. 20 (23), 6187–6199. doi: 10.1158/1078-0432.CCR-14-1030
Yu, X., Zheng, Y., Zhu, X., Gao, X., Wang, C., Sheng, Y., et al. (2018). Osteopontin promotes hepatocellular carcinoma progression via the PI3K/AKT/Twist signaling pathway. Oncol. Lett. 16 (4), 5299–5308. doi: 10.3892/ol.2018.9281
Zhang, Y., Du, W., Chen, Z., Xiang, C. (2017). Upregulation of PD-L1 by SPP1 mediates macrophage polarization and facilitates immune escape in lung adenocarcinoma. Exp. Cell Res. 359 (2), 449–457. doi: 10.1016/j.yexcr.2017.08.028
Zhou, S. L., Zhou, Z. J., Hu, Z. Q., Huang, X. W., Wang, Z., Chen, E. B., et al. (2016). Tumor-associated neutrophils recruit macrophages and T-regulatory cells to promote progression of hepatocellular carcinoma and resistance to sorafenib. Gastroenterology 1501646–1501658 (7), e1617. doi: 10.1053/j.gastro.2016.02.040
Zhou, Z. J., Xin, H. Y., Li, J., Hu, Z. Q., Luo, C. B., Zhou, S. L. (2019). Intratumoral plasmacytoid dendritic cells as a poor prognostic factor for hepatocellular carcinoma following curative resection. Cancer Immunol. Immunother. 68 (8), 1223–1233. doi: 10.1007/s00262-019-02355-3
Zhu, A. X., Finn, R. S., Edeline, J., Cattan, S., Ogasawara, S., Palmer, D., et al. (2018). Pembrolizumab in patients with advanced hepatocellular carcinoma previously treated with sorafenib (KEYNOTE-224): a non-randomised, open-label phase 2 trial. Lancet Oncol. 19 (7), 940–952. doi: 10.1016/S1470-2045(18)30351-6
Keywords: hepatocellular carcinoma, immune-related genes, ceRNA network, prognostic model, The Cancer Genome Atlas (TCGA)
Citation: Wang W-J, Wang H, Hua T-Y, Song W, Zhu J, Wang J-J, Huang Y-Q and Ding Z-L (2020) Establishment of a Prognostic Model Using Immune-Related Genes in Patients With Hepatocellular Carcinoma. Front. Genet. 11:55. doi: 10.3389/fgene.2020.00055
Received: 13 November 2019; Accepted: 17 January 2020;
Published: 25 February 2020.
Edited by:
Yun Xiao, Harbin Medical University, ChinaReviewed by:
Mohit Kumar Jolly, Indian Institute of Science (IISc), IndiaShaogui Wan, Gannan Medical University, China
Copyright © 2020 Wang, Wang, Hua, Song, Zhu, Wang, Huang and Ding. 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: Yue-Qing Huang, aHVhbmd5cV9zekAxNjMuY29t; Zhi-Liang Ding, c3pfbmpteTE5ODVAMTYzLmNvbQ==
†These authors have contributed equally to this work