Expression And Biological Interaction Network Of RHOC For Hepatic Carcinoma With Metastasis In PBMC Samples

Introduction

Hepatic carcinoma is the third leading cause among all of the cancer-related deaths globally.¹ Although its morbidity and mortality have decreased recently due to some related advances in cancer diagnosis and treatment, majority of the patients remains incurable with an unsatisfactory long-term prognosis once their tumor have become metastatic.^2–4 The existing imaging diagnosis and serum biomarkers detection are highly used for hepatic carcinoma. However, the typically sensitive and specific diagnostic method for the hepatic carcinoma with metastasis are still lacking, especially for the small metastasis or micro-metastasis. Thus, it is essential to explore new effectively diagnostic markers.

Recently, it has been generally accepted that the diversity of biological characteristics determines the invasive and metastatic potentials of hepatic carcinoma, especially for the alteration of gene expression underlying this process.^5,6 Therefore, growing attentions are paid on the analysis of differentially expressed genes (DEGs). Available information indicates that invasion and metastasis is to a large extent attributable to the ability of cell migration.^7,8 Ras homolog family member C (RHOC), a member of the Ras superfamily of small guanosine triphosphatases (GTPases), was regarded as a key molecular which was involved in the process of cell migration.⁹ Increasing evidence presented that up-regulated expression level of RHOC has been significantly linked to increased invasion and metastasis of hepatic carcinoma.^10–16 Besides the RHOC mediate Rho/ROCK signaling pathway which was involved in the process of cell migration,¹⁷ some other biological pathways were also reported to play a role in the hepatic carcinoma with metastasis, such as transforming growth factor-β1 (TGF-β1) induced signaling pathway,¹³ autocrine motility factor signaling pathway,¹⁴ and vascular endothelial growth factor (VEGF) induced angiogenesis signaling pathway.¹⁸ It indicated that there was a complex gene regulation network of RHOC, being significantly linked to the hepatic carcinoma with metastasis.

In order to shed light on these points above, RNA-seq was performed to delineate the profile of DEGs between hepatic carcinoma with metastasis and hepatic carcinoma with solitary tumor in peripheral blood mononuclear cell (PBMC) samples. It was worth noting that PBMC samples had two advantages when compared to the hepatic carcinoma tissues which were used as the main detection sample type in almost all published related studies.^10–16 The first one was that its component was simple. It would decrease the intra-tumor heterogeneity, making the result more reliable and consistent.¹⁹ The second one was that PBMC was not only easily accessible but also minimally invasive, so that it was the ideal detection sample for monitoring tumor progression. Our results showed that through RNA-seq, 469 DEGs between hepatic carcinoma with metastasis and hepatic carcinoma with solitary tumor in PBMC samples were identified. Then, RHOC, which was dysregulated in the tissue samples of hepatic carcinoma^10–16 and had a potential association with tumor metastasis, was selected for further validation by using q-RT-PCR. In addition, the gene regulation network of RHOC linked to hepatic carcinoma with metastasis in PBMC were also drawn by using multi-dimensional bio-informatic analysis methods. Thus, to our knowledge, our results firstly revealed that RHOC in PBMC samples had a potential to be a new easily accessible and minimally invasive diagnostic and minoring bio-marker for hepatic carcinoma with metastasis.

Methods And Materials

Patients

Thirty-three blood samples were collected from patients who visited Zhongda Hospital Affiliated Southeast University (Nanjing, China). Among them, 15 samples were collected from the patients who were diagnosed as hepatic carcinoma with solitary tumor; 18 samples were collected from the patients who were diagnosed as hepatic carcinoma with metastasis. All of them did not receive any treatment and wrote informed consents. This study got ethics approval from the Ethics Committee of Zhongda Hospital Affiliated Southeast University and was conducted in accordance with the Declaration of Helsinki. In this study, the hepatic carcinoma with single node and without intrahepatic or extrahepatic metastasis was defined as a solitary tumor, the hepatic carcinoma with intrahepatic or extrahepatic metastasis was defined as hepatic carcinoma with metastasis.

PBMC Isolating And RNA Extraction

Ficoll-PaqueTM PREMIUM was used to isolate PBMCs, and TRIZOL (Invitrogen, Carlsbad, CA) was used to extract RNA.^19,41 RNA quality was accessed using NanoDrop ND-1000 (Thermo Fisher Scientific, Waltham, MA), and its integrity was determined by RNA integrity number (RIN; Agilent 2100 RIN Beta Version Software).

RNA-Seq And Bio-Informatic Analysis

Twenty-three samples (13 samples of hepatic carcinomas with metastasis were regarded as the cases, and 10 samples of hepatic carcinomas with solitary tumor were regarded as the controls) were chosen for RNA-seq. Their clinical information was showed in our previous study.¹⁹ The preparation of poly-(A) enriched RNA sequencing libraries and the procedures of filtering, mapping, alignment, and DGEs analysis were described in our previous study.¹⁹ The P-value≤0.05 and an absolute value of log2 FC>2 were treated as the thresholds to judge the significance of DEGs.

For the gene expression network analysis, PANTHER (protein annotation through evolutionary relationship) classification system (http://www.pantherdb.org/) were used to perform GO and pathway classification analysis.⁴² STRING v 11.0 (http://string-db.org/cgi/input.pl) was used for constructing PPI network and performing Reactome pathway enrichment analysis.⁴³ cBio Cancer Genomics Portal (http://cbioportal.org), an open-access resource for interactive exploration of multidimensional cancer genomics data sets, currently containing 225 cancer studies,⁴⁴ was used to analyze the types and frequency of gene alterations and biological interaction network with neighboring genes in this study. The LinkFinder module of LinkedOmics database (http://www.linkedomics.org/login.php), a web-based platform for analyzing 32 TCGA cancer-associated multi-dimensional datasets, was used to select the co-expression genes of RHOC in TCGA database.⁴⁵ StarBase (http://starbase.sysu.edu.cn/) is designed for decoding the interaction networks of lncRNAs, miRNAs, ceRNAs, RNA-binding proteins (RBPs) and mRNAs supported by multi-dimensional sequencing data, including cross-linking immunoprecipitation (CLIP)-seq, degradome-seq and RNA-RNA interactome data.^46,47 StarBase v3.0 was updated in 2018. It identified more than 1.1 million miRNA-ncRNA, 2.5 million miRNA-mRNA, 2.1 million RBP-RNA, and 1.5 million RNA–RNA interactions. And its Pan-Cancer Analysis Platform is designed for decoding Pan-Cancer Networks of lncRNAs, miRNAs, pseudogenes, snoRNAs, RBPs and all protein-coding genes by analyzing their expression profiles across 32 cancer types (~10,000 RNA-seq and ~9,900 miRNA-seq samples) integrated from TCGA project. In this study, the RNA-ceRNA module and RNA-RNA module of StarBase v3.0 were used to confirm the co-expression genes with RHOC in TCGA database. Pan-Cancer module was used to analyze overall survival for LIHC. When the FDR or P-value was less than 0.05, the result was significant.

q-RT-PCR

Ten PBMC samples (5 hepatic carcinomas with metastasis were regarded as the cases, and 5 hepatic carcinomas with solitary tumor were regarded as the controls) were used for qRT-PCR. The primers were synthesized by Takara Bio, Inc and showed in Table 1. The procedure of q-RT-PCR was the same as that was previously reported.¹⁹ Two percent agarose gel electrophoresis and dissociation curve were used to access the quality of its amplification products (Figures S3 and S4).

Table 1 Sequences Of Primers

Statistical Analysis

The data were analyzed with MATLAB^® (version 2010b). Statistical differences were examined by a t-test. The relationship between the expression level of RHOC and T staging was analyzed by Spearman test. When P-values were less than 0.05, it was regarded as significant.

Results

DEGs Profiling For Hepatic Carcinoma With Metastasis Vs Hepatic Carcinoma With Solitary Tumor In PBMC Samples

Through RNA-seq, 469 DEGs were identified in PBMC samples of hepatic carcinoma with metastasis and hepatic carcinoma with solitary tumor. Among them, 352 DEGs were up-regulated, and 117 DEGs were down-regulated (Figure 1 and Table S1). Their functional classification was performed by PANTHER Classification System. As Figure 2 shown, these DEGs were the main component of cell (gene ontology (GO): 0005623), organelle (GO: 0043226), extra-cellular region (GO: 0005576), membrane (GO: 0016020), and protein-containing complex (GO: 0032991). They played a role in cellular process (GO: 0009987), such as binding (GO: 0005488), catalytic activity (GO: 0003824), etc. Their protein–protein interaction (PPI) network performed by STRING V11.0 is shown in Figure 3 (Number of nodes: 454; Number of edges: 1362; Average node degree: 6; Avg. local clustering coefficient: 0.398; Expected number of edges: 728; PPI enrichment P-value: <1.0E−16) and the Reactome pathway enrichment analysis was shown in Table S2.

Figure 1 Vocal plot of DEGs identified in PBMC samples for hepatic carcinoma with metastasis vs hepatic carcinoma with solitary tumor. Notes: 469 DEGs were identified by RNA-seq. Among them, 352 genes were up-regulated, and 117 genes were down-regulated.

Figure 2 Functional classification analysis of DEGs for hepatic carcinoma with metastasis in PBMC samples (PANTHER). Notes: (A) Classification biological process; (B) classification cellular component; (C) classification pathway; (D) classification molecular function; (E) classification pathway P00031. Red triangles indicated the items RHOC were involved in.

Figure 3 PPI network of DEGs for hepatic carcinoma with metastasis in PBMC samples (STRING). Notes: Number of nodes: 454; number of edges: 1362; average node degree: 6; Avg. local clustering coefficient: 0.398; expected number of edges: 728; PPI enrichment P-value: <1.0E−16. Thicker line indicates increased edge confidence.

Previous studies showed that RHOC was significantly involved in the process of metastatic hepatic carcinoma.^10–16 Our RNA-seq assay got the same result, which presented that the expression of RHOC was significantly higher in the PBMC samples of hepatic carcinoma with metastasis than in those of hepatic carcinoma with solitary tumor (Log2 fold change (FC): −2.42498, P-value=0.005; Table S1). And through functional classification analysis, RHOC was found to be involved in Ras Pathway (P04393), Rho (P00860) of Inflammation mediated by chemokine and cytokine signaling pathway (P00031), heterotrimeric G-protein signaling pathway-Gq alpha and Go alpha-mediated pathway (P00027), cytoskeletal regulation by Rho GTPase (P00016), Axon guidance mediated by semaphorins (P00007) or Slit/Robo (P00008), angiotensin II-stimulated signaling through G proteins and beta-arrestin (P05911), and Angiogenesis (P00005) (Figure 2C and E). In addition, the analysis of overall survival for liver hepatocellular carcinoma (LIHC) of cancer genome atlas (TCGA) database was performed by StarBase v3.0. The result showed that RHOC was significantly associated with poor prognosis of hepatic carcinoma (Figure 4, P-value=0.0025).

Figure 4 Overall survival for LIHC in TCGA (StarBase). Notes: Cancer number was 369. Coefficient was 0.55. Hazard ratio was 1.73. Log-rank P-value was 0.0025.

Validation Of RHOC By q-RT-PCR

In order to further validate the expression of RHOC, another 10 PBMC samples were collected and used to perform q-RT-qPCR. Their baseline characteristics are shown in Table 2. The result showed that RHOC was significantly higher in PBMC samples of hepatic carcinoma with metastasis than in those of hepatic carcinoma with solitary tumor (7.15±0.88 vs 8.66±0.62, P-value=0.014, Figure 5), which was the same with the results of RNA-seq. In addition, a significantly negative association was observed between the expression level of RHOC and T staging of hepatic carcinoma with metastasis (Spearman coefficient r=−0.949, P-value=0.014).

Table 2 The Baseline Characteristics Of The Samples That Were Used To Perform qRT-PCR

Figure 5 The differential expression of RHOC detected by qRT-PCR. Notes: The hepatic carcinoma with metastasis 7.15±0.88 vs hepatic carcinoma with solitary tumor 8.66±0.62, P=0.014. Data were represented as mean ± standard deviation.

Functional Network Analysis Of RHOC For Hepatic Carcinoma With Metastasis In PBMC Samples Based On Bio-Informatic Method

Frequency And Type Of RHOC Alterations In TCGA Database

Firstly, cBioPortal platform was used to describe the types and frequency of RHOC alterations in LIHC of TCGA database. As shown in Figure 6A, RHOC was altered in 24 of 377(6.4%) LIHC patients. These alterations were mRNA upregulation in 23 cases (6.1%) and mutation in 1 case (0.3%) (Table S3). Thus, the high expression is the most common type of RHOC alterations in LIHC, which altered in 12.5% of 8 hepatocellular carcinoma plus intrahepatic cholangiocarcinoma and 5.13% of 22 hepatocellular carcinoma (Figure 6B). Besides, the expression level of RHOC was significantly negatively associated with gene methylation (Figure 6C), indicating that gene methylation alteration affected its gene expression level.

Figure 6 RHOC alterations in TCAG LIHC (cBioPortal). Notes: (A) The OncoPrint provided an overview of genomic alterations in RHOC affecting individual patients (columns) in LIHC from the TCGA. The different types of genetic alterations were highlighted in different colors. *RHOC was altered in about 6% of LIHC patients. (B) RHOC alterations in different cancer types in LIHC of TCGA database. RHOC mRNA high altered in 12.5% of 8 hepatocellular carcinoma plus intrahepatic cholangiocarcinoma and 5.13% of 22 hepatocellular carcinoma. Mutation altered in 0.23% of 1 hepatocellular carcinoma. (C) The association between RHOC gene methylation and gene expression. P-value for Spearman test was 1.86E−7; P-value for Pearson test was 1.359E−6.

Biological Interaction Network Of RHOC With Its Neighboring Genes

Next, the tab Network in cBioPortal was used to identify RHOC-neighboring genes and determine the biological interaction network of them. When the alteration frequency was more than 20%, the neighbor gene was selected. The results were shown in Figure S1 and Table S3. RHOC-neighboring genes, including protein tyrosine kinase 2 gene (PTK2), phosphatidylinositol 4-phosphate 5-kinase type-1 alpha gene (PIP5K1A), Rho GTPase activating protein 39 gene (ARHGAP39), and Rho guanine nucleotide exchange factor 11 gene (ARHGEF11), were selected and their biological interaction network was determined. All of them presented up-regulated alterations in LIHC. Among them, PIP5K1A, ARHGAP39, and ARHGEF11 controlled the reaction that changed the state of RHOC. PANTHER GO-Slim and Pathway classification analysis showed that PTK2 was the component of cell junction (GO: 0030054), and involved in biological adhesion (GO: 0022610), biological regulation (GO: 0065007), multicellular organismal process (GO: 0032501), angiogenesis (P00005), cholecystokinin receptor (CCKR) signaling map (P06959), Gonadotropin-releasing hormone receptor pathway (P06664), integrin signaling pathway (P00034), and VEGF signaling pathway (P00056); PIP5K1A and ARHGAP39 were involved in catalytic activity (GO: 0003824); ARHGAP39 was involved in molecular function regulator (GO: 0098772); PTK2 and PIP5K1A were involved in metabolic process (GO: 0008152) and cellular process (GO: 0009987). However, none of RHOC-neighboring genes was found to be altered in our RNA-seq assay.

Biological Interaction Network Of RHOC With Its Co-Expression Genes

10,031 co-expression genes of RHOC were identified in LIHC of TCAG database by LinkedOmic platform, of which 232 genes were found to be altered in our RNA-seq assay. Through combination with the analysis performed by StarBase v3.0 platform, 15 co-expression genes were identified as the competing endogenous RNAs (ceRNAs) with RHOC and 4 co-expression genes presented RNA-RNA interactions with RHOC. The 15 co-expression genes were DDB1- and CUL4-associated factor 12 gene (DCAF12), peptidyl-prolyl cis-trans isomerase gene (FKBP8), sodium channel subunit beta-1 gene (SCN1B), forkhead box protein O3 gene (FOXO3), GAS2-like protein 1 gene (GAS2L1), tumor-specific transplantation antigen 3 gene (TSTA3), myosin regulatory light polypeptide 9 gene (MYL9), mitogen-activated protein kinase kinase 3 gene (MAP2K3), tensin-1 gene (TNS1), tetraspanin-9 gene (TSPAN9), EH domain-containing protein 3 gene (EHD3), activating signal cointegrator 1 complex subunit 2 gene (ASCC2), protein-glutamine gamma-glutamyltransferase 2 gene (TGM2), platelet endothelial aggregation receptor 1 gene (PEAR1), and Bcl-2-like protein 1 gene (BCL2L1). The 4 interactional genes with RHOC were glutathione peroxidase 1 gene (GPX1), FKBP8, transcription factor Dp-1 gene (TFDP1), and apolipoprotein L2 gene (APOL2). Their frequency and type of alterations and expression levels in LIHC of TCAG database were shown in Table S4 and Figure S2. And STRING v11.0 was used to present their biological interaction network (Figure 7, Tables 3 and 4). The results showed that the interaction between MYL9 and RHOC had the highest edge evidence (Score: 0.958), and through Reactome pathway enrichment analysis, they were found to be significantly involved in RHO GTPases activate ROCKs Pathway (HSA-5627117), RHO GTPases activate citron kinase (CIT) Pathway (HSA-5625900), RHO GTPases activate protein kinases Ns (PKNs) Pathway (HSA-5625740), and Sema4D induced cell migration and growth-cone collapse Pathway (HSA-416572) (False discovery rates (FDRs) were 0.0257, 0.0257, 0.0414, and 0.0257).

Table 3 The Interactions Of The Co-Expression Genes With RHOC For Hepatic Carcinoma With Metastasis In The Biological Network (STRING)

Table 4 The Reactome Pathways Enrichments In The Biological Network Of The Co-Expression Genes With RHOC For Hepatic Carcinoma With Metastasis

Figure 7 Biological interaction network of RHOC with its co-expression genes (STRING). Notes: Number of nodes: 19; number of edges: 8; average node degree: 0.842; Avg. local clustering coefficient: 0.342; expected number of edges: 2; PPI enrichment P-value: 0.000516. The disconnected nodes were hided. Thicker line indicated increased edge confidence. Nodes in green were involved in RHO GTPases activate ROCKs Pathway (HSA-5627117). Nodes in blue were involved in RHO GTPases activate CIT Pathway (HSA-5625900). Nodes in red were involved in Sema4D induced cell migration and growth-cone collapse Pathway (HSA-416572). Nodes in purple were involved in RHO GTPases activate PKNs Pathway (HSA-5625740). Nodes in yellow were involved in Intrinsic Pathway for Apoptosis (HSA-109606).

Discussion

Although there were some advances in hepatic carcinoma diagnosis and treatment, majority of the patients remains incurable with an unsatisfactory long-term prognosis, especially for hepatic carcinoma with metastasis.^2–4 Existing clinical imaging diagnosis and serum markers lack adequate sensitivity and specificity for small metastasis or micro-metastasis, which lead to diagnostic delays or omissions to threaten people’s life. Therefore, seeking effectively diagnostic and course monitoring biomarkers is still essential for hepatic carcinoma with metastasis.

Recently, increasing related studies have focused on this point. And it is generally accepted that transcellular migration of tumor cells through host structures is crucial in cancer invasion and metastasis. RHOC has been reported to play an important role in cell migration. It has attracted interest with its increased expression being linked to increased invasion, metastasis, and poor prognosis in various cancers, such as breast cancer, melanoma, pancreatic cancer, colon cancer, bladder cancer, non-small cell lung carcinoma, etc.^20–22 Recently, RHOC has been also reported to be significantly involved in the invasion and metastasis of hepatic carcinoma tissues.^10–16 However, the hepatic carcinoma tissue which was used to evaluate RHOC expression level in majority of the studies was not seemed to be the ideal detection sample type for clinical course supervision. We downloaded the data from TCGA (Table S5). Spearman rank correlation test was performed, and the result showed that the expression level of RHOC was significantly associated with the tumor purity of the tissue samples (Spearman coefficient was −0.138. P-value=0.011). It presented that the components of tumor tissue affected the expression level of RHOC and the heterogeneity of tumor tissues could cause inconsistencies among different studies. Thus, in this study, PBMC samples were used to identify the DEGs of hepatic carcinoma with metastasis, because it was not only being easily accessible and minimally invasive but also could benefit the reliability of the results.¹⁹ Our results of DEG screening assay by RNA-seq and validation assay by qRT-PCR showed that the expression level of RHOC in PBMCs was significantly higher in hepatic carcinoma with metastasis than in hepatic carcinoma with solitary tumor, which was consistent with previous studies.^10–16 Overall survival analysis was performed in LIHC patients of TCAG database, indicating that up-regulated RHOC was significantly linked to the poor prognosis of hepatic carcinoma. Therefore, our results suggested that RHOC in PBMC samples could be a potentially minimally invasive biomarker for the diagnosis and clinical course supervision of hepatic carcinoma with metastasis.

In order to better understand the role of RHOC in hepatic carcinoma with metastasis, the multidimensionally functional and biological interaction network analysis was performed. Eighteen co-expression genes with RHOC were identified, including 15 ceRNAs and 4 interactional genes with RHOC. Among them, TGM2 were reported to be dysregulated in hepatic carcinoma with metastasis by previous studies. Yamaguchi H et al presented that TGM2 was upregulated in the early intrahepatic metastasis of hepatic carcinoma and might contribute to early hepatic carcinoma metastasis through signaling pathways unrelated to epithelial–mesenchymal transition (EMT) and integrin signaling.²³ PPI network was drawn to present the interactions between these co-expression genes and RHOC. Among them, the interaction between MYL9 and RHOC had the highest edge evidence (Score: 0.958). MYL9 was myosin regulatory subunit that played an important role in regulation of both smooth muscle and non-muscle cell contractile activity via its phosphorylation and implicated in cytokinesis, receptor capping, and cell locomotion. Reactome pathway enrichment analysis showed that MYL9 and RHOC were significantly involved in RHO GTPases activate ROCKs Pathway (HSA-5627117), RHO GTPases activate CIT Pathway (HSA-5625900), RHO GTPases activate PKNs Pathway (HSA-5625740), and Sema4D induced cell migration and growth-cone collapse Pathway (HSA-416572). The interactional network of these pathways is shown in Figure 8. RHO GTPases activate ROCKs Pathway (HSA-5627117) which was reported to play an important role in the hepatic carcinoma with metastasis was involved in Sema4D induced cell migration and growth-cone collapse Pathway (HSA-416572). Sema4D-mediated plexinB1 activation activated RHOA/RHOB/RHOC and its downstream effector ROCK.^24–27 ROCK then phosphorylated myosin regulatory light chains (MRLC), including MYL9 and MYL12B, to induce actomyosin stress fiber contraction and to direct the assembly of focal adhesion complexes and integrin-mediated adhesion.^26–30 The phosphorylation of MRLC also could be promoted by RHO GTPases activate PKNs Pathway (HSA-5625740) and RHO GTPases activate CIT Pathway (HSA-5625900). RHO GTPases RHOA, RHOB, RHOC and Ras-related C3 botulinum toxin substrate 1 (RAC1) bound to PKN1, PKN2, and PKN3, bringing them in proximity to the PIP3-activated co-activator 3-phosphoinositide-dependent protein kinase 1 (PDPK1).^31,32 PDPK1 phosphorylated PKNs on a highly conserved threonine residue in the kinase activation loop.^31,32 Activated PKN1 (p-T774-PKN1) phosphorylated protein phosphatase 1 regulatory subunit 14A (PPP1R14A), which could subsequently phosphorylate on threonine T38 binds MLCP (myosin light chain phosphatase complex) and inhibit its catalytic activity to increase MRLC phosphorylation.^33–38 RHO GTPases RHOA, RHOB, RHOC, and RAC1 also could bind to CIT, which might be activated through autophosphorylation to phosphorylate the MRLC.^39,40

Figure 8 The interactional network of Reactome pathways that MYL9 and RHOC involved in. Notes: RHO GTPases activate ROCKs Pathway (HSA-5627117), RHO GTPases activate PKNs Pathway (HSA-5625740), and RHO GTPases activate CIT Pathway (HSA-5625900) were involved in Sema4D induced cell migration and growth-cone collapse Pathway (HSA-416572), through phosphorylating MRLC. Row indicated promotion, and “T” indicated inhibition.

Conclusion

To our knowledge, this study first presented that RHOC in PBMCs could be a potentially minimally invasive biomarker for the diagnosis and clinical course supervision of hepatic carcinoma with metastasis. In addition, the PPI network showed that MYL9 and RHOC had the highest edge evidence, and were involved in some cell migration-related pathways. Although this biological network needed to be further studied in the future, it would lay a foundation for the follow-up studies of the role of RHOC in tumor invasion and metastasis.