Downregulation of ZC3H13 by miR-362-3p/miR-425-5p is associated with a poor prognosis and adverse outcomes in hepatocellular carcinoma

Pan cancer analysis of ZC3H13

In this study, the expression of ZC3H13 in different tumors was initially assessed by the Oncomine database. Among them, ZC3H13 was highly expressed in five cancers, including colorectal cancer, kidney cancer, lymphoma, melanoma and sarcoma. But downregulated in four cancers, brain and CNS cancer, liver cancer, lung cancer, and lymphoma (Figure 1A). Next, to make our analysis more reliable, the expression of ZC3H13 in different cancers was further validated by the GEPIA database. The results showed that ZC3H13 expression was significantly increased in patients with DLBC, ESCA, LAML, PAAD, SKCM, STAD and THYM (Figure 1B). In LIHC, ZC3H13 was significantly decreased (Figure 1C). Taken together, both database results confirmed that ZC3H13 is downregulated in HCC, suggesting that abnormal ZC3H13 may drive the carcinogenesis of HCC.

Figure 1. Expression analysis for ZC3H13 in multiple cancers. (A) mRNA expression levels of ZC3H13 in various types of cancer (Oncomine). The threshold was designed with the following parameters: fold change = 2 and P value = 0.01. The cell number represents the number of datasets that meets the threshold. The color intensity (RED OR BLUE) is directly proportional to the significance level of upregulation or downregulation. (B, C) ZC3H13 expression in TCGA (BRCA, KIRC, LIHC, LUAD, STAD) tissues compared with the corresponding TCGA and GTEx normal tissues. *p value < 0.05.

ZC3H13 can evaluate the prognosis of HCC patients

Next, the prognostic impact of different expression levels of ZC3H13 on patients with BRCA, KIRC, LIHC, LUAD and STAD were analyzed. (Figure 2). The two prognostic indicators were overall survival (OS) and progression-free survival (RFS). If OS is used as an evaluation criterion, ZC3H13 can be used to assess the prognosis of patients with KIRC (P = 5.2E-07) and LIHC (P = 0.013). Again, with RFS assessment, only ZC3H13 is eligible for assessing prognosis in patients with LIHC (P = 0.046). No statistical significance was observed for ZC3H13 in predicting the prognosis in patients with other cancer types. In conclusion, by synthesizing differential expression and prognostic analysis results, ZC3H13 has met the screening criteria and can be used as a candidate factor to assess the prognosis of HCC patients.

Figure 2. Overall survival (OS) and disease-free survival (RFS) analysis for ZC3H13 in various human cancers determined by the GEPIA database. The OS and RFS plots of ZC3H13 in BRCA (A, B). The OS and RFS plots of ZC3H13 in KIRC (C, D). The OS and RFS plots of ZC3H13 in LIHC (E, F). The OS and RFS plots of ZC3H13 in LUAD (G, H). The OS and RFS plots of ZC3H13 in STAD (I, J).

MiR-362-3p/miR-425-5p is an upstream miRNA of ZC3H13

The ceRNA theory is used to describe the mechanism of regulating gene expression. Theoretically, ZC3H13 would be regulated by some ncRNAs; therefore, the upstream miRNAs to which ZC3H13 might be bound were first queried, and 12 miRNAs were identified (Table 1). Cytoscape software was used to establish the miRNA-ZC3H13 regulatory network (Figure 3A). Based on the mechanism by which miRNAs regulate target gene expression, a negative regulatory relationship should exist between miRNAs and ZC3H13; thus, correlation analysis was performed. ZC3H13 was inhibited by upstream miR-362-3p or miR-425-5p (Figure 3B, 3E). No negative correlation was found between ZC3H13 and the other nine predicted miRNAs. miRNAs are upregulated in liver cancer patients and function as oncogenes in terms of prognostic value. The above results are consistent with bioinformatics analysis, supporting the ceRNA hypothesis. Finally, based on the correlation analysis and prognosis analysis results of miR-362-3p/miR-425-5p and ZC3H13, the status of miR-362-3p/miR-425-5p as the upstream regulator of ZC3H13 was determined.

Figure 3. Mir-362-3p/mir-425-5p downregulates ZC3H13 in liver cancer. (A) The miRNA-ZC3H13 regulatory network established by Cytoscape software. (B) ZC3H13 is negatively correlated with miR-362-3p in LIHC analyzed by the starBase database. (C) Mir-362-3p is overexpressed in LIHC analyzed by UALCAN. (D) The prognostic value of miR-362-3p in HCC assessed by Kaplan–Meier plotter. (E) ZC3H13 is negatively correlated with miR-425-5p in LIHC analyzed by the starBase database. (F) Mir-425-5p is overexpressed in LIHC analyzed by UALCAN. (G) The prognostic value of miR-425-5p in HCC assessed by Kaplan–Meier plotter. (H, I) MiR-362-3p/miR-425-5p mimic and inhibitor results detected by qPCR. (J) MiR-362-3p/miR-425-5p mimic and inhibitor results detected by western blot. **P<0.01.

Huh7 and Hep3B cells were transfected with miR-362-3p/miR-425-5p, and the effect on ZC3H13 expression levels was tested. The downregulation of ZC3H13 expression was regulated by miR-362-3p/miR-425-5p mimics and qRT–PCR demonstrated that this reduction could be reversed by miR-362-3p/miR-425-5p inhibitors (Figure 3H, 3I). Furthermore, the protein levels of ZC3H13 were detected by western blotting. And the results were consistent with qRT–PCR results demonstrating that ZC3H13 protein levels were regulated by miR-362-3p/miR-425-5p mimics and inhibitors (Figure 3J). Strikingly, ZC3H13 was identified as a direct target of miR-362-3p/miR-425-5p, further corroborating the RNA-seq results (Figure 3).

Determining the function of ZC3H13 in HCC

To further explore the biological function and pathway analysis of the candidate factor ZC3H13 in HCC, TCGA sequencing data stored on LinkedOmics online tool were utilized. Red dots represent genes that are positively correlated with ZC3H13, and green dots represent negative correlations (Figure 4A, P < 0.05). The genes that were ranked in the top 50 positively and negatively associated with ZC3H13 were presented in a heat map in Supplementary Figure 1. ZC3H13 is involved in biological functions such as phosphoinositide metabolism, Chagas disease, propionate metabolism and transcriptional dysregulation in cancer. Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis showed that these genes were mainly enriched in the JAK-STAT signaling pathway and Hippo signaling pathway.

Figure 4. ZC3H13 coexpression genes in LIHC (LinkedOmics). (A) Global ZC3H13 highly correlated genes identified by the Pearson test in LIHC. Red and green dots represent positively and negatively significantly correlated genes with ZC3H13, respectively. (B) KEGG pathways of ZC3H13-correlated genes. (C) Gene Ontology of ZC3H13-correlated genes. KEGG: Kyoto Encyclopedia of Genes and Genomes; LIHC: Liver hepatocellular carcinoma.

ZC3H13 positively correlates with immune cell infiltration in HCC

Studies have identified zinc finger proteins in other diseases, and they are increasingly recognized as drug targets [14]. As a typical zinc finger protein, ZC3H13 has a finger domain. Tumor-infiltrating immune cells are represented by six types, namely B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils and dendritic cells. Clarifying tumor-infiltrating immune cell types may help to improve the survival prediction of HCC patients increasing the therapeutic effect. Therefore, we investigated the correlation between ZC3H13 and tumor-infiltrating immune cells (Figure 5). The results revealed that the expression of ZC3H13 in HCC patients was significantly and positively correlated with four types of immune cell infiltration, including CD4+ T cells (P = 2.32E-07), macrophages (P = 2.96E-03), neutrophils (P = 1.51E-09) and dendritic cells (P = 1.07E-03). Among these immune cells, tumor-infiltrating neutrophils showed the strongest correlation.

Figure 5. (A) The relationship of immune cell infiltration with ZC3H13 levels in HCC. (B–G) The infiltration level of various immune cells under different expression levels of ZC3H13 in HCC. The correlation of ZC3H13 expression level with B cell (B), CD8+T cell (C), CD4+T cell (D), macrophage (E), neutrophil (F), or dendritic cell (G) infiltration level in HCC.

Correlation analysis of ZC3H13 and immune cell biomarkers in HCC

Continuing the correlation analysis in the previous section, this study analyzed the correlation between ZC3H13 and biomarkers expressed by immune cells in patients with hepatocellular carcinoma. The biomarkers expressed on immune cells are listed in Table 2. For example, CD4 is a CD4+ T cell biomarker. Correlation analysis showed that ZC3H13 was significantly positively correlated with CD4+ T cell biomarkers (R = 0.25, P = 1.60E-06), M1 macrophage biomarkers (R = 0.4, P = 1.20E-15), M2 macrophage biomarkers (R = 0.25, P = 1.10E-06), neutrophil biomarkers (R = 0.25, P = 7.50E-07), and dendritic cell biomarkers (R = 0.26, P = 4.00E-07). These findings prompted us to consider that ZC3H13 may lead to different prognostic outcomes of patients by influencing the immune cells infiltrating hepatocellular carcinoma patients.

Relationship between ZC3H13 and immune checkpoints in HCC

PD1/PD-L1 and CTLA-4 are known to be key immune checkpoints, and they are responsible for tumor immune escape. We evaluated the correlation of ZC3H13 with these immune checkpoints in patients with hepatocellular carcinoma using bioinformatics. The analysis showed that the expression of ZC3H13 was significantly and positively correlated with PD-L1 in HCC (R = 0.437, P = 1.73E-17), which was adjusted by purity (Supplementary Figure 2). Consistent with the TIMER data analysis results, GEPIA database analysis found that ZC3H13 was positively correlated with PD-L1 (R = 0.33, P = 6.3E-11). In addition, no correlation was seen between ZC3H13 and PD1 and CTLA-4.

Predicted upstream regulators of ZC3H13

StarBase (http://starbase.sysu.edu.cn/) searches were performed for microRNA targets using high-throughput CLIP-Seq experimental data and degradome experimental data, providing various visual interfaces to explore microRNA targets. The miRNA-mRNA module was used to analyze upstream miRNAs [23]. StarBase database involves PITA, RNA22, miRmap, microT, miRanda, PicTar, and TargetScan target gene prediction programs. Only target miRNAs validated by more than two prediction programs were considered as miRNAs meeting the screening criteria. The miRNAs screened above were considered candidate miRNAs of ZC3H13. And starBase was used to analyze the correlation between miR-362-3p/miR-425-5p and ZC3H13 in liver cancer patients. Comparison of miR-362-3p/miR-425-5p expression differences between hepatocellular carcinoma patients and normal was also done through this website.

Oncomine analysis

The Oncomine database is a large oncology gene microarray database. Because it covers 65 gene microarray datasets, large-scale gene expression data are recognized and used by everyone. This study used it to analyze gene expression differences. It can also have functional modules such as identifying abnormal values and predicting co-expressed genes [24]. In this study, we analyzed the differences in expression of m6A-related genes in HCC patients and healthy individuals based on mRNA levels using the Oncomine database. Analysis was based on the following criteria: fold change = 2 and P value = 0.01.

GEPIA database analysis

GEPIA (http://gepia.cancer-pku.cn/) integrates TCGA cancer data with GTEx normal tissue data. By revealing driver genes, alleles, differential expression or oncogenic factors, it can provide key clues for our analysis of novel cancer targets and markers. Therefore, it is often used to address differential gene expression in cancer biology [25]. The expression of ZC3H13 in different types of human tumors was detected by GEPIA. In addition, the relationship of ZC3H13 with immune cell surface markers and immune checkpoints in patients with hepatocellular carcinoma was evaluated. The selection criteria were statistically significant with R > 0.1 and a p value < 0.05.

Kaplan–Meier plotter analysis

The Kaplan–Meier plotter (http://kmplot.com) database is processed by a PostgreSQL server that integrates both gene expression and clinical data [26]. Gene expression data and recurrence-free and overall survival information were downloaded from GEO, EGA, and TCGA. To analyze the prognostic value of ZC3H13 and upstream miRNAs, samples of hepatocellular carcinoma patients were divided into two groups with high and low expression. The two cohorts were compared by Kaplan-Meier survival plots, and hazard ratios and log-rank P values were calculated using 95% confidence intervals. Survival time was selected as 60 months.

TIMER database analysis

The TIMER database (https://cistrome.shinyapps.io/timer/) stores a comprehensive resource of immune infiltrations of different cancer types. It is one of the common databases used to analyze tumor immune infiltrates. This study involved six immune cells containing B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendrites [27]. We analyzed the correlation between ZC3H13 and immune cell infiltration levels or immune checkpoint expression levels in liver cancer patients. A p value < 0.05 was considered statistically significant.

Gene set enrichment analysis

LinkedOmics (http://www.linkedomics.org) has over a billion data points. Specifically, it contains the TCGA dataset covering 32 cancer types and stores multi-omic and clinical data on more than 10,000 patients [28]. The “LinkInterpreter” module can translate big data about ZC3H13 into biological functions and enumerate the most likely signaling pathways involved through its inherent pathway analysis and network analysis. It was used for Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis by Gene Set Enrichment Analysis (GSEA). In the “LinkFinder” module of LinkedOmics, we performed a statistical analysis of ZC3H13 co-expression using the Pearson test, and the results are presented in the form of volcano plots. The rank criterion was a false discovery rate (FDR) < 0.05, and there were 500 simulations.

UALCAN database analysis

UALCAN (http://ualcan.path.uab.edu/index.html) is also based on the TCGA database of relevant cancer data. It is often used for online analysis and mining of cancer data [29]. It can perform commands such as expression profiling, survival analysis, correlation analysis, methylation levels and even target gene prediction. Gene expression of miR-362-3p/miR-425-5p in the LIHC TCGA dataset was explored using UALCAN. P < 0.05 indicated statistical significance.

Cell culture, cell transfection, RNA isolation and quantitative real-time PCR (qPCR) analysis

The HCC cell lines used in this study were purchased from the American Type Culture Collection and included Hep3B and Huh7 cells. All cancer cells survived in an environment based on high glucose Dulbecco’s modified Eagle’s medium (DMEM; Thermo Scientific) supplemented with 10% fetal bovine serum (FBS, Gibco), 0.1 mmol/L MEM nonessential amino acids (NEAA; Invitrogen) and 1% L-glutamine (Invitrogen). All the cell lines were cultured under growth conditions of 5% CO₂ in an incubator at 37° C and 100% humidity. The cells were transfected at a density of 2 × 10⁶/ml with the indicated miRNA mimics or inhibitors using Lipofectamine® 3000 (Invitrogen; Thermo Fisher Scientific, Inc.) as indicated by the manufacturer. After 48 hours, the cells were used for subsequent experiments.

Total RNA was extracted from HCC cells using TRIzol® Reagent (Invitrogen; Thermo Fisher Scientific, Inc.) following the manufacturer’s instructions. First, total RNA (1 μg) was reverse transcribed into cDNA using the PrimeScript RT kit with gDNA Eraser (Takara Biotechnology Co., Ltd.) (Experimental conditions: 37° C for 15 minutes and 85° C for 5 seconds). Subsequently, qPCR was performed using TB Green Fast qPCR mix (Takara Biotechnology Co., Ltd.). The following primers were used for qPCR: ZC3H13 forward, 5’- AAAGGAGGGTTTCACCAGAAGTG-3’ and reverse, 5’- CGCTTCGGAGATTTGCTAGAC-3’; and GAPDH forward, 5’- GAAATCCCATCACCATCTTCCAGG-3’ and reverse, 5’- GAGCCCCAGCCTTCTCCATG- 3’. qPCR was performed under the following thermal cycling conditions: initial denaturation at 95° C for 2 min, followed by 40 cycles of 95° C for 20 s, 58° C for 20 s, and 72° C for 20 s. The mRNA expression levels were quantified using the 2-ΔΔCt method and normalized to the internal reference gene GAPDH.

Western blotting

Cells were collected and washed twice with PBS. The cells were then lysed with RIPA lysis buffer containing protein inhibitors on ice for 30 minutes. The cell lysate was then centrifuged at 13 000 rpm for 30 minutes. Supernatants containing protein were collected, and the protein was quantitated using a Pierce BCA Protein Assay Kit (Thermo Fisher Scientific, Waltham, MA). SDS-PAGE involves the principle of separating proteins by electrophoresis according to their size. After electrophoresis, the proteins are transferred to polyvinylidene difluoride (PVDF) membranes. The membranes were blocked with 5% skim milk and shaken at room temperature for 30 minutes. The primary antibody targeting the protein of interest was incubated with the membrane overnight at 4° C with shaking. The next day, an HRP-labeled secondary antibody was used to detect the protein of interest. Finally, the bands were visualized using Pierce ECL (Thermo Fisher Scientific, Waltham, MA). All of the images were captured using the Molecular Imager ChemiDoc XRS+ (Bio–Rad, Hercules, CA).

Statistical analysis

Statistical analyses were automatically calculated from the aforementioned online databases and were divided into the following four aspects: (1) Comparison of mRNA expression using Student’s t test; (2) Correlation of gene expression using Spearman’s correlation to assess or infiltrate immune cells and Pearson’s correlation method; (3) Fisher’s exact test to measure gene enrichment; (4) a p value < 0.05 or a log rank p value < 0.05 was considered statistically significant.

Data availability

The datasets presented in this study can be found in publicly available databases. The names of the databases and their websites have been marked in the article.