In silico identification of EP400 and TIA1 as critical transcription factors involved in human hepatocellular carcinoma relapse.

Hepatocellular carcinoma (HCC) is the second leading cause of cancer-associated mortality worldwide. Transcription factors (TFs) are crucial proteins that regulate gene expression during cancer progression; however, the roles of TFs in HCC relapse remain unclear. To identify the TFs that drive HCC relapse, the present study constructed co-expression network and identified the Tan module the most relevant to HCC relapse. Numerous hub TFs (highly connected) were subsequently obtained from the Tan module according to the intra-module connectivity and the protein-protein interaction network connectivity. Next, E1A-binding protein p400 (EP400) and TIA1 cytotoxic granule associated RNA binding protein (TIA1) were identified as hub TFs differentially connected between the relapsed and non-relapsed subnetworks. In addition, zinc finger protein 143 (ZNF143) and Yin Yang 1 (YY1) were also identified by using the plugin iRegulon in Cytoscape as master upstream regulatory elements, which could potentially regulate expression of the genes and TFs of the Tan module, respectively. The Kaplan-Meier (KM) curves obtained from KMplot and Gene Expression Profiling Interactive Analysis tools confirmed that the high expression of EP400 and TIA1 were significantly associated with shorter relapse-free survival and disease-free survival of patients with HCC. Furthermore, the KM curves from the UALCAN database demonstrated that high EP400 expression significantly reduced the overall survival of patients with HCC. EP400 and TIA1 may therefore serve as potential prognostic and therapeutic biomarkers. Copyright: © Hong et al.

Introduction

Hepatocellular carcinoma (HCC) is one of the most common types of cancer and the second leading cause of cancer-associated mortality worldwide (1). Its incidence is increasing in numerous countries (2). Progression of HCC is characterized by abnormal cell differentiation, fast infiltrating growth, early metastasis, high-grade malignancy and poor prognosis (3). Liver transplantation (LT) is considered to be one of the major treatment options for HCC (4), as not only it eliminates the tumor but could also cure the underlying liver disease. However, the high relapse rate of HCC following LT, which is estimated between 15 and 20%, remains an important clinical challenge (5). It is therefore crucial to determine the underlying mechanisms of HCC relapse, in order to increase the overall survival of patients with HCC.

Transcription factors (TFs) serve crucial roles in the regulation of tumor progression (6–8). The study of TFs has improved our understanding of the mechanisms underlying the dysregulation of gene expression in cancers. For example, multitudinous compelling evidence have recently showed that HIF-1 plays important roles in many critical aspects of HCC, including tumorigenesis, progression, and metastasis (9–13). Several forkhead box proteins (14) and zinc finger proteins (ZNFs) (15–21) have also been reported to serve crucial roles in HCC. Recently, a study demonstrated that ZNF687 overexpression promotes HCC recurrence (22); however, the TFs associated with HCC relapse remain unknown (22).

Weighted correlation network analysis (WGCNA) has emerged as a powerful technique for multi-gene analysis. This approach is designed to uncover networks and critical genes associated with some phenotypes of interest. WGCNA has been widely used to detect the co-expressed modules (23,24), driver genes (24–27) and driver TFs (28) associated with a disease. In the present study, WGCNA was used to build TFs co-expression network and to investigate critical TFs that may drive HCC relapse. The results from this study may serve at understanding the role of TFs as diagnostic markers of HCC or as therapeutic targets for HCC treatment.

Materials and methods

RNA sequencing (RNA-seq) datasets and preprocessing of data

The RNA-seq data and clinical magnifestations data of the HCC dataset SRP040998 (29) were downloaded from recount2, which is an online RNA-seq resource (30). The dataset consists of tumor samples and matched adjacent normal samples from 21 patients with HCC. Among these patients, nine presented recurrent liver tumors or remote metastasis in the 24 months following orthotopic LT (OLT), whereas the remaining 12 patients were tumor free following OLT (29). Principal component analysis clustered the analyzed samples into the tumor and adjacent groups with only one exception per group. These two samples (sample numbers SRR1220147 and SRR1220148), which were from the same patient, may have been mislabeled and were therefore excluded from this study. Genes with low counts may represent a bias of sequencing. In order to minimize the false positive in differential expression analysis and to speed up WGCNA analysis, only the expressed genes in terms with total counts ≥10 in the samples were kept.

The TF list was obtained from the tftargets (version 1.3) package (https://github.com/slowkow/tftargets), which covers the TFs from Encyclopedia of DNA Elements (31), integrated TF platform (32), RegulatoryCircuits (33), RegulatoryNetworks (34,35), Transcriptional Regulatory Element Database (36) and TRRUST (37) datasets.

Differential expression analysis

Differential gene expression analysis was performed using DESeq2 package (version 1.20.0) (38). A gene was defined as a differentially expressed gene (DEG) when the false discovery rate (FDR) adjusted P-value between the tumor and adjacent groups was ≤0.05 (FDR ≤0.05) and the fold change (FC) is at least 2 times higher or lower (|log2FC| ≥1). However, a transcription factor is defined as a differentially expressed TF (DET) for the cutoff values FDR ≤0.05 and |log2FC| ≥0.6.

WGCNA network construction

A signed WGCNA network between the tumor and matched adjacent tissues was constructed based on the biweight midcorrelation using any gene that was expressed at the total counts value of 10 or higher in ≥90% of samples. Different from an unsigned WGCNA network, which uses the absolute value of the Pearson correlation as an unsigned co-expression similarity measure, the similarity between genes in a signed WGCNA network reflects the sign of the correlation of their expression profiles. Therefore, highly connected hub genes in a signed networks may upregulate adjacent genes since they are positively correlated with them, while in unsigned networks, highly connected hub genes may activate or repress their neighboring genes (39). The count values were normalized by variance stabilizing transformation using DESeq2 package (38). In order to achieve a scale-free topology, soft power parameter was selected based on the criterion of approximate scale-free topology with mininal scale-free fit (SFT) index R2>0.85 and used to derive a pair-wise distance matrix for selected genes using the topological overlap measure. The dynamic hybrid cut method was used to detect clusters of co-expressed genes using R functions in the WGCNA package (version 1.66) (40).

Identification of biological pathways and processes associated with HCC relapse

The biological pathways related to the genes in the module of interest, including the Kyoto Encyclopedia of Genes and Genomes (KEGG) (41), were provided by g:Profiler (42–44).

Analysis of hub TFs in the module of interest

Hub gene is defined as an abbreviation of ‘highly connected gene’, representing a small proportion of nodes with maximal information exchange with other nodes in a network (45). In the present study, hub TFs in the module of interest were defined by module connectivity that was measured by signed eigengene-based connectivity with the cut-off value of absolute kME >0.15. kME, representing the module connectivity of gene k, was determined as the Pearson's correlation coefficient between gene expression values and the module eigengene. In addition, all genes from the module of interest were uploaded to the Search Tool for the Retrieval of Interacting Genes (version 11.0) database (46), by choosing a confidence >0.4 in order to construct a protein-protein interaction (PPI) network. In the PPI network, genes with a connectivity degree >4 (edges) were defined as hub nodes, and TFs with a connectivity degree >4 were defined as hub TFs (47,48).

Differentially connected TFs in the Tan model between relapsed and non-relapsed subnetworks

To identify which TFs were differently regulated in the networks, the differences in connectivity (DiffK) were compared between the relapsed vs. non-relapsed subnetworks according to the following formula: DiffK(i)=K1(i)-K2(i), where K1(i) and K2(i) indicate the connectivity of the gene (i) in the relapsed subnetwork and in the non-relapsed subnetwork, respectively. The connectivity of K1 and K2 were calculated using R functions in the WGCNA package (version 1.66). To facilitate the comparison between the connectivity measures of each network, standardization was carried out in each network by dividing each TF's connectivity by the maximum connectivity of TFs-genes co-expression sub-network according to the following formula: K1(i)=K1(i)/max(K1).

The difference between the connectivity values of two subnetworks was defined as DiffK. DiffK values ranged from −1 to 1. A DiffK value >0 suggested that the TFs were more highly connected in the relapsed subnetwork compared with the non-relapsed subnetwork, whereas a DiffK value <0 indicated that the TFs were more highly connected in the non-relapsed subnetwork compared with the relapsed subnetwork. TFs were defined as differentially connected when the absolute value of DiffK was >0.4 (28).

Upstream TFs regulating gene expression in the module of interest

TFs regulating gene expression in the module of interest were analyzed using the plugin iRegulon (v1.3) (49) in Cytoscape network (version 3.7.1) (50). iRegulon is a computational tool that can identify the upstream TFs and predict direct target genes in a set of human, mouse and Drosophila genes. iRegulon uses >9,000 known position weight matrices from various sources and different species and link them to candidate binding TFs using a ‘motif2TF’ procedure. This allows to link motifs of TFs from other species to candidate human TFs. Predicted upstream TFs are rated and grouped according to the Normalized Enrichment Score (NES) (28,51).

Kaplan-Meier (KM) survival analysis

To evaluate the prognostic potential of hub TFs, UALCAN (52), Gene Expression Profiling Interactive Analysis (GEPIA) (53) and KMplot (54) tools were used to obtain the KM curves for the overall survival (OS), disease-free survival (DFS) and relapse-free survival (RFS) of patients with HCC, respectively, for the genes of interest. These online tools provides the data on the effect of genes on cancer survival, including HCC survival.

Correlation analysis

Correlation analysis for binary and continuous variables was performed using standard screening binary trait and standard screening numeric trait functions in the WGCNA package (40), respectively. The output of these functions include q-values of the correlations calculated from the P-values using an optimised false discovery rate approach.

Results

Identification of modules associated with HCC recurrence using signed WGCNA

Once the two mislabeled samples and lowly-expressed genes were excluded, the remaining 26,324 genes from 58,037 genes detected in HCC samples were grouped into 32 modules according to their expression profiles via hierarchical clustering (Fig. 1). To guarantee a scale-free topology (minimal SFT R2>0.85), the soft power parameter β-value was set at 12 (Fig. 1A). The pairwise correlation was converted into adjacency matrix of connection strengths through soft-thresholding approach, in order to construct a dissimilarity matrix based on topological overlap measure (TOM) and identify gene modules through a dynamic tree-cutting algorithm. Each module was assigned to a corresponding color (Fig. 1B). The module eigengenes were calculated by the first principal component to represent each module. With module eigengenes, Tan [Spearman correlation, r=−0.58; P-value adjusted, (padj)=1×10−4] and Steelblue (r=0.57; padj=2×10−4) modules were found to be highly associated with HCC recurrence. The Steelblue module contained only two TFs. Although these two TFs, poly(rC) binding protein 2 (55) and eukaryotic translation initiation factor 4B (56), have been reported to participate in HCC development, the present study did not detect differences in their expression between the tumor and adjacent normal tissues. Subsequently, the present study focused on the Tan module, which eigengene levels were also associated with primary tumor grade (Spearman correlation, r=−0.61; padj=4×10−5) and α-fetoprotein (AFP) (Pearson correlation, r=−0.41; padj=9×10−3; Fig. 2).

Figure 1.

Gene co-expression network identified by weighted gene co-expression network analysis. (A) Network topology for different soft-thresholding powers. Numbers in the plots indicate the corresponding soft thresholding powers. The approximate scale-free topology can be attained at the soft-thresholding power of 12. (B) Gene dendrogram obtained by clustering the dissimilarity based on topological overlap with the corresponding module colors indicated by the color row. Each colored row represents a color-coded module which contains a group of co-expressed genes. A total of 32 modules were identified.

Figure 2.

Relationships between module eigengenes and variables of hepatocellular carcinoma. Each row in the table corresponds to a module and each column corresponds to a clinical trait. The module name is shown on the left side of each cell. Numbers in the table corresponds to the correlations of the corresponding module eigengenes and clinical traits, with the P-values printed below the correlations. Intensity and direction of correlations are indicated on the right side of the heatmap (red, positively correlated; green, negatively correlated).

The Tan module consisted of 615 genes, including 55 transcription factors. The results form KEGG pathway analysis demonstrated that the Fanconi Anemia pathway was enriched in the Tan module (adjusted P-values, 1.782×10−2).

Hub TFs in the module of interest

Hub genes were first screened from the Tan module based on the eigengene-based connectivity. By using the cut-off value of absolute kME>0.15, 17 TFs with high connectivity were identified in the Tan module. Among them, T-cell-restricted intracellular antigen-1 (TIA1) had the highest connectivity. Intraflagellar transport protein 80 homolog (IFT80), mediator of DNA damage checkpoint 1 (MDC1) and zinc finger protein 260 (ZNF260) were differentially expressed between the tumor and adjacent tissues (Table I).

Table I.

Expression differences of the hub transcription factors between the tumor and adjacent tissues.

	Tumor vs. adjacent
Ensembl ID	Log2 fold change	P-value	Adjusted P-value	kME	Gene symbol
ENSG00000116001	0.28	0.04	0.08	0.314	TIA1
ENSG00000179912	0.33	0.02	0.04	−0.281	R3HDM2
ENSG00000163666	0.29	0.15	0.24	0.264	HESX1
ENSG00000124459	0.46	0.01	0.02	−0.260	ZNF45
ENSG00000185670	0.28	0.04	0.07	−0.256	ZBTB3
ENSG00000148200	0.14	0.40	0.51	0.253	NR6A1
ENSG00000167081	0.39	0.04	0.09	0.239	PBX3
ENSG00000162086	0.12	0.43	0.54	−0.200	ZNF75A
ENSG00000089335	0.44	0.01	0.02	0.197	ZNF302
ENSG00000177463	0.33	0.04	0.08	0.192	NR2C2
ENSG00000133111	0.55	<0.01	<0.01	0.189	RFXAP
ENSG00000100307	0.34	0.07	0.13	0.186	CBX7
ENSG00000254004	0.62	<0.01	<0.01	−0.172	ZNF260
ENSG00000173153	0.16	0.21	0.32	0.169	ESRRA
ENSG00000137337	0.65	<0.01	<0.01	−0.166	MDC1
ENSG00000068885	0.75	<0.01	<0.01	0.155	IFT80
ENSG00000166432	0.08	0.76	0.83	0.151	ZMAT1

Differential expressions were analyzed using R package DESeq2 (version 1.20.0). kME indicates the eigengene-based connectivity of gene co-expression network.

Hub genes were also screened from the Tan module according to the connectivity of genes in the PPI network (Table SI). By using the cutoff of confidence >0.4 and the connectivity degree of 4 (node/edge), nine TFs, including EP400, were identified as the hub TFs (Table SI).

These two approaches identified 23 hub TFs in total, including TIA1, IFT80 and MDC1 that were identified by both approaches.

Analysis of the differentially connected TFs between the relapsed and non-relapsed subnetworks in the module of interest

Genes that are differentially co-expressed between different sample groups are more likely to be regulators, and may therefore explain differences between phenotypes. Differential network analysis was performed between the relapsed and non-relapsed subnetworks in the Tan module. The analysis identified TIA1 and nuclear receptor subfamily 2 group C member 2 (NR2C2) as differentially connected TFs which were more highly connected in the relapsed subnetwork compared with the non-relapsed subnetwork. Furuthermore, EP400 and five other TFs were differentially connected TFs that were more highly connected in the non-relapsed subnetwork compared with the relapsed subnetwork (Table II).

Table II.

Differentially connected transcription factors between the relapsed and non-relapsed subnetworks and their expression differences between the tumor and adjacent normal tissues.

Ensembl ID	DiffK	Log2 fold change	P-value	Adjusted P-value	Gene symbol
ENSG00000117569	−0.58	0.60	<0.01	<0.01	PTBP2
ENSG00000267680	−0.55	0.47	0.01	0.02	ZNF224
ENSG00000076108	−0.54	0.39	<0.01	0.01	BAZ2A
ENSG00000183495	−0.47	0.41	<0.01	0.01	EP400
ENSG00000185946	−0.45	0.49	0.01	0.02	RNPC3
ENSG00000186908	−0.40	0.13	0.24	0.35	ZDHHC17
ENSG00000116001	0.40	0.28	0.04	0.08	TIA1
ENSG00000177463	0.45	0.33	0.04	0.08	NR2C2

Differential expressions were analyzed using R package DESeq2 (version 1.20.0). DiffK indicates the differential connectivity between the relapsed and non-relapsed subnetwork.

Master regulatory factors for the module of interest

It is regarded that genes co-expressed or participating in the same biological process may be regulated by the same or similar TFs (51). In order to gain insight into the upstream regulators of gene expression in the module of interest, iRegulon was used to search the user-defined space for motifs enriched around the transcription start site of the genes in the Tan module. The most enriched TF motif was ZNF143 with NES 6.39. There were 42 target genes of ZNF143 in the Tan module. Among them, TRIM74 was negatively correlated with ZNF143 expression level, and 15 were positively correlated with ZNF143 expression level, including seven TFs such as MDC1 and ZNF260 (Table III).

Table III.

Correlation of zinc finger protein 143 with its target genes in the Tan module.

Ensembl ID	Correlation coefficient	P-value	q-value	TFs	Gene symbol
ENSG00000089335	0.43	0.005	0.014	Yes	ZNF302
ENSG00000109118	0.37	0.019	0.033	Yes	PHF12
ENSG00000124459	0.40	0.010	0.022	Yes	ZNF45
ENSG00000125945	0.39	0.012	0.024	No	ZNF436
ENSG00000127483	0.51	0.001	0.004	No	HP1BP3
ENSG00000127957	0.38	0.015	0.029	No	PMS2P3
ENSG00000131115	0.50	0.001	0.006	Yes	ZNF227
ENSG00000137337	0.41	0.009	0.021	Yes	MDC1
ENSG00000155428	−0.39	0.014	0.027	No	TRIM74
ENSG00000159905	0.38	0.015	0.028	No	ZNF221
ENSG00000175787	0.42	0.007	0.017	Yes	ZNF169
ENSG00000196597	0.33	0.035	0.050	No	ZNF782
ENSG00000198707	0.47	0.002	0.009	No	CEP290
ENSG00000254004	0.62	<0.001	0.001	Yes	ZNF260
ENSG00000256294	0.40	0.010	0.022	No	ZNF225
ENSG00000263002	0.50	0.001	0.006	No	ZNF234

Pearson correction coefficient was calculated using R package (version 3.5.1). q-values are adjusted P-values, which were calculated from the P-values using an optimized false discovery rate approach using R function in the weighted correlation network analysis package. TFs, transcription factors.

iRegulon predicted 30 TFs for ZNF143. Among them, 4 were negatively correlated with ZNF143 in the gene expression levels (Table IV), and YY1 together with another 6 TFs were positively correlated with ZNF143 expression level.

Table IV.

Correlation of zinc finger protein 143 with its upstream transcription factors.

Ensembl ID	Correlation coefficient	P-value	q-value	Gene symbol
ENSG00000078399	0.45	0.003	0.011	HOXA9
ENSG00000100811	0.44	0.004	0.013	YY1
ENSG00000105866	0.49	0.001	0.007	SP4
ENSG00000106689	−0.54	<0.001	0.003	LHX2
ENSG00000147133	0.52	0.001	0.004	TAF1
ENSG00000164002	0.37	0.019	0.033	EXO5
ENSG00000167182	0.52	0.001	0.004	SP2
ENSG00000168066	−0.40	0.010	0.022	SF1
ENSG00000172845	0.40	0.011	0.024	SP3
ENSG00000177374	−0.47	0.002	0.009	HIC1
ENSG00000185630	−0.45	0.004	0.012	PBX1

Pearson correction coefficient was calculated using R package (version 3.5.1). q-values are adjusted P-values, which were calculated from the P-values using an optimised false discovery rate approach using R function in the weighted correlation network analysis package.

The most enriched TF motifs were also analyzed in the dataset consisting of 55 TFs in the Tan module, where the most enriched TF motifs were YY1 with NES 5.54. YY1 had 36 target TFs. Among these target TFs, 12 were positively correlated with YY1 expression level (Table V), including MDC1 and ZNF260.

Table V.

Correlation of Yin Yang 1 with its target transcription factors.

Ensembl ID	Correlation coefficient	P-value	q-value	Gene symbol
ENSG00000005339	0.46	0.003	0.012	CREBBP
ENSG00000076108	0.36	0.023	0.048	BAZ2A
ENSG00000089335	0.39	0.013	0.034	ZNF302
ENSG00000106261	0.37	0.018	0.042	ZKSCAN1
ENSG00000117569	0.37	0.019	0.043	PTBP2
ENSG00000124177	0.53	<0.001	0.004	CHD6
ENSG00000131051	0.46	0.003	0.012	RBM39
ENSG00000131115	0.38	0.016	0.038	ZNF227
ENSG00000137337	0.52	0.001	0.004	MDC1
ENSG00000173575	0.37	0.019	0.043	CHD2
ENSG00000254004	0.46	0.003	0.013	ZNF260
ENSG00000166478	0.44	0.004	0.016	ZNF143

Pearson correction coefficient was calculated using R package (version 3.5.1). q-values are adjusted P-values, which were calculated from the P-values using an optimized false discovery rate approach using R function in the weighted correlation network analysis package.

iRegulon predicted 20 TFs for YY1. Among them, one was negatively correlated with YY expression level, and six were positively correlated with YY1 expression level (Table VI), including YY1 itself.

Table VI.

Correlation of Yin Yang 1 with its upstream transcription factors.

Ensembl ID	Correlation coefficient	P-value	q-value	Gene symbol
ENSG00000100811	1.00	<0.001	<0.001	YY1
ENSG00000105866	0.38	0.016	0.039	SP4
ENSG00000112592	0.39	0.014	0.035	TBP
ENSG00000147133	0.42	0.008	0.024	TAF1
ENSG00000167182	0.48	0.002	0.009	SP2
ENSG00000172845	0.65	<0.001	<0.001	SP3
ENSG00000177374	−0.58	<0.001	0.001	HIC1

Pearson correction coefficient was calculated using R package (version 3.5.1). q-values are adjusted P-values, which were calculated from the P-values using an optimized false discovery rate approach using R function in the weighted correlation network analysis package.

KM survival analysis

To validate the effects of the 23 hub TFs and two master regulatory factors associated with HCC relapse, KMplot and GEPIA were used to retrieve the RFS curves of 313 patients with HCC and the DFS curve of 362 patients with HCC, respectively (Table SII). In order to determine the prognostic relevance of these TFs in HCC, the OS of patients with HCC was determined for each TF according to its expression using UALCAN (Table SII), which contains data from 365 patients with HCC.

The results from KMplot and GEPIA demonstrated that TIA1, ZNF260 and EP400 expression levels were associated with the RFS and DFS of patients with HCC (Fig. 3A and B). Furthermore, ZNF260 and EP400 were also associated with the OS of patients with HCC, however, TIA1 was only trend associated with the OS of patients with HCC (P=0.09) according to UALCAN (Table SII). In addition, the results demonstrated that MDC1 was associated with the RFS and DFS of patients with HCC according to KMplot and UALCAN analyses. However, according to GEPIA, MDC1 was not associated with the DFS of patients with HCC (P=0.07) (Table SII).

Figure 3.

KM survival curves for TIA1, NR6A1, RFXAP, ZNF260, EP400 and PTBP2. Red color indicates high expression and grey/blue colors indicate low expression. (A) Effect of TIA1, NR6A1, RFXAP, ZNF260, EP400 and PTBP2 expression levels on RFS of patients with HCC (data from KMplot). (B) Effect of TIA1, NR6A1, RFXAP, ZNF260, EP400 and PTBP2 expression levels on DFS of patients with HCC (data from Gene Expression Profiling Interactive Analysis). DFS, disease-free survival; EP400, E1A binding protein p400; NR6A1, nuclear receptor subfamily 6 group A member 1; PTBP2, polypyrimidine tract binding protein 2; RFS, relapse-free survival; RFXAP, regulatory factor X associated protein; TIA1, TIA1 cytotoxic granule associated RNA binding protein; ZNF260, zinc finger protein 260; KM, Kaplan-Meier.

Two master regulatory factors, ZNF143 and YY1, were demonstrated to significantly influence the OS of patients with HCC, but not the RFS.

Discussion

Although liver resection is the most effective curative treatment for HCC, relapse remains frequent (57,58). Investigating the underlying mechanisms of HCC recurrence may therefore lead to the development of novel therapeutic strategies and prognostic biomarkers. Recent studies reported that miR-125b (59) and miR-1246 (60) could be considered as novel biomarkers for HCC relapse. miR-125b has been described as a tumor suppressor that induces cellular senescence and apoptosis in hepatocellular carcinogenesis by targeting sirtuin6 (61). Similarly, miR-1246 has been reported to promote the cell apoptosis of HCC cells (62), and enhance migration and invasion in HCC (63).

In addition to miRNAs, TFs may also serve important roles in HCC relapse. In order to identify the TFs that cooperatively drive HCC relapse, the present study constructed signed WGCNA gene co-expression network by using the tumor and matched adjacent normal samples of 20 patients with HCC. The Tan module was found to be associated with HCC relapse and HCC staging, and was enriched in the Fanconi Anemia pathway. Furthermore, 23 hub TFs were identified, including TIA1, IFT80 and MDC1, which were identified by both eigengene-based and PPI network approaches. Eight TFs, including TIA1 and EP400, were found to be differentially connected when comparing the non-relapsed and relapsed subnetworks. ZNF143 and YY1 were detected as master upstream regulator genes that could potentially regulate the expression of the genes and TFs of the Tan module, respectively. KM survival analysis demonstrated that TIA1, ZNF260 and EP400 expression levels were associated with the RFS and DFS of patients with HCC.

Among the 23 hub TFs, TIA1, EP400 and NR2C2 were differentially connected when comparing the non-relapsed and relapsed subnetworks. These TFs were less cooperative with other genes in the Tan module of patients with HCC relapse patients compared with patients with non-relapse, suggesting that they may serve important role in HCC relapse.

EP400 is a SWItch/sucrose non-fermentable DNA-dependent ATPase (64) that can alter chromatin structure during DNA double-strand break repair. EP400 has also been reported to be implicated in HCC (65). Similarly, the present study demonstrated that EP400 was associated with the RFS, DFS and OS of patients with HCC. Previous studies reported that preoperative AFP serum level (66) and half live of serum AFP (67), which are used to calculate the rate of serum AFP decline, could be considered as an early prognostic index of HCC relapse and patients survival following hepatic resection. The module eigengene provides the most appropriate synopsis of gene expression profiles of any given module. In the current study the eigengene levels of the Tan module were associated with HCC relapse and with AFP expression levels. EP400 was a critical TF in the Tan module with a PPI connectivity of 9 (Table SI), and its expression level was weakly correlated with AFP expression level. These results suggested that EP400 may serve a crucial role in HCC relapse.

TIA1 is a member of an RNA-binding protein family and possesses nucleolytic activity against cytotoxic lymphocyte target cells (68). In the present study, TIA1 was found to be associated with the RFS and DFS of patients with HCC. However, TIA1 was not associated with the OS of patients with HCC. A recent study reported that high TIA-1 expression level is associated with the poor survival rate of patients with HCC (69). Furthermore, TIA-1 can regulate IGF binding protein-3 (IGFBP3) at the post-transcriptional level in human HCC cells (70). IGFBP3 is the primary binding protein of IGF-I, and IGF-I has been reported to be involved in early HCC relapse (71). Therefore, TIA1-mediated IGFBP3 regulation may serve an important role in HCC relapse.

In addition to TIA1, MDC1 was found to be a hub TF in the present study. MDC1 was also differentially expressed between the tumor and adjacent normal tissues. MDC1 is involved in checkpoint activation and subsequent DNA repair following DNA damage (72,73). Previous studies reported MDC1 contributes to breast cancer (74) and pancreatic cancer (75); however, its role in HCC remains unknown. The result from the present study demonstrated that MDC1 reduced the OS and RFS of patients with HCC.

This study demonstrated that ZNF143 was a master upstream regulator gene that could potentially regulate expression of the genes in the Tan module. In addition, ZNF143 expression level was correlated with 16 of its target genes, including MDC1 and ZNF260. ZNF143 is a chromatin-looping factor that contributes to the architectural foundation of the genome by providing sequence specificity at promoters connected with distal regulatory elements (76). A previous study reported that ZNF143 expression is activated following DNA damage induced by etoposide, cisplatin and Adriamycin (77). Furthermore, ZNF143 is involved in cellular motility via the zinc finger E-box binding homeobox 1-cadherin-linked pathway in colon cancer cells (78). ZNF143 protein level is also correlated with clinical outcomes in patients with lung adenocarcinoma (79). A previous study demonstrated that ZNF143 activity inhibition by small molecules induced tumor regression in vitro and in vivo (80). Similarly, the present study demonstrated that ZNF143 expression level was associated with the OS of patients with HCC; however, ZNF143 had no influence on the RFS and DFS of patients with HCC (Table SII).

In the present study ZNF143, MDC1 and ZNF260, and YY1 itself, were all predicted targets of YY1. Furthermore, they were all positively correlated with YY1 expression level (Table V). YY1 belongs to the polycomb group of proteins; this type of protein may cause epigenetic remodeling of the chromatin and therefore dynamically regulate expressions of their target genes (81). YY1 overexpression observed in the majority of cancers has been correlated with poor prognosis of patients (82). Previous studies have demonstrated that YY1 is implicated in HCC (83–88). It has been reported that YY1 acts predominantly as an epigenetic modulator, influencing the activity and/or localization of epigenetic modifiers molecules, including DNA methylation transferases, histone deacetylases or non-coding RNAs (81). YY1 may therefore increase expressions of ZNF143 and MDC1, and dampen DNA repair pathways in HCC progression. Similarly, the present study demonstrated that YY1 expression level was associated with the OS of patients with HCC; however, YY1 was not associated with the RFS of patients with HCC.

The Fanconi anemia pathway is essential for the repair of DNA damage and is involved in three classic DNA repair pathways named homologous recombination, nucleotide excision repair and mutagenic translesion synthesis (89). A recent study demonstrated that genes from the Fanconi anemia/BReastCAncer pathway are involved in HCC chemoresistance (90). In the present study, the Tan module was enriched in the Fanconi Anemia pathway, and several hub TFs involved in DNA repair, including MDC1 and EP400, were associated with the RFS of patients with HCC, suggesting that dysfunction in DNA repair pathways may be important mechanism involved in HCC relapse.

The DNA repair pathways may play a role in HCC relapse. The hub TFs TIA1 and EP400 were differentially connected between the non-relapsed and relapsed subnetworks. TIA1 and EP400 may be considered as critical drivers for HCC relapse and serve therefore as promising targets of HCC relapse. However, further investigation is required to confirm these in silico results.