THRSP identified as a potential hepatocellular carcinoma marker by integrated bioinformatics analysis and experimental validation

Identification of crucial modules by WGCNA

The “WGCNA” package was used to group genes into modules by the average linkage hierarchical clustering. In this study, the soft powers β = 3 and 5 were selected as the soft-thresholding to ensure scale-free networks (Figure 1A, 1B), and 11 modules in the TCGA-LIHC (Figure 1C) and 8 modules in the GEO datasets (Figure 1D) were generated. The heatmaps (Figure 1E, 1F) of module-trait relationships were plotted to identify modules most significantly correlated with clinical features (normal and tumor). We found the brown modules in the TCGA-LIHC (containing 2057 co-expression genes) and GEO datasets (containing 2145 co-expression genes) had the highest association with tumor tissues (brown module in TCGA-LIHC: r = 0.69, p = 7e-60; brown module in GEO datasets: r = 0.86, p = 4e-56), which were selected as modules of interest for the subsequent analysis.

Figure 1. Identification of co-expression modules associated with the clinical traits in the TCGA-LIHC dataset and two GEO datasets. (A, B) Analysis of the scale-free fit index for various soft-thresholding powers (β). (C, D) Dendrogram of all genes in the TCGA-LIHC dataset or GEO datasets clustered based on the 1-TOM matrix. (E, F) Correlation between modules and clinical traits (normal and tumor). Each cell contains the corresponding correlation coefficient (the upper number) and the P-value (the lower number).

Identification of DEGs and differentially co-expressed genes

After normalization of the microarrays, 2705 differentially expressed genes (DEGs) between the HCC and normal tissues from the TCGA dataset (TCGA_diff) and 567 DEGs from the GEO datasets (GEO_diff) were screened by the “limma” package in R. Then, a Venn diagram was performed to examine the intersection among the DEGs and co-expressed genes of key modules. As shown in Figure 2A, 60 differentially co-expressed genes were finally obtained.

Figure 2. Identification of hub genes. (A) The Venn diagram for selection of the differential co-expression genes among DEG lists and co-expression modules. (B) PPI network of the intersection genes between DEG lists and co-expression modules. Each blue node represents a gene. Edges among nodes indicate interaction associations between genes. (C) Identification of the core genes from the PPI network by MCC algorithm. Darker colors refer to higher MCC sores. (D) The top three significant modules of the PPI network were evaluated in MOCDE. Pink nodes represent genes in corresponding modules.

PPI network and hub genes

The PPI network of the 60 differentially co-expressed genes was constructed in the STRING database (Figure 2B). As showed in Figure 2C, the core genes were extracted from the PPI network by the MCC algorithms via the CytoHubba plug-in. Meanwhile, the significant modules of the 60 differential co-expression genes were established using the MCODE plug-in in Cytoscape (Figure 2D). Combining the above two algorithms, AFP, IGF1, BCHE, ACSM3, LCAT, ACSL4, ACADS, ENO3, CYP1A2, THRSP, GSTZ1, and SRD5A2 were finally selected as hub genes. Among them, AFP and ACSL4 were up-regulated and the other 10 hub genes were down-regulated in HCC tissues (Figure 3). To evaluate the prognostic values of the hub genes in HCC patients, a 5-year survival analysis was performed by Kaplan-Meier plotter. As shown in Figure 4, higher expression of LCAT, ACSM3, IGF1, SRD5A2, THRSP, and ACADS was associated with a better 5-year overall survival (p<0.05).

Figure 3. Validation of expression levels of the 12 hub genes in HCC and normal tissues using GEPIA. *P<0.01 is considered statistically significant. Tumor tissues are shown in red, and normal liver tissues are shown in blue.

Figure 4. Overall survival (OS) Kaplan-Meier curves of the 12 hub genes in HCC patients. (A) AFP, p = 0.780. (B) IGF1, p = 0.008. (C) BCHE, p = 0.108. (D) ACSM3, p = 0.003. (E) LCAT, p <0.001. (F) ACSL4, p = 0.709. (G) ACADS, p = 0.045. (H) ENO3, p = 0.156. (I) CYP1A2, p = 0.353. (J) THRSP, p = 0.013. (K) GSTZ1, p = 0.089. (L) SRD5A2, p = 0.016.

THRSP was down-regulated in HCC tissues and cell lines

To verify the expression of THRSP in HCC tissues and cell lines, the quantitative real-time PCR (qRT-PCR) and immunohistochemistry (IHC) assays were performed in this study. The results of qRT-PCR verified that THRSP mRNA expression was remarkably lower in HCC tissues as well as in HCCLM3 and Huh-7 cell lines compared with normal groups (Figure 5A, 5B). The results of IHC indicated that THRSP was down-regulated in most of HCC samples (8/10) (Figure 5C).

Figure 5. The expression of THRSP at mRNA and protein levels. (A, B) The mRNA expression of THRSP in HCC tissues and cells compared with control examined by RT-qPCR. (C) The protein expression of THRSP in HCC tissues and the adjacent normal tissues examined by immunohistochemical (200×).

The correlation between THRSP expression and clinical characteristics

Based on data from Illumina HumanMethylation 450 platform, we found a significant negative correlation (r = −0.81, p <0.001) between THRSP mRNA expression and methylation (Supplementary Figure 1A). Four THRSP promoter CpG sites were shown in Supplementary Figure 1B. The Spearman correlation analysis demonstrated that CpG sites were negatively correlated with the expression of THRSP (Supplementary Figure 1C–1F). Wilcox test was used to compare the difference in THRSP expression or methylation in groups divided by age, family history of cancer, gender, grade, Ishak fibrosis score and race. The details were shown in Supplementary Figures 2, 3. Then, the HCC patients were dichotomized into two groups (“low” or “high”) based on their THRSP expression levels or THRSP methylation levels using the median values as the cutoff point. The chi-square test was used to evaluate the correlation of THRSP expression or THRSP methylation with clinical characteristics. As listed in Table 1, THRSP mRNA expression or methylation was closely associated with the clinical indicators including age, family history of cancer, Ishak fibrosis score and gender et al.

Relationship of THRSP with immune infiltration based on TIMER 2.0

As illustrated in the scatter plots (Supplementary Figure 4), THRSP expression was negatively correlated with infiltrating levels of B cells, CD4+ T cells, dendritic cells, and positively correlated with CD8+ T cells, but was not correlated with macrophages or neutrophils. In addition, we investigated the correlations of THRSP expression with the gene markers of various immune immune-infiltrating cells (including T cell, B cell, monocyte, neutrophil, dendritic cell, tumor-associated macrophage and different types of functional T cells) in HCC with the TIMER 2.0 database. As listed in Supplementary Table 1, THRSP expression was negatively associated with most of the gene markers.

Silencing of THRSP promoted HCC progression

To explore the functional roles of THRSP in HCC, three candidate siRNAs of THRSP (THRSP-siRNAs) were transfected into HCCLM3 and Huh-7 cells. qRT-PCR was performed to evaluate the inhibition efficiency of the three THRSP-siRNAs. Three THRSP-siRNAs significantly inhibited THRSP expression compared with the negative control (Figure 6A). The si-THRSP-3 was selected for the following experiments due to its robust silencing efficiency. The CCK-8 assays indicated that silencing of THRSP could promote proliferation of HCCLM3 and Huh7 cells (Figure 6B, 6C). The flow cytometry assays indicated that the si-THRSP group presented a significantly higher percentage of HCCLM3 cells in the S phase and the G2/M phase, but a lower percentage of cells in the G0/G1 phase compared with the control group (Figure 6D–6F). Transwell and wound scratch assays indicated that silencing of THRSP also promoted migration and invasion of HCC cells (Figure 7A–7C).

Figure 6. The effect of THRSP expression on proliferation and cell cycle of HCC cells. (A) The efficiency of the THRSP silencing determined by qRT-PCR. (B, C) The proliferation of HCCLM3 and Huh-7 cells examined by CCK-8. (D–F) The cell cycle assay detected by flow cytometry. NC group: black bars; si-THRSP group: grey bars. *p<0.05, **p<0.01, ***p<0.001.

Figure 7. Transwell and wound healing assays. (A, B) The migration and invasion of HCCLM3 and Huh-7 cells detected by Transwell assays (magnification: 200×). (C) The migration ability of HCCLM3 and Huh-7 cells examined by wound healing assays (magnification: 40×).

Silencing of THRSP inhibited apoptosis of HCC cells

In this study, cells were strained with Annexin V/PI and subjected to flow cytometry to determine the apoptotic cells. The results indicated that silencing of THRSP could inhibit HCC cell apoptosis (Figure 8A, 8B). qRT-PCR and Western blotting assays were performed to detect the cell apoptosis-related molecules, including bax, bcl-2 and caspase 3. The results indicated that silencing of THRSP significantly reduced the expression of bax and caspase 3, while enhanced the expression of bcl-2 (Figure 8C, 8D).

Figure 8. The effect of THRSP expression on cell apoptosis and channel regulation in HCC. (A, B) The apoptosis of HCCLM3 cells detected by flow cytometry. (C–F) The expression of apoptosis-related molecules (bax, bcl-2 and caspase 3) and MAPK/NF-κB pathway-related molecules (erk, p-erk, p38, p-p38, p65 and p-p65) examined by qRT-PCR or Western blotting assays. (G) Each of the protein bands. The dividing lines (dashed lines) indicated that the grouping of images were from different parts of the same gel. (H–J) Invasion and migration rates were analyzed when treated with AG-126, SB203580, SP600125 and BAY-11-7082 inhibitors.

THRSP regulated HCC cell progression by modulating MAPK/NF-κB signaling pathway

To further understand the molecular mechanism by which si-THRSP promoted the migration and invasion of HCC cells, we explored the potential signaling pathways including the NF-κB and MAPK signaling pathways by Transwell and wound healing assays. As shown in Figure 8H–8J, compared with negative control, the migration and invasion of HCC cells transfected with si-THRSP or si-NC were inhibited after being treated with BAY-11-7082 (NF-κB inhibitor), AG-126 (ERK1/2 inhibitor), SB203580 (p38 MAPK inhibitor) and SP600125 (JNK inhibitor). And there was no difference between the si-THRSP group and the si-NC group when the cells were treated with BAY-11-7082, AG-126 and SB203580, indicating that the cell migration and invasion caused by the siRNA-induced silencing of THRSP might depend on the NF-κB, ERK1/2 and p38 MAPK signaling pathways. qRT-PCR and Western blotting assays were performed to assess the expression of MAPK/NF-κB pathway-related molecules (p65, p-p65, p38, p-p38, erk1/2 and p-erk1/2). The results showed that the silencing of THRSP increased the phosphorylation of ERK1/2, p38 MAPK and p65, and enhanced the expression of p65, p38, and erk1/2 (Figure 8E, 8F). Each of the protein bands were showed in Figure 8G.

Datasets preprocessing

Two gene expression profiles GSE84005 and GSE121248 were obtained from the GEO database. The GSE84005 dataset, including 38 tumor tissues and paired 38 normal tissues from HCC patients, was based on the GPL5175 platform (Affymetrix Human Exon 1.0 ST Array). The GSE121248 dataset, including 70 HCC tissue samples and 37 adjacent non-tumor tissue samples, was based on the GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array). Two datasets were merged by Perl 5.3 (available online: http://www.perl.org/) to increase the sample size. The merged dataset was batch-normalized by “limma” and “sva” packages in R 3.6.3 (https://www.r-project.org/) to remove batch effects.

Besides, the RNA-sequencing (RNA-seq) data of 373 HCC and 49 normal samples and the corresponding clinical information was downloaded from the TCGA database. As recommended by the package “edgeR” in R, genes with low read counts (count per million (cpm) ≤ 1) were omitted. Gene expression was calculated and normalized to RPKM (reads per kilobase per million) values using function “rpkm” in the “edgeR” package. Moreover, DNA promoter methylation data (Methylation 450k, including 430 samples) was downloaded from the TCGA database via the UCSC Xena browser (https://xenabrowser.net/).

WGCNA analysis

The “WGCNA” package in R was applied to construct co-expression networks and to explore the key modules of highly relevant genes by cluster analysis for relating modules to sample traits. In this study, the gene expression profiles of TCGA-LIHC and the merged dataset of GSE84005 and GSE121248 were respectively used to construct WGCNA. Briefly: an adjacency matrix was created by Pearson’s correlations between each of the gene pairs. Next, the adjacency matrix was utilized to erect a scale-free co-expression network based on the soft threshold power β which was selected using the pickSoftThreshold function [46]. Subsequently, the adjacency matrix was converted into a topological overlap matrix (TOM) as well as the corresponding dissimilarity (1-TOM). Then, module identification was conducted using the dynamic tree Cut approach by average linkage hierarchical clustering based on the TOM-based dissimilarity measure with the parameters of minModuleSize of 50, deepSplit value of 2, and mergeCutHeight of 0.25 for the genes dendrogram. Afterward, the correlation between module eigengenes (MEs) and the clinical trait information was calculated by the module-trait relationship analysis of WGCNA to identify the clinically significant modules in a co-expression network. Finally, modules with a high correlation coefficient were selected for further analysis.

Identification of differentially co-expressed genes

The “limma” package in R was used to filter the DEGs between the HCC and normal samples in the TCGA-LIHC and the merged dataset of GSE84005 and GSE121248. The adjusted P-value (adj. P) < 0.05 and |log2foldchange (FC)| >1 was set as the criteria of DEGs. Then, the overlapping genes between DEG lists and co-expression genes from significant modules were screened out using the “VennDiagram” package in R.

Identification of hub genes

The PPI network was constructed in the STRING (http://string-db.org) and visualized by the Cytoscape software (Cytoscape_v3.8.0, https://cytoscape.org/) [47]. The maximal clique centrality (MCC) analysis was performed to extract the candidate hub genes with the top20 MCC values in the PPI network using the cytoHubba plug-in [48]. Meanwhile, the Molecular Complex Detection (MCODE) [49] plug-in of Cytoscape was implemented to find significant PPI modules with degree cut-off ≥2, node score cut-off ≥0.2, K-core ≥2, and max. depth =100. Finally, the overlapping genes obtained from the MCC analysis and MCODE analysis were regarded as hub genes.

The gene expression and prognostic analysis of the hub genes

The differential expression analysis of the hub genes between HCC and normal tissues was performed based on the GEPIA2 database. The 5-year overall survival (OS) analysis of these hub genes was performed using the “survival” and “survminer” packages in R [50].

The correlation between gene expression and clinical features

According to the median value of gene expression or methylation, the HCC patients were divided into the low- or high- group. The chi-square test was utilized to investigate the correlation of gene expression as well as methylation with clinical characteristics. The correlation between DNA methylation and gene expression in HCC samples from the TCGA database was examined using the Spearman correlation coefficient and visualized by “ggplot2” and “ggpubr” packages in R.

Immune infiltrate analysis

The online tool TIMER 2.0 (https://cistrome.shinyapps.io/timer/) was used for immune infiltrate analysis [51]. The abundance of six types of immune cells (including CD8+ T cells, CD4+ T cells, B cells, neutrophils, macrophages, and dendritic cells) were computed by TIMER algorithm. In addition, we exploited the correlation between gene expression and the gene markers of different kinds of immune cells. The immune gene markers of interest used in this study were referred to the previous studies [52–56].

Cell culture

The human normal liver cell line (LO2) and the human hepatoma cell lines (HCCLM3 and HUH-7) were purchased from Procell Life Science and Technology Co., Ltd. (Wuhan, China). The cells were cultured in DMEM (Servicebio Technology Co., Ltd, Wuhan, China) supplemented with 10% Fetal Bovine Serum (FBS, G-CLONE, Beijing, China) and maintained in a incubator with 5% CO2-humidified atmosphere at 37° C.

Total RNA extraction and qPCR

The quantitative real-time polymerase chain reaction (qRT-PCR) was employed to detect the expressions of THRSP. GAPDH was served as a reference gene. Total RNA was extracted from tissues and cells by the TRIzol reagent (G-CLONE, Beijing, China). The expression of THRSP was determined by the SweScript RT I First Strand cDNA Synthesis Kit with gDNA Remover and the SYBR Green qPCR Master Mix (High ROX) (Servicebio Technology Co., Wuhan, China) according to the manufacturer’s protocol. The qRT-PCR was performed on the StepOne Plus Real-Time PCR Systems. The primers used in this study were as follows: THRSP forward: 5’-CAGGTGCTAACCAAGCGTTAC-3’, THRSP reverse: 5’-CAGAAGGCTGGGGATCATCA-3’; GAPDH forward: 5’-GGACCTGACCTGCCGTCTAG-3’, GAPDH reverse: 5’-GTAGCCCAGGATGCCCTTGA-3’.

IHC analysis

10 pairs of HCC tissues and paracancerous tissues were fixed in formalin, dehydrated, and embedded in paraffin. The paraffin sections were deparaffinized for antigen retrieval and treated with 3% hydrogen peroxide for blocking peroxidase activity, with 3% bovine serum albumin (BSA) for serum sealing. Afterward, the paraffin sections were incubated with primary THRSP antibody (Guangzhou Alexan Biotech Co., Ltd., China) overnight at 4° C, and then with HRP-conjugated secondary antibody for 50 min at room temperature. 3,3′-Diaminobenzidine (DAB) liquid substrate was used for staining and the hematoxylin solution was used for nucleus counterstaining. Finally, after dehydration and mounting, a microscope was used to acquire images of the staining of tissues.

Transfection of small interfering RNA

Three short interfering RNAs (siRNAs) of THRSP (T1: 5’-ACACCTACTTCACCATGCT-3’; T2: 5’-CCAGGAAATGACGGGACAA-3’; T3: 5’-CATGCACCTCACCGAGAAA-3’) and negative control siRNA (si-NC) were purchased from RiboBio Co., Ltd. (Guangzhou, China). HCCLM3 and Huh-7 cells (2×10⁵ per well) were inoculated on 24-well plates for 24h and then transfected with 25pmol of the RNA duplex according to the manufacturer’s protocol of Lipofectamine 2000 (Invitrogen, Grand Island, NY, USA). After 24h, the transfected cells were harvested for the following experiments.

CCK-8 assay

The transfected HCCLM3 and Huh-7 cells were seeded into 96-well plates (3×10³ cells/well) and incubated for 24h, 48h and 72h. At each time point, 10μl of CCK-8 reagent (Guangzhou Alexan Biotech Co., Ltd., China) was dripped into each well and the cells were cultured for an extra 4h. A micro-plate reader was used to detect the absorbance at 450nm to evaluate cell proliferation.

Transwell and scratch wound healing assays

Transwell assays were used to determine the migration and invasion ability of the HCC cells. The transfected HCCLM3 cells and Huh-7 cells in 200μl serum-free DMEM (5×104 cells/well) were plated into the upper chamber and 600μl of complete medium was added in the lower chamber. For cell invasion assay, the Transwell chamber was coated with Matrigel (Beijing Solarbio Science and Technology Co., Ltd., China). After incubating at 37° C for 24h, the remaining cells in the upper chamber were removed by a cotton swab. The cells that have invaded to the lower surface of the filter were fixed with 4% (v/v) neutral formaldehyde solution (Servicebio, Wuhan, China) for 30min and stained with 0.1% Crystal violet for 30min. Finally, the cells in ten random microscope fields of each filter were counted.

The transfected HCCLM3 and Huh-7 cells (5×10⁵ cells per well) were seeded into 24-well plates. When the cell confluence reaches 90%, a scratch wound was created using a sterilized pipette tip (200μl) on confluent cells. The images of wounds were acquired with a phase-contrast light (40×) at 0h and 24h. The heal area of each scratch wound was determined by ImageJ.

To further understand the molecular mechanism, the potential signaling pathways including the NF-κB and MAPK signaling pathways were detected by Transwell and wound healing assays. Briefly: post of transfection for 48h, the cells were harvested and pre-treated with BAY-11-7082 (5μM, NF-κB inhibitor), AG-126 (10μM, ERK1/2 inhibitor), SB203580 (10μM, p38 MAPK inhibitor) and SP600125 (50μM, JNK inhibitor) for 2h in serum-free DMEM medium. Then, the cells were cultured in the DMEM medium containing inhibitors for 24 hours after being seeded into Transwell chambers or culture wells. The cells treated with dimethyl sulfoxide (DMSO) were used as the negative control.

Flow cytometry

For the cell cycle assays, after 48h of transfection, the HCCLM3 cells were trypsinized, and washed with cold phosphate-buffered saline (PBS) and then fixed in 70% cold ethanol at 4° C for 24 h. After centrifuging and washing, the cells were stained with 500μl PI buffer (50μg/mL, containing RNase, Beijing Leagene Biotechnology Co., Ltd., China) at 37° C in the dark for 30 min. The cell cycle distribution was determined by the Flow cytometer after PI staining. The apoptosis analysis was performed following the instruction of Annexin-V Apoptosis Detection kit (Jiangsu KeyGEN BioTECH Co., Ltd., China). Briefly, 72h post-transfection, the trypsinized cells were washed with cold PBS twice and re-suspended in binding buffer, and then stained with Annexin V-FITC and PI at room temperature for 5-15 min in the dark. The apoptotic cells were analyzed by flow cytometry within an hour.

Western blotting

The total protein was extracted with the nucleoprotein and cytoplasmic protein extraction kit (Jiangsu Keygen Biotech Co., Ltd. China) and quantified by the BCA protein assay kit (Beijing Bomaide Gene Technology Co., Ltd., China). After that, the protein solution was subjected to electrophoresis, detached via SDS-PAGE (12% gels), and transferred to PVDF membranes. The PVDF membranes were blocked with 5% skim milk powder in TBST (Tris buffered saline with 0.5% Tween 20) for 2 h, and then incubated with primary antibodies overnight at 4° C. After washing with TBST for 5 times, the PVDF membranes were incubated with secondary antibody at room temperature for an hour. The protein bands were exposed by enhanced chemiluminescent (ECL) substrate kit (Labgic Technology Co., Ltd. Hefei, China) and analyzed by ImageJ software. The Anti-Bcl-2 antibody, Anti-Caspase-3 antibody, Anti-Bax antibody, Anti-NF-κB p65 antibody and Anti-NF-κB p65 (phospho S536) antibody were purchased from Abcam. The p38 MAPK, Phospho-p38 MAPK (Thr180/Tyr182), P44/42 MAPK (Erk1/2), and Phospho-p44/42 MAPK (Erk1/2) (Thr202/Tyr204) were purchased from Cell Signaling Technology. β-actin was severed as the internal control and anti-β-actin was purchased from Labgic Technology Co., Ltd.

Ethics approval

The manuscript has been approved by the ethics committee of the affiliated hospital of Southwest Medical University (Ethical ID: KY2021142).