Transcriptome changes in ERGIC3-knockdown hepatocellular carcinoma cells: ERGIC3 is a novel immune function related gene

Cell culture

Human hepatocellular carcinoma SMMC-7721 cells, obtained from the Shanghai Cell Bank, Chinese Academy of Sciences, were cultured in complete RPMI 1640 medium (Hyclone, USA) containing 1% penicillin/streptomycin and10% fetal bovine serum (Jiangsu Enmoasai Biotechnology Co., LTD, China) at 5% CO₂, 37 °C and saturated humidity and subcultured for 2–3 days.

Experimental grouping and liposome transfection

HCC cells were randomly divided into control group and ERGIC3i group, each group had three biological replicates. One day before transfection, 7.5  × 10⁵ cells per well were seeded in a six-well plate. When the confluence of cells was about 70–80%, Negative Control #1 siRNA and ERGIC3 siRNA (Thermo Fisher Scientific, Waltham, MA, USA) were transfected by Lipofectamine 3000 (Invitrogen, Waltham, MA, USA) into the control group and ERGIC3i group, respectively. The transfection procedures were performed according to the instructions. After 72 h of culture, cell samples were collected for subsequent experiments.

Total RNA extraction, Q-RT-PCR and library establishment

Total RNA was extracted by TRIzol method, the integrity of total RNA was detected by Gel Doc XR Gel imaging system (Bio-Rad, Hercules, CA, USA), the concentration and purity of total RNA were determined by Nanodrop 2000 ultra-micro spectrophotometer (Thermo Scientific, USA). The cDNA synthesis was performed by according to its specification. Quantitative Real-Time PCR (qRT-PCR) was performed with CFX96 Touch™ fluorescent quantitative PCR instrument (Bio-Rad). The PCR reaction procedure was as follows: 98 °C 30s, 98 °C 5s, 60 °C 5s, 40 cycles. The relative gene expression was calculated by 2^−ΔΔCT (Schmittgen & Livak, 2008). The primers sequence were as follows: GAPDH, forward: 5′GGCCTCCAAGGAGTAAGACC3′, reverse: 5′AGGGGAGATTCAGTGTGGTG3′; ERGIC3, forward: 5′GGAGAGGTACTGAGGACAAATCA3′, reverse: 5′AGCTCATAGAGGACGAAGACTC3′.

Cells in each group were collected, and an appropriate amount of RNAiso Plus lysate (TakaRa Biotechnology, Beijing, China) was added and transported to Shenzhen BGI Co., Ltd on dry ice. Total RNA extraction and quality control, mRNA enrichment, library construction and quality control were completed by Shenzhen BGI Co., Ltd.

Western blot

The cells were collected, added with RIPA lysis buffer (Solarbio Lifesciences, Beijing, China) and protease inhibitor (BBI, Shanghai, China), and lysed on ice for 10–15 min. The total proteins were obtained by centrifugation, and the protein concentration was determined by BCA method. The proteins were boiled and denaturated, and then separated by SDS-PAGE electrophoresis. The primary antibody (GAPDH (ProteinTech Group, Wuhan, China), 1:10000 dilution, ERGIC3 (Abcam, Cambridge, UK), 1:1000 dilution) was incubated overnight at 4 °C, and the secondary antibody was incubated at 37 °C for 2 h, ECL illuminant was added, and exposure and development were performed.

Sequencing analysis

Firstly, clean reads were obtained by filtering out low-quality, adaptor-polluted and high content of unknow base (N) reads with SOAPnuke v1.5.2, a filtering software independently developed by BTU. After reads filtering, we mapped those clean reads onto human reference genome GRCh38 using HISAT2 v2.0.4, followed with novel gene prediction. Finally, we identified differentially expressed genes between samples and did clustering analysis and functional annotations.

Bioinformatics analysis of differentially expressed genes

The correlation analysis between ERGIC3 expression and the level of immune cell infiltration

The association between ERGIC3 expression and the abundance of tumor immune infiltrates (B cells, CD8+ T cells, CD4+ T cells, macrophage, neutrophil and myeloid dendritic cells) was analyzed via TIMER (http://timer.cistrome.org/) (Li et al., 2020).

Statistical analysis

The GraphPad Prism 8.3.0 was used for statistical analysis of the qRT-PCR and Western Blot data. The t-test was used to compare the mean of two independent samples for intra-group analysis. One-way analysis of variance was used for multi-group comparison. P < 0.05 was considered statistically significant.

The high expression of ERGIC3 in HCC is associated with survival prognosis

Studies have shown that the abnormally high expression of ERGIC3 is associated with the proliferation and growth of HCC cells and epithelial-mesenchymal transition (EMT) (Zhang et al., 2013). Further, GEPIA (http://gepia.cancer-pku.cn/detail.php) database (—log₂FC ∣cutoff = 1, p-value cutoff = 0:01) and UALCAN (http://ualcan.path.uab.edu/analysis.html) online analytical tools were selected to identify the difference of ERGIC3 expression, survival prognosis and its relationship with clinical features between hepatocellular carcinoma patients and normal controls. The results showed that ERGIC3 was highly expressed in HCC (Fig. 1A), and the high expression level of ERGIC3 was associated with poor overall survival rate of HCC patients. The five-year survival rate was about 27.5% (Fig. 1B). In addition, ERGIC3 is more expressed in 61–80 years old, early stage, intermediate T grades, and tumors with No regional lymph node metastasis (Fig. 1C).

ERGIC3 was successfully knocked down in human hepatoma SMMC-7721 cells

In order to further reveal the molecular mechanism of ERGIC3 promoting the growth of liver cancer cells, ERGIC3 was knocked down in SMMC-7721 cells by RNAi technology. The results of Q-RT-PCR and Western blot analysis demonstrated that ERGIC3 levels was significantly decreased in SMMC-7721 cells after transfection with ERGIC3 siRNA for 24 h, 48 h and 72 h, respectively (Fig. 2). Then transcriptome sequencing technology was used to obtain the differentially expressed genes between the control group and the ERGIC3i group. Finally, the possible molecular mechanisms were analyzed by bioinformatics.

Sequencing data analysis

RNA sequencing of 6 samples were performed by the Illumina HiSeq platform, and the average output of each sample was 6.90 GB of data. The original transcriptome sequencing data have been submitted to NCBI, the accession number is PRJNA776533. The Q20(%) of clean reads was greater than 98%, and the Q30(%) of clean reads was greater than 95%. Moreover, the average comparison rate of samples for genomes was 85.42% as well as the average comparison rate of the gene set was 71.87%. The above results indicated that it was high sequencing quality (Table 1), and the data could be next step.

In order to further validate the accuracy of the sequencing results, we randomly selected five genes for Q-RT-PCR analysis. The expression trend of these genes was basically consistent with the transcriptome sequencing results (Fig. 3 and Table 2), which indirectly indicates that the transcriptome sequencing data have a certain credibility.

Screening and analysis of differentially expressed genes

The differentially expressed genes (DEGs) were defined with FDR<= 0.001 and the difference ratio was more than two fold change. A total of 210 differentially expressed genes were screened, including 192 annotated genes and 18 unknown genes. The expression of each differential gene was plotted on the horizontal and vertical axes according to the values between different groups, namely, as well as the scatter plot of the differential genes was obtained (Fig. 4A). According to the distribution of DEGs, the expression levels of 176 genes were up-regulated and 34 genes were down-regulated. These results indicated that ERGIC3 knockdown resulted in differential expression of many genes in SMMC-7721 cells, suggesting that the promotion of ERGIC3 to HCC cells growth is a complicated process. At the same time, a total of 20 intergroup differential genes were obtained by Venn analysis (Fig. 4B), among which nine genes were up-regulated and 11 genes were down-regulated. In order to more visually display the differentially expressed genes between the control group and the ERGIC3i group, a heatmap was made for each group of DEG (Fig. 4C).

Go and Pathway enrichment analysis of differentially expressed genes

Gene Ontology (GO) and KEGG database were used to annotate the functions and pathways involved in the differentially expressed genes.

It (Fig. 5A) showed that, in the biological process, the most DEGs were involved in biological regulation, metabolic process, signaling, localization, cellular component organization or biogenesis, immune system process, locomotion, biological adhesion, growth. In the cellular component, the DEGs were enriched in macromolecular complex, membrane-enclosed lumen, cell junction, cytoskeleton. Significantly, in the molecular function, the DEGs involved in molecular transducer activity, transporter activity, nucleic acid binding transcription factor activity, transcription factor activity, protein binding and chemoattractant activity.

In the form of directed acyclic graph (DAG), the distribution of enrichment items of different genes at three levels of GO is concentrated, which are biological process (Table S1), cellular component (Table S2) and molecular function (Table S3). The pathway with FDR <=0.01 was considered as significant enrichment. Differentially expressed genes were enriched in the biological process group in regulation of leukocyte migration, positive regulation of leukocyte migration, regulation of signaling, cell surface receptor signaling pathway, positive regulation of leukocyte chemotaxis. The integral component of membrane, extracellular region, extracellular region part, cell periphery, vesicle are significantly enriched in cellular component. Growth factor receptor binding, vascular endothelial growth factor receptor binding, cytokine receptor binding, collagen binding, receptor binding are significant enrichment items in molecular function group.

The significantly enriched pathways were analyzed according to KEGG functional annotation results of differentially expressed genes. The pathway with FDR <= 0.01 was regarded as significant enrichment, and only the top 20 items were shown (Fig. 5B). The enrichment pathway mainly included the PI3K-Akt signaling pathway, focal adhesion, cytokine-cytokine receptor interaction, microRNAs in cancer, phagosome, regulation of actin cytoskeleton, natural killer cell mediated cytotoxicity, chronic myeloid leukemia, NOD-like receptor signaling pathway, Jak-STAT signaling pathway, NF-kappa B signaling pathway. The KEGG functional significance enrichment entries of differentially expressed genes were consistent with the positive regulation of leukocyte migration and chemotaxis, cell surface receptor signaling pathway of biological process components and the growth factor receptor binding, vascular endothelial growth factor receptor binding, cytokine receptor binding, collagen binding of molecular functional components in the GO functional significance enrichment items. It indicated that the related genes regulated by ERGIC3 played an important role in PI3K-Akt, NOD-like, Jak-STAT, NF-kappa B and other protein kinase-coupled receptors mediate signal transduction pathways, tumor nonspecific immune response, collagen-integrin receptor-actin axis, and miRNA pathways. Interestingly, in the top 10 of KEGG pathways of the DEGs, 8 pathways were link to immune system and infectious diseases. The other two were cytokine-cytokine receptor interaction, NF-kappa B signaling pathway, which also played enssential roles in the immune response (Table 3).

Protein–protein interaction networks of DEGs

Proteins usually interact with each other to form complexes that perform specific functions. It is the basis of PPI software and DEGs with interaction usually have similar functions. In order to explore the interaction between 210 differentially expressed genes, a PPI network (Fig. 6A) was constructed using the STRING protein interaction database and Cytoscape software, which revealed 120 nodes and 303 edges in total. Furthermore, MCODE plug-in was used for module analysis, and modules with node number greater than four and MCODE score greater than four were selected to obtain two important modules (Figs. 6B and 6C). It showed that, from the proteins related to modules, the proteins were mainly involved in tumor immune response, integrin gene family and growth factors. Subsequently, in the up-regulated genes, ten hub genes with high connectivity were selected to analysize (Fig. 6D and Table 4). Among these genes, IL6, CD44, and CXCL8/IL8 were the most significant genes with connectivity degrees of 38, 32, 28, respectively. CXCL8 and CD44 were the core in the regulatory network.

Survival analysis of hub genes

In order to investigate the prognostic value of the hub genes, Kaplan–Meier was used to analyze the survival data of the 10 hub genes. The results showed that the high level of CXCL8/IL8 was associated with poor overall survival. Patients with high expression of TLR4 and VEGFC had longer overall survival than those with low expression (Fig. 7), while the association of the other seven genes was not significant changes.

ERGIC3 may be an immune-related key gene

Considering that ERGIC3 expression caused significant changes in the interleukin gene family, we speculated that ERGIC3 might be an immune-related key gene. Then we use the TIMER (http://timer.cistrome.org/) algorithm (Li et al., 2020) to evaluate the correlation between ERGIC3 expression and the level of immune cell infiltration. The results demonstrated that the expression of ERGIC3 was positively correlated with the infiltration levels of B cells (r = 0.311, P < 0.05), macrophage (r = 0.238, P < 0.05) and myeloid dendritic cells (r = 0.436, P < 0.05). The expression of ERGIC3 was weakly positively associated with infiltration levels of CD4⁺ T cells (r = 0.177, P < 0.05), but not CD8⁺ T cells (r = 0.051, P > 0.05) and neutrophil (r = 0.027, P > 0.05). Moreover, there was a weak positive correlation between the expression of ERGIC3 and tumor purity (r = 0.172, P < 0.05) (Fig. 8). These results suggested that ERGIC3 was closely related to immune cell infiltration of HCC.