Combined identification of ARID1A, CSMD1, and SENP3 as effective prognostic biomarkers for hepatocellular carcinoma

Genetic mutation landscape of HCC in the discovery population

Patient clinical characteristics are presented in Table 1. Briefly, all subjects were males with an average age of 48.33 ± 9.75 years. All patients had a history of HBV infection. The sample collection, sequencing, and data analysis are shown in Figure 1A. To identify the landscape of genetic mutations of HCC patients with HBV infection, WES was performed in 36 pairs of HCC samples; a total of 4231 somatic SNVs, 192 somatic indels, and 12 somatic CNVs were detected. The mean number of SNVs, indels, and CNVs in each patient was 117.5, 5.3, and 0.3, respectively (Figure 1B). According to the frequency of sample mutations, 8 genes (TP53, MUC16, CTNNB1, TTN, ARID1A, PCLO, NBPF10, and CSMD1) were identified with a high mutation frequency (> 15%) in HCCs (Figure 1C). TP53 is the most frequently occurring gene (61.1%). Compared to the most frequently mutated genes in HCC from the TCGA database, NBPF10 and CSMD1 were highly mutated genes in HCCs in the Chinese population. Furthermore, KEGG pathway analysis (P < 0.05, Figure 1D) identified representative HCC proteins and showed that the calcium signaling pathway, axon guidance, circadian entrainment, and nicotine addiction pathways were significantly enriched in HCC samples.

Prognosis-related alterations in the discovery population

To evaluate the prognostic impact of mutated genes HCC patients were divided into three groups: recurrence within 6 months after undergoing LT (group A, n = 10), non-recurrence for > 1 year after undergoing LT (group B, n = 10), and recurrence > 1 year after undergoing LT (group C, n = 8). The top 10 most frequently mutated genes of each group (A, B, and C) are presented in Figure 2A–2C. Compared to patients without tumor recurrence > 1 year after LT (Group B), patients with tumor recurrence within 6 months had significantly different gene landscapes. Only TP53 and CTNNB1 genes were highly and frequently mutated between the two groups. In addition, the types of variations were complicated in patients with early recurrence. These results demonstrated that significant genetic heterogeneity affected HCC prognosis.

To further identify the representative genes associated with HCC prognosis, alterations in the five most frequently mutated genes (TP53, CTNNB1, KMT2C, ARID1A, and PCLO) were detected between groups A and B. CTNNB1, KMT2C, ARID1A, and PCLO were significantly associated with HCC prognosis (likelihood ratio P = 0.018 for CTNNB1 and likelihood ratio P = 0.03 for KMT2C, ARID1A, and PCLO; Table 2). Furthermore, we also compared the expression of these five genes between tumor samples (n = 39) and non-tumor samples (n = 15) based on the expression dataset GSE65372. These findings revealed that KMT2C and ARID1A were significantly and differentially expressed between the two groups (adjusted P-values = 0.046 and 0.039, respectively, Table 2). Therefore, it could be deduced that KMT2C and ARID1A are significantly associated with HCC prognosis.

New HCC driver genes identified in the discovery population

Herein, we attempted to discover new driver genes of HCC in patients with early and late tumor recurrence, using the CHASM software to identify the driver genes in HCC. Finally, except for the known driver genes of HCC, five new driver genes of HCC, including MAP4K3, COX5B, ACTN3, CFTR, and LRRC7, were associated with the early recurrence of HCC (recurrence within 6 months) and 2 genes, PRKCG and SENP3, were considered new driver genes associated with the late recurrence of HCC (recurrence after 1 year) (P < 0.005, false discovery rate (FDR) ≤ 0.1, Table 3). These results showed that the driver genes differed significantly between patients with early and late tumor recurrence.

ARID1A, CSMD1, and SENP3 are associated with poor prognosis in HCC from TCGA database

As described above WES and comparative analysis identified 11 novel genes (KMT2C, ARID1A, NBPF10, CSMD1, MAP4K3, COX5B, ACTN3, CFTR, LRRC7, PRKCG, and SENP3) that were significantly associated with HCC in our Chinese population. To validate these findings, the expression (Supplementary Figure 1) and survival (Supplementary Figure 2) patterns of these genes were observed in the TCGA database. Compared to the control samples, three genes (ARID1A, CAMD1, and SENP3) showed significantly high expression in HCC patients (P < 0.05, Supplementary Figure 1B, 1E, 1J) and were associated with poor prognosis (P < 0.05, Supplementary Figure 2B, 2E, 2J). According to the TCGA database, the data of KMT2C expression were insufficient, while the remaining seven genes did not show any association with HCC prognosis (P > 0.05, Supplementary Figure 2A, 2C, 2D, 2F–2I).

ARID1A, CSMD1, and SENP3 are effective prognostic biomarkers for HBV-infected HCC patients in the replication population

Herein, we focused on ARID1A, CSMD1, and SENP3 during tissue microarray analysis of 180 HCC samples paired with tumor and paracarcinoma tissues as a replication population. These three genes showed significantly high expression in HCC patients in the TCGA database. The clinical characteristics of these patients are presented in Table 4. The IHC analysis of ARID1A, CSMD1, and SENP3 proteins showed that the expression level of all the three genes was significantly different between tumor and paracarcinoma tissues (P < 0.05, Figure 3G–3I). The ARID1A expression was significantly increased in tumor tissues (P < 0.05, Figure 3A, 3B) consistent with the TCGA database. Surprisingly, the expression of CSMD1 and SENP3 was significantly decreased in tumor tissues (P-value < 0.05, Figure 3C, 3D for CSMD1, and Figure 3E, 3F for SENP3).

Furthermore, ROC curve analysis assessed the prognostic power of the three genes based on their expression levels: intensity × positive expression rate. Only SENP3 had significant diagnostic power (AUC = 0.609, P = 0.013, Table 5 and Figure 4), while ARID1A and CSMD1 showed weak diagnostic power (AUC = 0.489 for ARID1A and 0.573 for CSMD1; P > 0.05; Table 5 and Figure 4). Cox regression analysis evaluated the prognostic power of the three genes (Table 6) and showed that SENP3 as a single gene, the combination of two genes (ARID1A and SENP3; CSMD1 and SENP3), and the combination of three genes (ARID1A, CSMD1, and SENP3) could be used as prognostic biomarkers for HCC (P < 0.05, Table 6 and Figure 5). Among these, the combination of CSMD1 and SENP3 genes was the optimal prognostic biomarker for HCC in our replication population (P = 0.006, Table 6 and Figure 5). Moreover, the expression level of SENP3 ≤ 1.2 and CSMD1 < 1.5 was significantly associated with poor HCC prognosis. Based on these results, ARID1A, CSMD1, and SENP3, especially the combination of CSMD1 and SENP3, are effective prognostic biomarkers for HCC individuals in the Chinese population.

CSMD1 prevented HCC by suppressing cell invasion

Finally, to verify the function of ARID1A, SENP3, and CSMD1 genes, we performed cell proliferation and invasion assays in the hepatic Hep3B cell line. In the EdU assay, after the suppression of ARID1A, SENP3, and CSMD1 genes by siRNAs, respectively, for 48 h, the cell proliferation rate was slightly increased by suppression of CSMD1 (Figure 6A, 6D, P > 0.05), slightly reduced by suppression of SENP3 (Figure 6A, 6C, P > 0.05), and almost unaltered by suppression of ARID1A (Figure 6A, 6B, P > 0.05). On the other hand, transwell invasion assay revealed that the cell invasion rate was unaltered by the suppression of SENP3 (Figure 7A, 7C, P > 0.05) but significantly increased by the suppression of CSMD1 (Figure 7A, 7D, P < 0.05), which was consistent with the results of tissue microarray analysis. Also, the suppression of ARID1A elevated the cell invasion rate (Figure 7A, 7B, P < 0.05), which was inconsistent with the results of tissue microarray analysis and TCGA database. Therefore, based on these functional assays, we deduced that CSMD1 prevented HCC by suppressing cell invasion.

Clinical samples

For the discovery population, a total of 36 HCC samples were collected from patients who underwent liver transplant (LT) in our center between 2017 and 2018. All liver grafts were voluntarily donated after death, and informed consent was obtained from all recipients before LT. Tumor tissues obtained from each patient were fixed with 4% paraformaldehyde and embedded in paraffin. Additionally, 10 mL peripheral blood (PB) was collected from each case-matched patient for paired-analysis, as described previously [25, 26]. For the replication cohort, a tissue microarray of 180 HCC samples with HBV infection paired with tumor and paracarcinoma tissues (HLivH180Su08 and HLivH180Su15, Shanghai Outdo Biotech Co. Ltd, Shanghai, China) was performed for immunohistochemical (IHC) analysis. HCC samples were diagnosed by histological analysis.

RNAi, cell culture, and transfection

Hep3B cell line was purchased from ATCC (American Type Culture Collection, Rockville, MD, USA) and cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS) at 37° C in a humidified incubator with 5% CO₂. ARID1A, CSMD1, and SENP3 shRNAs and siRNC were designed and synthesized by Guangzhou RioboBio (Guangzhou, Guangdong, China). The target sequences are listed in Table 7. For transfection, the cells were seeded at a density of 1×10⁵ cells/well for 24 h. Then, the cells were transfected with siRNA fragments and negative control siRNAs.

Whole-exome sequencing (WES) and discovery population analysis

Genomic DNA was extracted from tumor tissues and matched to the PB samples. Genomic DNA libraries were prepared using the protocols provided by Illumina HiSeq2000 platform (Genetron Health Co., Ltd). WES was performed using a TruSeq Exome Enrichment kit (Illumina). The captured DNA libraries were sequenced on an Illumina HiSeq2000 Genome Analyzer. WES provided a 200× average coverage for tumor and PB samples. After raw data (FASTQ) were collected, quality control was performed. The BAM files were processed for local realignment, duplicate removal, and base quality recalibration using GATK (https://software.broadinstitute.org/gatk/) and Picard Tools (http://broadinstitute.github.io/picard/index.html). Mutations were annotated using ANNOVAR [27]. The tumor tissues and normal PB sequencing reads were compared to identify germline single nucleotide variants (SNVs)/insertions-deletions (indels) and somatic SNVs/indels/copy number variants (CNVs)/structural variations (SVs).

Pathway enrichment analysis

Pathway enrichment analyses for genes harboring SNVs, SVs, or differentially expressed genes (DEGs) were performed using Kyoto Encyclopaedia of Genes and Genomes (KEGG) canonical pathways and DAVID Bioinformatics Resources 6.7 database (Annotation, Visualisation, and Integrated Discovery; https://david-d.ncifcrf.gov/). P < 0.05 indicated statistical significance.

GEO dataset retrieval for HCC analysis

The HCC expression dataset GSE65372 [28] was downloaded from GEO (Gene Expression Omnibus) [29] and analyzed using the GEO2R online analyzer. Next, we identified DEGs and compared the cases (tumor samples, n = 39) and controls (non-tumor samples, n = 15) using the GEO2R online analyzer. Adjusted P-values < 0.05 were considered significant, using the Benjamini–Hochberg procedure.

Identification of new driver genes

In a recent study, we used CHASM software (https://wiki.chasmsoftware.org/) to identify the driver genes in HCC, as described previously [30]. Briefly, CHASM predicts the functional significance of somatic missense mutations, using a Random Forest classifier trained with 49 predictive features. This method classifies the predictive data and provides the corresponding scores. Furthermore, the CHASM method is used to test the hypothesis between the scores and the passenger genes in the random forest training set; finally, the P-value is adjusted by Benjamini–Hochberg correction.

TCGA dataset retrieval for HCC analysis

Gene expression and survival data of somatic mutations were downloaded from the TCGA database [31] and analyzed using the UALCAN analyzer (Analyze, Integrate, Discover; http://ualcan.path.uab.edu/index.html). DEGs and survival patterns with P-values < 0.05 between HCC patients and controls in the TCGA database were identified using UALCAN.

IHC analysis for the replication population

The tissue microarray (HLivH180Su08 and HLivH180Su15, Shanghai Outdo Biotech Co. Ltd, China) was incubated in a dry oven at 63° C for 1 h. Subsequent IHC analysis was performed, as described previously [32]. Briefly, the slides were incubated with primary antibodies (ARID1A (ab182560), 1:10000; CSMD1 (ab198906), 1:4000; and SENP3 (ab247139), 1:400, Abcam, UK) overnight at 4° C in a humidified chamber. DAB (3,3’-diaminobenzidine) was used as a chromogenic substrate, and the sections were counterstained with hematoxylin. The IHC staining intensity was classified separately for the nucleus/nuclear membrane and cytoplasm and graded as strong (value = 3), moderate (value = 2), weak (value = 1), or absent (value = 0). In every tissue specimen, three areas with different staining intensities were selected. In each area, 100 cells were observed, and the percentage of positive cells was calculated as X₁%. The other two values were recorded as X₂% and X₃%, respectively. The average of X₁%, X₂%, and X₃% comprised the final staining rate. Scanning and evaluation were performed using Aperio Scanning software (Aperio XT, Leica, Germany).

EdU assay

Hep3B cells (ATCC^® HB-8064) were cultured in a 6-well plate and treated with 100 μL media containing 20 μM EdU. After continuous incubation at 37° C with 5% CO₂ for 24 h and 48 h, the cells were fixed with 4% paraformaldehyde for 15 min and incubated with 0.5% Triton X-100 in phosphate-buffered saline (PBS) for 10–15 min. Fluorescence microscopy was employed to acquire and analyze the images. All the experiments were performed at least three times, and the data were represented as mean ± standard deviation (SD).

Transwell invasion assay

Transwell assays were performed using polyethylene terephthalate-based migration chambers and BD BioCoat Matrigel Invasion Chambers (Becton Dickinson Labware, USA). Hep3B cells were seeded on Matrigel-coated transwell inserts with 200 μL of serum-free medium. The lower chamber was filled with 500 μL medium containing 10% FBS. After incubation for 24 h, the cells remaining on the upper surface of transwell inserts were wiped with cotton wool. The invaded cells were stained with crystal violet for 10 min. Images were captured, and the cell number was counted. All the experiments were performed at least three times, and the data are expressed as mean ± SD.

Statistical analysis

Statistical analysis was performed using SPSS version 17.0 software package (SPSS Inc., Chicago, IL, USA) and GraphPad Prism 7.0 (GraphPad Software, San Diego, CA, USA). Clinical data were presented as mean ± SD. A paired-samples t-test was performed to test the difference between tumor and paracarcinoma tissues for the replication population. Independent samples t-test was performed to test the cell proliferation and invasion rates between siRNA (siARID1A, siSENP3, and siCSMD1) and siRNC groups. Cox regression analysis was used to evaluate the prognostic value adjusted by patient sex, age, and AJCC stage (version VII). We divided the samples into two groups based on the expression level (intensity times positive expression rate) and compared the accumulated survival rate between the two groups for each gene and combination of two or three genes. A receiver operating characteristic (ROC) curve analysis was performed, and the area under the ROC curve (AUC) was calculated to evaluate the prognostic power. The Kaplan–Meier method was used to assess the survival of patients and compared using the log-rank test. P-value < 0.05 indicated statistical significance.

Ethics approval

The clinical HCC samples complied with the Declaration of Helsinki 1975, revised in 2008. This study was approved by the appropriate local institutional review boards on human research at the Huazhong University of Science and Technology (IRB Number: S104). Written consent was obtained from all subjects before participation in the study.

Availability of data and materials

The data have been deposited with links to BioProject accession number PRJNA607376 in the NCBI BioProject database (https://www.ncbi.nlm.nih.gov/bioproject/).