Bioinformatic Analysis of Hepatocellular Carcinoma Cell Lines to the Efficacy of Nimotuzumab

Introduction

Primary liver cancers, of which hepatocellular carcinoma (HCC) is the main type, are the third leading cause of cancer-related mortality worldwide.^1,2 It is estimated that new cases of HCC reached 906,000, with mortality as high as 830,000 annually.^2,3 The incidence of HCC has increased over the last years, and over half of the global incidence and mortality of HCC occurs in China, attributed to HBV infection.⁴ Radical therapeutic options for HCC patients are radical surgical resection, local ablation or liver transplantation, but these are confined to a small proportion of patients with early-stage HCC. HCC patients at an advanced stage will receive the systemic treatment, including chemotherapy, targeted therapy and immunotherapy, as a means of prolonging life, although most current drugs show limited efficacy.^5,6 Due to the high morbidity and high mortality of HCC, novel therapy is an urgent need in HCC treatment.

The epidermal growth factor receptor (EGFR) is a well-known classic onco-therapeutic target in various tumors. It is a cell-surface receptor belonging to ErbB family. It plays a vital role as a regulator of pathways involved in cancer pathogenesis and progression. Hence, specific EGFR inhibition is one of the key targets for cancer therapy. Anti-EGFR monoclonal antibodies (mAbs) and tyrosine kinase inhibitors (TKIs) are two major approaches that have demonstrated benefits for targeting EGFR. mAbs, including cetuximab, panitumumab, necitumumab and nimotuzumab, which ligand-bind directly to the extracellular domain of EGFR, block the activation of downstream signaling pathways.⁷ TKIs, such as gefitinib, afatinib, osimertinib, could internalize into cells and inhibit the phosphorylated activation by binding with the intracellular tyrosine kinase site of EGFR.⁸ However, these anti-EGFR agents have not been proved in a positive evaluation in HCC despite the expression of EGFR being commonly high.^9–12 Some studies showed that the complex signaling pathways are the primary cause of drug resistance in HCC cells.^13,14

Nimotuzumab is a humanized IgG1 isotype anti-EGFR mAb with two ligands, while other anti-EGFR mAbs only have one, developed at the Center of Molecular Immunology in Cuba.¹⁵ Nimotuzumab was approved for the treatment of glioma, esophageal cancer, squamous cell carcinoma of the head and neck and pancreatic cancer in Cuba and some of these conditions in China.^16–18 Nimotuzumab’s data on HCC-related clinical trials were not available. However, it was observed that nimotuzumab was highly effective in a patient with advanced HCC who progressed after multi-line and multi-course treatment in our hospital. Another patient in the same situation was reported as well.¹⁹ It was of interest to us, and we aimed to explore the efficacy and resistance mechanism of nimotuzumab in HCC.

We selected six HCC cell lines for nimotuzumab susceptibility experiments and tested the EGFR protein on cell surface. Then we conducted RNA sequencing of all HCC cell lines. After that, bioinformatic analysis was performed to compare the differences in gene expression between the sensitive and the resistant cell lines.

Materials and Methods

Cell Lines and Cell Culture

The human HCC cell lines Hep 3B2.1–7, Li-7, PLC/PRF/5, SK-HEP-1, SNU-182, and SNU-387 were obtained from the National Collection of Authenticated Cell Cultures in Shanghai. After resuscitation, all the cell lines were cultured in DMEM medium supplemented with 10% fetal bovine serum (Biological Industries) containing 1% double resistance in a 5% humidified CO₂ incubator at 37°C. All cell lines were authenticated by short tandem repeat analysis. These cells were observed under a microscope every day. Then they were passaged and cultured when grown to mid-log phase (OD600 nm 0.4–0.6). After third passage, cells were used for subsequent experiments.

Flow Cytometry Analysis

Flow cytometry analyzed the expressions of EGFR protein on cell surface of the cell lines. After separation, cultivation, and passage, the third generation of cell lines were digested by 0.05% pancreatin and the obtained cell suspension was centrifuged in 1000 rpm for 5 min. The cell density was adjusted to 1×10⁶/mL, after which 1mL cell was added in the labeled flow cytometry tubes. After centrifugation the supernatant was discarded. The cells were washed twice with PBS. Cells were stained with fluorochrome-conjugated monoclonal antibodies: EGFR and IgG2b (BD Pharmingen, USA) as an isotype control in accordance with the manufacturer’s specifications. These antibodies were proven by their manufacturers to be cross-reactive with mouse antigens. Samples were mixed well and incubated at 4°C for 10–15 min in the dark, and they were washed twice with PBS prior to analysis. Thousands of fluorescent cells were accumulated and detected by flow cytometer (Canto II, BD Biosciences, USA).

Cell Proliferation Assay (CCK-8)

The cells from the cell lines were seeded at a density of 5×10⁴ cells/mL in 100 μL/well (5000 cells/well) in 96-well plates. Then 200 μL of PBS were added into wells in the plate which were then placed in the incubator. Cell proliferation assays were classified into two groups: control and nimotuzumab. The negative controls (C0) contained medium that was not supplemented with a drug. Nimotuzumab were amended with four different concentrations (C1-C4): 25 μg/mL, 50 μg/mL, 100 μg/mL, and 200 μg/mL. The mediums were supplemented with different concentrations of nimotuzumab. After the cell attached to the wall, the solutions with different drugs were added into the well (three replicated per group) and continued to be cultured. Supernatants were discarded after 72 h. 90 μL of fresh growth medium and 10 μL CCK-8 solution were added to each well. Then plates were incubated in an incubator for 1–4 h and the color change in the reaction solution was noted by visual observation. The optical density (OD) of each well was measured at 450 and 600 nm. The experiments were independently repeated at least 3 times.

RNA-Sequencing (RNA-Seq)

Total RNA from cell lines were extracted and purified using an RNA extraction kit (Qiagen). Then mRNAs were enriched by polyA selection. Library preparations were conducted using an Illumina TruSeq RNA sample prep kit (Illumina, Beijing, China) following the manufacturer’s instructions and subjected to sequencing on a HiSeq 4000 platform (Novogene).

Bioinformatic Analysis

RNA-seq raw reads were processed using fastp to filter out sequencing adapters, low-quality reads.²⁰ Clean reads were aligned to the human reference genome (GRCh38.91) by HISAT2.²¹ The relative abundances of genes were measured in FPKM (Fragments per kilobase per million mapped reads) using StringTie.²²

Differentially expressed genes (DEGs) between groups were evaluated with the R (version 4.0.2) package “DESeq2”, “edgeR”, “limma” using FPKM data. Genes were nominated to be differentially expressed with two cut-off values (fold-change > 2 and p <0.02).²³ Gene ontology and KEGG enrichment analyses were performed using the “clusterProfiler” package in R.^24–26 Gene Set Enrichment Analysis (GSEA) was conducted using GSEA 4.1.0 software with gene sets annotated in the Molecular Signatures Database (MSigDB v.7.2).^27,28 The cut-off values for GSEA results were fold-change >1.8, PCT score >0.3 and false discovery rate (FDR) q-value <0.1. Heatmaps were generated in R package “pheatmap”.

The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO) Data Analysis

Raw counts of RNA-sequencing data and corresponding clinical information from TCGA-LIHC were obtained from The Cancer Genome Atlas (TCGA) dataset (https://portal.gdc.cancer.gov/) in May 2021, in which the method of acquisition and application complied with the guidelines and policies. All the above analysis methods and R package were implemented by R foundation for statistical computing and ggplot2.

All microarray data from GSE102079 was submitted to the GEO database (http://www.ncbi.nih.gov/geo). The raw data were downloaded as MINiML files. The extracted data were normalized and processed by log2 transformation. The microarray data were normalized using the “preprocessCore” package in R software. Probes were converted to gene symbols according to the platform annotation information of the normalized data.

Statistical Analysis

Data were analyzed and graphed using GraphPad Prism 9. Data were presented as mean ± SEM. Experiments that involved a single comparison between two groups were analyzed with t-tests, whereas differences between multiple groups were analyzed with one-way analysis of variance (ANOVA). A p-value <0.05 was regarded as statistically significant.

Results

Expressions of EGFR Protein Have No Correspondence with Efficacy of Nimotuzumab

To explore whether the expression of EGFR protein contacted with the anti-EGFR agency, nimotuzumab, we analyzed the expressions of EGFR protein on every cell face of the six cell lines by flow cytometry. Results from flow cytometry (Figure 1A) revealed that all six HCC cell lines expressed EGFR protein relative to isotype. There were no significant differences of expression level among these cell lines. Hep 3B2.1–7 showed the highest mean fluorescence intensity (MFI) and had a significant difference with other liver cell lines besides PLC/PRF/5. PLC/PRF/5 had a higher MFI than the rest of the four cell lines but these did not differ significantly from each other. In conclusion, Hep 3B2.1–7 had a higher expression of EGFR protein than other liver cell lines and EGFR protein of other liver cell lines did not show a significant difference.

Figure 1 Expressions of EGFR protein have no correspondence with efficacy of nimotuzumab. (A) Representative flow cytometric histograms showed the expressions of EGFR protein on cells surface of six HCC cell lines (Hep 3B2.1–7, Li-7, PLC/PRF/5, SK-HEP-1, SNU-182, SNU-387). Mean fluorescence intensity (MFI) of flow cytometric histograms (mean ± SEM of 3 independent experiments). Statistical differences were determined by one-way ANOVA and Student’s t-test, with *p<0.05, **p<0.01, ***p<0.0001 against isotype; #p<0.05, ##p<0.01, ###p<0.0001 against Hep 3B2.1–7. (B) Cell proliferation was evaluated using the CCK-8 assay. Statistical differences were determined by Student’s t-test, with *p<0.05, **p<0.01, ***p<0.0001 against Hep 3B2.1–7 C0; #p<0.05, ##p<0.01, ###p<0.0001 against PLC/PRF/5 C0 (C0 = negative control).

We then examined the impact of nimotuzumab on HCC cell proliferation by CCK-8 assay. The cells were divided into two sets, treated for 72 h with either control or nimotuzumab with different concentrations. Cell proliferation assay (Figure 1B) results demonstrated that nimotuzumab did not inhibit HCC cells proliferation at all in Li-7, SK-HEP-1, SNU-182 and SNU-387 cell lines at every concentration. Hep 3B2.1–7 and PLC/PRF/5 cell lines were inhibited by nimotuzumab at every concentration. With the concentration increase, these two cell lines proliferated less. Hep 3B2.1–7 cell line proliferation ratio was 75.92% and PLC/PRF/5 cell line was 71.43% on nimotuzumab 200 ug/mL.

According to the EGFR protein expression above, the sensitive cell line Hep 3B2.1–7 had a relatively higher expression on EGFR protein. However, PLC/PRF/5 did not show a significant difference with other resistant cell lines in EGFR protein. So, these results showed no correspondence between the expression of EGFR proteins and drug efficacy of nimotuzumab which was an anti-EGFR agent.

Differentially Expressed Genes Between Resistant and Sensitive HCC Cell Lines

To explore the differences between resistant HCC cell lines (Li-7, SK-HEP-1, SNU-182 and SNU-387) and sensitive cell lines (Hep 3B2.1–7 and PLC/PRF/5) on a genetic level, we analyzed RNA-sequences from the six HCC cell lines. After processing reads, we removed low-quality sequences, and the clean read mapped more than 95% of the raw data. Protein-coding genes and non-coding RNAs were all included.

We separated the six HCC cell lines into two groups, resistant and sensitive. The resistant group was chosen as reference. After that, we performed DESeq2, edgeR and limma voom in R to get differentially expressed genes (DEGs) between two groups. Results from three ways were intersected. A total of 480 DEGs (p<0.05, fold change ≥2) were detected with 114 up-regulated genes and 366 down-regulated genes in the sensitive group (Figure 2).

Figure 2 DEGs between sensitive and resistant HCC cell lines were analyzed by DESeq2, edgeR and limma voom packages in R: (A–C) Heatmap of DEGs and (D–F) Volcano plots of DEGs in three ways. Intersection of three ways of DEGs: (G) Heatmap and (H and I) Venn diagram. A total of 114 up-regulated genes (p<0.05, fold change ≥2) and 366 down-regulated genes (p<0.05, fold change ≥2) were detected in sensitive group (resistant group as reference). The color-ratio bar and plot indicate intensity of gene up-regulation (red) and down-regulation (blue).

Enrichment in Differentially Expressed Genes

Gene ontology (GO) enrichment analysis was performed to classify the DEGs function. According to the database annotation, DEGs can be categorized into 220 functional groups consisting of three domains: “biological process,” “cellular component,” and “molecular function.” Top 10 enrichment groups are shown in Figure 3A. The GO enrichment results of all DEGs are significantly enriched in “extracellular matrix organization (GO:0030198)”, “extracellular structure organization (GO:0043062)” and “urogenital system development (GO: 0001655).”

Figure 3 GO enrichment analysis and KEGG pathway enrichment of DEGs. (A) Enrichment analysis of gene ontology terms. (B) Enrichment analysis of KEGG pathway in up-regulated genes and down-regulated genes. Top gene sets in GSEA in sensitive (C) and resistant (D) groups.

Kyoto Encyclopedia of Genes and Genomes (KEGG) database was used as an alternative approach to categorize gene function with emphasis on biochemical pathways. In total, 48 up-DEGs and 115 down-DEGs were collected by KEGG (Figure 3B). DEGs were assigned to 9 and 5 KEGG pathways in up and down genes. Most of the up-regulated DEGs were enriched in “Cytokine-cytokine receptor interaction” (7, 14.58%), “Complement and coagulation cascades” (7, 14.58%). Down-regulated DEGs are mainly enriched in “Focal adhesion” (10, 8.70%), followed by “cGMP-PKG signaling pathway” (9, 7.83%).

Gene Set Enrichment Analysis in Resistant and Sensitive Groups

Gene Set Enrichment Analysis (GSEA) is a conventional approach to identify pathways related to gene expression. After the filters, we chose 10 terms in resistant and sensitive groups which were associated with our study topic in enriched terms (supplement 1). Then, the top 50 enriched terms in each group were analyzed by GSEA automatically (supplement 2). After that, we compared the top automatic terms and our chosen terms and picked four similar terms in each group from these two ways for further analysis. The results revealed that sensitivity was associated with proliferation, mitochondria functions, liver specific genes, relative lower hypoxia (Figure 3C), while resistance was related to collagen expression, EMT phenotype, stem cell phenotype, and TGF beta signaling (Figure 3D).

Interaction of EGFR-Related Gene Sets and Gene Expression in Resistant and Sensitive Groups

Nimotuzumab was a classical anti-EGFR monoclonal antibody, so we wanted to explore the interaction of EGFR-related gene sets and gene expression in resistant and sensitive cell lines. According to the GSEA analysis, we found EGFR tyrosine kinase inhibitor resistance and Kobayashi EGFR signaling 24HR, two EGFR-related gene sets, revealed to sensitive group genes besides the traditional EGF/EGFR signaling pathway gene set. After that, we analyzed gene expression of sensitive and resistant cell lines in these three pathways (Figure 4A–C). However, these cell lines could not be clustered in two groups by cluster parameters.

Figure 4 Interaction of EGFR-related gene sets and Hoshida subclass of HCC with gene expression in sensitive and resistant groups. (A) EGF/EGFR signaling pathway. (B) Kobayashi EGFR signaling 24HR. (C) EGFR tyrosine kinase inhibitor resistance. (D) Hoshida subclass of HCC.

Sensitivity of Nimotuzumab Might Associate with Hoshida Subclass of HCC

Hoshida subclass of HCC was a class with different AFP level and mechanism of carcinoma.²⁹ Since one enriched gene term of sensitive cell lines from GESA above was Hoshida subclass S3, we tested the gene expression of six cell lines in Hoshida subclass genes (Figure 4D). It was interesting that two sensitive cell lines showed a significant high expression in S2 & S3 subclass genes and low expression in S1 subclass genes. Meanwhile, four resistant cell lines revealed a high expression in S1 subclass genes and low expression in S2 & S3 subclass genes. According to these results, we inferred that differences between mechanisms of three Hoshida subclasses might play important roles in the efficacy of nimotuzumab to HCC.

Differentially Expressed Genes Between Hoshida S1 and S2&S3 in Human Tumor Samples

We found that Hoshida subclass may be a key to the efficacy of nimotuzumab to HCC. Thus, we wanted to explore the DEGs in real human tumor samples between Hoshida S1 and S2 & S3. Data of array from GSE102079 which were from the Gene Expression Omnibus (GEO). Data of RNA sequences were from The Cancer Genome Atlas (TCGA) dataset (TCGA-LIHC). We filtered all data and kept HCC samples. After that, we divided all tumor samples to Hoshida S1, S2 and S3 according to gene expression in every subclass gene set. Counts in S1/S2/S3 subclass were 42/37/73 in GSE102079 and 93/112/169 in TCGA-LIHC.

We got DEGs between S1 versus S2 and S1 versus S3 (Figure 5). Array data from GSE102079 were calculated by limma. RNA-seq data from TCGA-LIHC also was analyzed from intersection of DESeq2, edgeR and limma voom. Then we intersected DEGs of S1 versus S2 with DEGs of S1 versus S3. After that, we compared DEGs between GSE102079 and TCGA-LIHC. Thus, we got 20 down-related genes after these procedures (supplement 3). We listed all these 20 gene-related pathways and EGFR-related pathways and found four genes had the same pathways with EGFR. These four genes were COL1A1, COL1A2, COL3A1 and MMP9, which might relate to efficacy of nimotuzumab in HCC.

Figure 5 DEGs between Hoshida S1 and S2&S3 in GSE102079 and TCGA-LIHC. (A) Venn diagram in up-regulated genes. (B) Venn diagram in down-regulated genes. (C) DEGs were analyzed from GSE102079. (D) DEGs were analyzed from TCGA-LIHC.

Discussion

Hepatocellular carcinoma (HCC) continues to be a cancer with rising incidence, high mortality, and recurrence rate. The therapeutic effects on HCC are not satisfactory currently. EGFR is an important factor in HCC and previous studies have demonstrated that EGFR is involved in various HCC mechanisms including cell proliferation, late-stage metastasis, treatment and drug resistance.^30–35 However, anti-EGFR agencies do not show ideal results in HCC involving extremely complicated pathways.^10,12 Lin et al. reported that the anti-EGFR drug resistance in HCC is caused by the promoted interaction of EGFR with mTORC2.¹³ Zhang et al. found that CD317 activates EGFR in resistant progress as well.¹⁴ To explore the resistance mechanism in depth and apply the possibility of EGFR targeted drugs, our research analyzed six HCC cell lines (Hep 3B2.1–7, Li-7, PLC/PRF/5, SK-HEP-1, SNU-182, SNU-387) for sensitivity to nimotuzumab. Then we compared the RNA-seq differences between sensitive and resistant group cell lines. We enriched DEGs in GO and KEGG and performed GSEA in each group. Genes in the two groups did not show obvious different expressions in EGFR pathways, while Hoshida subclass of HCC seemed to associate with efficacy of nimotuzumab. Therefore, we analyzed genes in real human tumor samples which were from TCGA (TCGA-LIHC) and GEO (GSE102079). We found that COL1A1, COL1A2, COL3A1, and MMP9 were the focus DEGs of S1 and S2&S3 related to EGFR.

We discovered that all HCC cell lines had positive expression of EGFR protein on cell surface. It was also reported that EGFR was expressed in most HCC previously.^36,37 Expression of EGFR protein had no correspondence with efficacy of nimotuzumab. Higher EGFR expression did not equal to better effect. It was also observed in triple-negative breast cancer (TNBC) that TNBC tumors express EGFR but the clinical efficacy of anti-EGFR therapy in TNBC is low.³⁸ That might relate to the RTK family receptors which are generally thought to be activated by ligand binding promoting homodimerization. Dimer ligand may consist of two same RTKs (e.g., EGFR/EGFR) or two different RTKs (e.g., EGFR/FGFR) on the cell surface.³⁹ Therefore, EGFR expression alone could not reflect the real situation of EGFR on the cell surface. Moreover, our research also found that EGFR-related pathway gene sets did not express differences between sensitive and resistant cell lines. These results revealed that EGFR might connect with other pathways in cancer cell proliferation, so other survival pathways existed under anti-EGFR drug. Therefore, further study of this aspect is needed.

After GO and KEGG enrichment were used in DEGs between resistant and sensitive groups, we performed GSEA in each resistant and sensitive group. At that time, we found

Hoshida subclass of HCC might be a key to efficacy of nimotuzumab. Two sensitive group cell lines showed a high expression in S2 & S3 genes, which proved that these two cell lines Hep 3B2.1–7 and PLC/PRF/5 were Hoshida S2. Resistant group showed a high expression in S1 genes that also proved that SNU-182, SNU-387 and SK-HEP-1 were Hoshida S1.²⁹ Li-7 was uncertain in past studies. Different subclasses had different mechanisms of HCC progress in that S1 reflected aberrant activation of the WNT signaling pathway, S2 was characterized by proliferation as well as MYC and AKT activation, and S3 was associated with hepatocyte differentiation.²⁹ Different mechanisms of HCC progress might lead to different efficacy of nimotuzumab. We should explore the advanced association between the subclass mechanisms and anti-EGFR drug.

After analyzing the DEGs in real human tumor samples between S1 and S2 & S3, we got four genes, COL1A1, COL1A2, COL3A1, and MMP9. COL1A1, COL1A2, and COL3A1 are all extracellular matrix (ECM) genes and MMP9 is the downstream gene of EGFR. In addition to functioning as structural scaffold elements, ECM proteins can serve a wide array of signaling functions especially during fibrogenesis and fibrosis resolution which results in liberation and activation of growth factor signaling cascades and cancer development in injured livers.⁴⁰ MMP9 is involved in cancer progression and pre-metastatic niche formation in HCC.⁴¹ There are many studies that show some targeted factors like cancer-associated fibroblasts (CAFs) and M2 macrophages can influence MMP9 expression to promote HCC development.^42,43 The detailed relation of anti-EGFR treatment and ECM and MMP9 should be researched in the future.

HCC cell lines might have a difference in vivo and our study were not concluded using large amounts of data. These limitations may affect our study results. However, the results of our study may provide a new direction of anti-EGFR drug in HCC. We could collect more human HCC tissues to explore the connection of genes and nimotuzumab especially in Hoshida subclass in the future. Meanwhile, ECM genes and MMP9 should be researched in relation with efficacy of nimotuzumab in HCC.

Conclusions

Li-7, SK-HEP-1, SNU-182, and SNU-387 cell lines were resistant to nimotuzumab, while Hep 3B2.1–7 and PLC/PRF/5 cell lines were sensitive. The efficacy of nimotuzumab did not show direct relevance with EGFR protein expression and EGFR-related pathway. However, efficacy could be associated with Hoshida subclass of HCC in that S2 & S3 showed better therapeutic effect than S1. After we compared DEGs between S1 and S2 & S3 related to EGFR, three ECM genes (COL1A1, COL1A2, COL3A1) and MMP9 were focused on, which might play important roles in curative effect of nimotuzumab in HCC.