This article is mentioned in:

Introduction

Circadian rhythms control the 24 h cycle of specific metabolic functions required by living beings, ensuring an efficient metabolic homeostasis (1,2). In humans, the circadian rhythms are controlled by a master pacemaker situated in the suprachiasmatic nuclei of the hypothalamus, which is synchronized to the photoperiod (1,2). The molecular clock of transcription involves a translational feedback mechanism of genes, including clock circadian regulator (CLOCK), period circadian clock (PER), aryl hydrocarbon receptor nuclear translocator like (BMAL1) and cryptochrome circadian clock (CRY), able to regulate a number of physiological properties, including body temperature, melatonin secretion, hormone secretion, blood pressure and the sleep-wake cycle (3).

It is well-known that disturbances in the circadian rhythm may cause the development of diseases, including major depressive disorder, seasonal affective disorder, schizophrenia, bipolar disorder (4–9), stress, desynchronosis (9), anxiety disorder, diabetes (10), obesity, diseases associated with aging (11), genome instability (12) and cancer (13,14). Prior studies have demonstrated the association between circadian rhythm alterations and the development of breast (15) and prostate cancer (16), B-cell lymphoma (17), non-small cell lung (18), testicular (19) and ovarian cancer (20).

The liver has a central and unique metabolic function in maintaining energy homeostasis via glycolysis and gluconeogenesis associated with fatty acid metabolism (biosynthesis/beta oxidation) (21). Rhythmic fluctuations have been identified in hepatic metabolic functions with a 24-h periodicity (22). Previous studies have demonstrated that liver cancer initiation may be due to alterations in circadian rhythmic genes, including PER3 (23) and CRY genes, and casein kinases (24). Additionally, it has been revealed that the dysregulation of metallothionein-1 (MT-1), MT-2 and metal transcription factor-1 are involved in the alterations to circadian rhythms present in liver cancer (25). Furthermore, the presence of hepatitis viral infection, which is a cause of liver cancer, has been revealed to cause dysregulation of the expression of circadian genes (26). However, it is crucial to underline that liver cancer is additionally caused by chronic exposure to toxic chemicals, hepatosteatosis, type 2 diabetes and obesity (27–29). The underlying molecular mechanism of circadian clock disruption in non-viral liver cancer remains unknown.

Our previous study analyzed the transcriptome of a human line of hepatoblastoma cells without viral infection (HepG2), compared with normal hepatocytes, and the gene expression data, obtained from the liver tissues of patients with hepatitis C virus (HCV), HCV-associated cirrhosis and liver cancer with HCV-associated cirrhosis, using the publicly available E-MTAB-950 database (30). Despite the HepG2 cell line being revealed as a misidentified hepatocellular carcinoma cell line and subsequently identified as hepatoblastoma by the International Cell Line Authentication Committee (iclac.org/databases/cross-contaminations), the outcomes of our studies are validated because they focused on a liver cancer cell line without viral infection. Our previous analysis enabled the identification of specific clusters of genes for the various stages of liver cancer, and allowed the isolation of a network of 26 HUB genes that specifically control critical metabolic functions, independent of the viral infection (30). All 26 HUB genes were revealed to encode intrinsically disordered proteins (IDPs), thus they exhibited multifunctional behaviors and were involved in metabolic cellular control.

The present study aimed to identify whether these HUB genes (and associated proteins) were components of the network that controls the human circadian rhythm. Therefore, the network of genes involved in circadian rhythms was extracted from the human interactome to identify the nodes with high centrality and interactions with the most studied circadian gene, CLOCK. Furthermore, the association between the circadian network and 26 HUB genes in the HepG2 network was evaluated. In this way, we identified the genes linking the two main networks and evaluated if the proteins coded from them have high disorder propensities. Finally, comparing the deregulated microRNA (miRNA/miR) in HepG2 cells with miRNAs in normal hepatocytes, the target genes were predicted in order to evaluate if there were miRNAs involved in the network between HepG2 cells and circadian networks.

Materials and methods

Network analysis

A network of 71 genes was extracted, from the human interactome compiled from various databases, including Pathway Commons (31), Biological General Repository for interaction Datasets (BioGRID) (32), Human Protein Reference Database (33), ConsensusPathDB (34), Database of Interacting Proteins (35) and the Breast Cancer Information Core and Michigan Molecular Interactions (36), and identified to be involved in circadian rhythms. Only the connected component of these 71 seed networks were considered for statistical analysis using different tools, including Network-Analyzer (37), the Database for Annotation, Visualization and Integrated Discovery (DAVID) (38) and the Biological Networks Gene Ontology tool (39). Using the Cytoscape 3.5 package (www.cytoscape.org), statistical analysis was performed to evaluate the following three measures of centrality: i) The degree, which indicates the number of interactions of a particular node with other nodes in the network; ii) the betweenness, which evaluates the importance of a node in the network and how the other interactions in the network are controlled by this node (40); iii) the closeness centrality of a node, which measures the speed of information flow through this node to reachable nodes in the network and ranges from 0–1 (41). However, the power law is used to predict the HUB nodes that have functions in the network. The power law details the functional association between two quantities, where one quantity varies as a power of another; based on the power law distribution degree, a network may be defined as scale-free indicating that ‘riches get richer’ (42–44).

Other topological analysis, including average characteristic path length, network density, centralization and heterogeneity were evaluated using the Cytoscape 3.5. package (41,45). The characteristic path length is calculated by identifying the shortest path between all pairs of nodes, adding them and dividing by the total number of pairs. This indicates the number of steps required to get from one member of the network to another (41,45). The density of a network is defined as a ratio of the number of edges to the number of possible edges (46); whereas, the centralization produces rankings identify the most important nodes in a network model (47). In particular, networks with topologies resembling a star have a centralization close to 1, whereas decentralized networks are characterized by having a centralization of ~0.47. Furthermore, the network heterogeneity, evaluated using the Cytoscape 3.5. package, reflects the tendency of a network to contain HUB nodes (48). Finally, a cluster analysis was performed, which groups similar objects to form clusters; therefore, objects in the same cluster are more similar to each other, compared with those in other clusters. In particular, the overlapping clusters were calculated on the basis of cohesiveness quality functions (49).

Disorder propensity analysis

The associated protein sequences corresponding to the HUB genes common between circadian and liver cancer networks were extracted from the UniProt database (www.uniprot.org). To assess the proportion of residues involved in intrinsic disorder, the DisProt tool (50) was used to subdivide the sequences into three major groups extracted on the basis of similar contents of disorder (10–15%, 16–50% and >50%).

Target genes prediction of miRNAs

A list of miRNAs (51) that have been identified as dysregulated in HepG2 cells, compared with human normal hepatocytes were selected. Furthermore, predictions of miRNA complementarity to 3′ untranslated regions (UTRs) in mRNAs were performed using three commonly used tools for target prediction: TargetScan Human 6.2 (www.targetscan.org) (52), PITA (53) and miRanda (www.microrna.org) (54). This analysis was based on identifying conserved sites that match the seed region of each miRNA, corresponding to the position between nucleotides 2 and 8in mature miRNAs. A list of putative targets for each miRNA was obtained and those predicted from 2/3 tools were selected for functional annotation analysis of pathways, which was performed using the DAVID program and by selecting the more significantly enriched pathways, with a number of genes >60 and P<0.05 (38).

Results and discussion

Human circadian network. The network of genes involved in circadian rhythms, on the basis of seed nodes, was extracted from the human interactome, and included the following genes: Aryl Hydrocarbon Receptor Nuclear Translocator Like (ARNTL), casein kinase (CSNK)1E, inter-α-trypsin inhibitor heavy chain family member 5 (ITIH5), replication factor C subunit 3 (RFC3), WD repeat domain 41, PER1, CSNK1D, syntrophinβ2 (SNTB2), acyl-CoA thioesterase 13, chondroitin sulfate N-acetylgalactosaminyltransferase 1, PER2, ARNTL2, PDZ domain containing ring finger 3, growth arrest specific 2, fibronectin leucine rich transmembrane protein 1, PER3, neuronal PAS domain protein 2, low density lipoprotein receptor, zinc finger protein (ZNF) 286A, G protein-coupled receptor (GPR)116, nuclear receptor subfamily 1 group D member 1, cyclin dependent kinase (CDK) L5, splicing factor proline and glutamine rich (SFPQ), adipogenesis regulatory factor, translocase of inner mitochondrial membrane 8A, basic helix-loop-helix family member E (BHLHE)40, 7-dehydrocholesterol reductase, solute carrier family 39 member 14, suppressor of cytokine signaling 2 (SOCS2), GPR6, BHLHE41, histone deacetylase 4 (HDAC4), HLF PAR bZIP transcription factor (HLF), solute carrier organic anion transporter family member 4A1, γ-secretase activating protein, BMAL1, methyl-CpG binding protein 2 (MECP2), ETS variant 5 (ETV5), Kruppel like factor 11 (KLF11), ZNF394, D-box binding PAR bZIP transcription factor, neurexin 1, TNFAIP3 interacting protein 2 (TNIP2), exocyst complex component 1 (EXOC1), extended synaptotagmin 1 (ESYT1), nuclear receptor subfamily 1 group D member 2, SH3 and multiple Ankyrin repeat domains 3 (SHANK3), zw10 kinetochore protein (ZW10), phospholipid scramblase 1 (PLSCR1), CLOCK, solute carrier family 2 member 1 (SLC2A1), hydrocretin receptor 2, 5-methyltetrahydrofolate-homocysteine methyltransferase (MTR), transferrin receptor (TFRC), casein kinase 1-α-1, synemin, sprouty RTK signaling antagonist 4 (SPRY4), ubiquitin specific peptidase 2, glycogen synthase kinase 3β (GSK3B), hypoxia inducible lipid droplet associated, Scmpolycomb group protein like 1, CRY1, nuclear factor interleukin 3 regulated, ATPase H+/K+ transporting α subunit, Ras homolog family member B, CRY2, insulin induced gene 1, unc-13 homolog A (UNC13A) and apolipoprotein L domain containing 1 (APOLD1). The human circadian network consists of 2151 nodes and 75821 interactions (Table I). The circadian network was identified to be highly centralized (0.235); a higher value of centralization indicates that the network is concentrated in the center with an overall integration towards the high degree nodes. The network density of the circadian network, which describes the proportion of potential connections in a network that are actual connections, as a measure of network effectiveness, was identified to be 0.033. In addition, the circadian network exhibited a high value of heterogeneity, which demonstrates its tendency to contain HUB nodes. The characteristic path length was identified to be 2.373, whereas the average number of neighbors was 70.5.

Table I. Statistical analysis of the genes involved in the network obtained for human circadian rhythms.

Table I.

Statistical analysis of the genes involved in the network obtained for human circadian rhythms.

Statistical analysis	Circadian network
Number of nodes	2151
Number of interactions	75821
Network density	0.033
Network centralization	0.235
Characteristic path length	2.373
Network heterogeneity	1.012
Neighbor average number	70.498

The human circadian network was demonstrated to follow the small-world rule (41), as the characteristic path length is very short. The nodes that exhibited a high centrality were small ubiquitin-like modifier 2, CDK2, heat shock protein 90A, p53, nuclear respiratory factor 1 and GSK3B. As CLOCK is one of the most studied circadian genes, its sub-network was extracted from the general network of circadian genes using the Cytoscape tool, which demonstrated that it contained 87 nodes with 86 direct interactions. The analysis of the present study demonstrated that CLOCK is associated with other genes including proliferating cell nuclear antigen, PER and sirtuin 1, 3 and 5.

Association between the circadian network and 26 HUBs in the HepG2 network

Four networks have previously been compared, each obtained from the differentially expressed genes in HepG2 cells and in liver tissues from patients with HCV, HCV-associated cirrhosis and liver cancer with HCV-associated cirrhosis, using the publicly available E-MTAB-950 with the entire human interactome as the background (30,55). The aim was to discriminate between liver cancer in the presence or absence of viral infection, and to identify the presence and the function of common or specific HUB nodes in the four networks. Although HepG2 cells were revealed to be misidentified as a hepatocellular carcinoma cell line instead a hepatoblastoma by International Cell Line Authentication Committee (iclac.org/databases/cross-contaminations), the results obtained using this cell line were used, as the study focused on a liver cancer cell line without viral infection. In the present study, it was evaluated whether these specific genes in HepG2 cells were components of the circadian rhythm, by identifying their presence in the circadian network. The results of the present study demonstrated that 20/26 HUB genes [CSNK2α1, SH2 domain containing (SRC), ubiquitin D, aurora kinase B (AURKB), cytoskeleton associated protein 5 (CKAP5), replication factor C subunit 4, cell division cycle 20 (CDC20), stratifin, minichromosomemaintenance complex component (MCM)6, checkpoint kinase 1, centromere protein A (CENPA), HLA-B, baculoviral IAP repeat containing 5 (BIRC5), MCM3, mitotic arrest deficient 2 like 1 (MAD2L1), MCM4, ZW10 interacting kinetochore protein (ZWINT), kinesin family member (KIF)2C, inner centromere protein (INCENP) and SPC24 NDC80 kinetochore complex component (SPC24)] were demonstrated to be components of the human circadian network (Fig. 1).

Figure 1. Schematic of the 20 HepG2 HUB genes identified to be involved in the circadian network. With respect to each other in the circadian rhythm network: Red, high degree; yellow, intermediate degree; green, low degree.

The aforementioned 20 HUB genes revealed high degree values in the circadian network, ranging between 287 and 77, indicating that these genes control a large number of metabolic functions and the flow of information via the circadian network. Furthermore, the results of the present study demonstrated that these 20 HUB genes interacted with 31 seed circadian genes, including CLOCK, PER1-3, CRY1-2, ARNTL2, CSNK1D, HDAC4, ZW10, CSNK1E, RFC3, MTR, SFPQ, ESYT1, transferrin receptor, GSK3B, EXOC1, SHANK3 and PLSCR1. In particular, our study has revealed that CLOCK is associated with HUB genes of the HepG2 network via CKAP5, which exhibits a high degree (217), a short path length value of 2.053 and a high value of stress centrality. Thus, CKAP5 interacts via a number of the shortest path-lengths of the network, which makes it a perfect link between circadian and HepG2 networks.

CKAP5 encodes a cytoskeleton-associated protein belonging to the TOG/XMAP215 family, and CKAP is also known as a colonic and hepatic tumor overexpressed gene protein (56,57). Its coded protein serves two distinct functions in spindle formation and in the protection of kinetochore microtubules, via the control of the de-polymerization process and centrosomal microtubule assembly (58,59). These two processes regulate the mitotic cell cycle via spindle formation (60) and the interaction between microtubules and the cell cortex for the directional cell movement. Notably, the results of the present study identified that CKAP5 is associated with three circadian genes (CSNK1E, CSNK1D and HDAC4) and with 10 HepG2 genes (SRC, ZWINT, AURKB, CDC20, CENPA, INCENP, MAD2L1, BIRC5, SPC24 and KIF2C; Fig. 2). Therefore, the two sub-networks of circadian genes and HepG2 genes connected via CKAP5 may be disturbed by any alteration associated with circadian and liver cancer genes. Due to the close associations between nodes, a putative perturbing stressor, for example an alteration of a circadian gene, can induce a perturbation to the global network, and thus, cancer progression (61).

Figure 2. Association between (A) the circadian sub-network and (B) HepG2 HUB-sub-network. Yellow represents CKAP5, the HepG2 HUB-genes and four circadian seed genes (CLOCK, CSNK1E, CSNK1D and HDAC4), which are directly connected via CKAP5. The other genes present in circadian and HepG2 HUB sub-networks are indicated in green and blue, respectively. CKAP5, cytoskeleton associated protein 5; CLOCK, clock circadian regulator; CSNK, casein kinase; HDAC4, histone deacetylase 4.

Structural analysis on common nodes between HepG2 and circadian networks

As all HUB genes identified in our network analysis encode for proteins, it is crucial to understand whether the genes possess specific structural features. Previously, it was demonstrated that the metabolic sub-network specific for liver cancer is formed only by IDPs (30,55). Structural flexibility and binding plasticity enable IDPs to interact with a broad range of molecular partners (30,55). Therefore, in the present study, disorder propensity and the number of molecular partners, with which these proteins can interact, was evaluated using DisProt and BioGRID tools, respectively. As presented in Table II, the proteins encoded by the genes common to the circadian and HepG2 networks belonged to the IDP family. In particular, 3, 11 and 7 proteins exhibited 15%, 16–50% and >50% ID regions (IDRs), respectively. Subsequently, the physical interactions between all proteins encoded by the genes common to the circadian and HepG2 network and other proteins with which they can interact were analyzed, and demonstrated that they are able to form between 17 and 667 interactions. This result supports our hypothesis that the flexibility of the disordered regions functions in the establishment a high number of interactions.

Table II. Molecular properties of the 20 HUB nodes specific for HepG2 that are common to the circadian network.

Table II.

Molecular properties of the 20 HUB nodes specific for HepG2 that are common to the circadian network.

Gene name	Protein code	Protein name	IDR	INT	SEQ
AURKB	Q96GD4	Aurora kinase B	++	268	344
BIRC5	O15392	Baculoviral IAP repeat-containingprotein 5	+++	137	142
CDC20	Q12834	Cell division cycle protein 20 homolog	++	595	499
CENPA	P49450	Histone H3-like centromeric protein A	++	129	140
CHEK1	O14757	Serine/threonine-protein kinase Chk1	++	161	476
CKAP5	Q14008	Cytoskeleton associated protein 5	+++	39	2032
CSNK2A1	P68400	Casein kinase II subunit alpha	++	599	391
HLA-B	P01889	HLA class I histocompatibility antigen, B-7 alpha chain	++	63	362
INCENP	Q9NQS7-INCE human	Inner centromere protein	+++	17	918
KIF2C	Q99661	Kinesin-like protein KIF2C	++	82	725
MAD2L1	Q13257	Mitotic spindle assembly check point protein MAD2A	+	48	205
MCM3	P25205	DNA replication licensing factor MCM 3	++	171	808
MCM4	P33991	DNA replication licensing factor MCM 4	++	152	863
MCM6	Q14566	DNA replication licensing factor MCM 6	++	160	821
RFC4	P35249	Replication factor C subunit 4	+	68	363
SFN	P31947	14-3-3 protein sigma	+++	305	248
SPC24	Q8NBT2	Kinetochore protein Spc24	+++	31	197
SRC	P12931	Proto-oncogene tyrosine-protein kinase Src	++	272	536
UBD	15205	Ubiquitin D	+++	667	165
ZWINT	95229	Zw10- interactor	+	55	277

[i] IDR, intrinsically disordered region; +, IDR <15%; ++, IDR 15–50%; +++, IDR >50%; INT, number of interactions with molecular partners (determined using BioGRID); SEQ, sequence length (number of amino acids); AURKB, Aurora kinase B; BIRC5, Baculoviral IAP repeat-containing protein 5; CDC20, Cell division cycle protein 20 homolog, CENPA, histone H3-like centromeric protein A; CHEK1, serine/threonine-protein kinase Chk1; CKAP5, cytoskeleton-associated protein 5; CSNK2A1, casein kinase II subunit-α; HLA-B, HLA class I histocompatibility antigen, B-7 α-chain; INCENP, inner centromere protein; KIF2C, Kinesin-like protein; MAD2L1, mitotic arrest deficient 2 like 1; MCM3, minichromosome maintenance complex component 3; RFC4, replication factor C subunit 4; SFN, 14–3–3 protein-σ; SPC24, kinetochore protein Spc24; SRC, proto-oncogene tyrosine-protein kinase Src; UBD, ubiquitin D; ZWINT, Zw10 interactor.

Association between miRNAs and genes of the human circadian and HepG2 networks

Following the identification of the associations between circadian and HepG2 networks, and demonstrating that genes common to the two networks encode IDPs, the presence of miRNAs involved in the sub-network between HepG2 and circadian genes was evaluated. A prior study identified the deregulated miRNAs in HepG2 cells compared with normal hepatocytes (51). A total of 11 downregulated (miR-146b-5p, miR-195, miR-122, miR-122a, miR-375, miR-885-5p, miR-768-5p, miR-101, miR-192, miR-194 and miR-215) and 2 upregulated miRNAs (miR-221 and miR-99b) were identified (51). Therefore, in the present study, the target genes of these 13 miRNAs were predicted using the three aforementioned tools. The results of the present study identified 5415 target genes belonging to different metabolic pathways including axon guidance, hippo signaling, endocytosis, phosphoinositide 3-kinase/protein kinase B (AKT) signaling, RAS signaling, RAP1 signaling and chemokine signaling.

To determine the target genes involved in the HepG2 and circadian networks, genes common to the identified 5,415 target genes and the circadian genes were identified. The results of the present study identified 28 genes (APOLD1, BHLHE40, BHLHE41, CLOCK, CRY1, CRY2, extended synaptotagmin 1, ETV5, HDAC4, HLF, ITIH5, KLF11, MECP2, MTR, nuclear receptor subfamily 2 group D member 2, PER1, PER2, PER3, RFC3, SLC2A1, SNTB2, SOCS2, SPRY4, TFRC, TNIP2, UNC13a and ZW10) that are common between circadian genes and the 5415 identified target genes. In addition, four genes (checkpoint kinase 1, KIF2C, MCM6 and MAD2L1), that associated the HepG2 network with the circadian network, were targeted by the following four miRNAs: miR-195-5p, miR-192-5p, miR-122-5p and miR-101-3p. The present study identified three genes (SNTB2, HLF and CRY2) that correlated four miRNAs (miR-195-5p, miR-192-5p, miR-122-5p and miR-101-3p) between them (Fig. 3).

Figure 3. Associations between genes and miRNAs that link HepG2 and circadian networks. In details, four genes (CHECK1, KIF2C, MCM6 and MAD2L1) are reported in red; four miRNAs (miR-195-5p, miR-192-5p, miR-122-5p and miR-101-3p) are reported in green and three genes (SNTB2, HLF and CRY2) that correlate the miRNAs are reported in orange. miR, microRNA; SNTB2, syntrophin β2; MCM6, minichromosome maintenance complex component 6; HLF, HLF PAR bZIP transcription factor; KIF2C, kinesin family member 2C; CRY2, cryptochrome circadian clock 2; MAD2L1, mitotic arrest deficient 2 like 1; CHECK1, checkpoint kinase 1.

Previous studies have suggested that liver cancer is associated with abnormalities in circadian rhythms (26,62) due to alterations in the expression of certain circadian genes in cancerous cells, induced by hypoxia (24), or the overexpression of the mammalian timeless protein (23), a protein that controls chromosome integrity, growth and development. The present study was, to the best of our knowledge, the first to identify a sub-set of HUB genes consisting of genes present in the HepG2 cell network, and involved in cancer progression, in common with human circadian rhythm genes. The results of the present study revealed the following: i) CLOCK is associated, via CKAP5, with the HUB genes of the HepG2 network; ii) CKAP5 is associated with three other circadian genes (CSNK1E, CSNK1D and HDAC4), and with 10HepG2 genes (SRC, ZWINT, AURKB, CDC20, CENPA, INCENP, MAD2L1, BIRC5, SPC24 and KIF2C); iii) the genes linking the circadian system and liver cancer codify for proteins that exhibit IDRs; iv) a sub-panel of seven genes and three miRNAs link human circadian rhythms with liver cancer in a single network.

Sassoni-Corsi et al (63) demonstrated that the liver operates as an exclusive endogenous metabolic reorganizer in tumor-bearing mice (63). Notably, associations between cancer and circadian genes are maintained while the pro-inflammatory response of the liver is altered and leads to the disturbance of AKT, AMP-activated protein kinase and sterol regulatory element binding protein signaling, which, in turn, affects glucose and lipid metabolism. These results demonstrate the requirement to study associations between the circadian rhythms in liver cancer and/or other types of cancer (63). In this context, the results of the present study indicated that further studies are required to determine how the structural perturbation of the HUB nodes in liver cancer may trigger significant and widespread sources of functional changes, which consequently may produce the distinct metabolic functions of cancer cells.

References