- Research article
- Open Access
Identification of a ZC3H12D-regulated competing endogenous RNA network for prognosis of lung adenocarcinoma at single-cell level
BMC Cancer volume 22, Article number: 115 (2022)
To identify hub genes from the competing endogenous RNA (ceRNA) network of lung adenocarcinoma (LUAD) and to explore their potential functions on prognosis of patients from a single-cell perspective.
We performed RNA-sequencing of LUAD to construct ceRNA regulatory network, integrating with public databases to identify the vital pathways related to patients’ prognosis and to reveal the expression level of hub genes under different conditions, the functional enrichment of co-expressed genes and their potential immune-related mechanisms.
ZC3H12D-hsa-miR-4443-ENST00000630242 axis was found to be related with LUAD. Lower ZC3H12D expression was significantly associated with shorter overall survival (OS) of patients (HR = 2.007, P < 0.05), and its expression was higher in early-stage patients, including T1 (P < 0.05) and N0 (P < 0.05). Additionally, ZC3H12D expression was higher in immune cells displayed by single-cell RNA-sequencing data, especially in Treg cells of lung cancer and CD8 T cells, B cells and CD4 T cells of LUAD. The functional enrichment analysis showed that the co-expressed genes mainly played a role in lymphocyte activation and cytokine-cytokine receptor interaction. In addition, ZC3H12D was associated with multiple immune cells and immune molecules, including immune checkpoints CTLA4, CD96 and TIGIT.
ZC3H12D-hsa-miR-4443-ENST00000630242 ceRNA network was identified in LUAD. ZC3H12D could affect prognosis of patients by regulating mRNA, miRNA, lncRNA, immune cells and immune molecules. Therefore, it may serve as a vital predictive marker and could be regarded as a potential therapeutic target for LUAD in the future.
Lung adenocarcinoma (LUAD), accounting for approximately 60% of non-small cell lung cancer (NSCLC), is the most common subtype of lung cancer with a high incidence worldwide [1, 2]. With the advancement of technological innovation, the mechanisms behind LUAD are gradually revealed. Growing evidence highlights the vital role of the competitive endogenous RNA (ceRNA) regulatory networks in LUAD. For instance, AC079160.1-miR-539-5p-CENPF axis may participate in hypoxia-induced tumor cell stemness of LUAD. Low expression of AC079160.1 and CENPF and high expression of miR-539-5p were correlated with hypoxia and stemness index, indicating a better prognosis . Furthermore, linc01833-miR-519e-3p-S100A4 axis might be associated with LUAD progression, and linc01833 overexpression can significantly improve proliferation and invasion ability of lung cancer cells as well as promote the epithelial-mesenchymal transformation process .
Although the function of ceRNA network concerning LUAD has been uncovered gradually, the role of single-cell RNA sequencing (scRNA-seq) in LUAD is still waiting to be explored. ScRNA-seq is a novel hot spot to demonstrate the underlying mechanisms of LUAD, uncovering new differentially expressed genes as well as the heterogeneity of immune response-related genes [5, 6]. Besides, scRNA-seq was applied to unravel the molecular and cellular reprogramming mechanisms in metastatic LUAD and the cell-cycle state of the circulating tumor cells of cerebrospinal fluid in LUAD patients with leptomeningeal metastases [7, 8]. When it comes to immunotherapy, scRNA-seq technology might locate the expression of hub genes in immune cells of LUAD to improve precision of immunotherapy in the future . Therefore, integrating ceRNA network with scRNA-seq data may provide a promising strategy for further understanding the underlying mechanisms of LUAD, and might collect more valuable knowledge to improve individual targeted treatment.
At present, there is no investigation on zinc finger CCCH-type containing 12D (ZC3H12D)-hsa-miR-4443-ENST00000630242 axis, which might play a critical role in LUAD. Anti-oncogene ZC3H12D, also called p34, is a member of the zinc finger CCCH-type protein family that is associated with gene expression such as IER3, TNF, IL-6, NF-κB and TLR, degrading inflammatory transcripts and attenuating macrophage response [10,11,12]. Previous studies have demonstrated that ZC3H12D could be targeted by miR-128-3p, involving in cell proliferation and migration in osteosarcoma . However, the functions of hsa-miR-4443 in different cancers are heterogeneous. It could promote the drug resistance to epirubicin in NSCLC ; While its overexpression of hsa-miR-4443 acts in a tumor-suppressive manner, decreas the invasiveness of hepatocellular carcinoma (HCC) , glioblastoma (GBM)  and colorectal cancer (CRC) . Next, long non-coding RNAs (lncRNAs) have been demonstrated their roles in ceRNA networks as regulators, participating in the regulation of various pathological processes related to cancers. Compared with the mRNA and miRNA, little is known about the function of lncRNA FAM30A. Highly expressed in B cells, FAM30A is correlated with the regulation of immune response and immunoglobulin genes [18, 19]. Therefore, our study aims to reveal the function of the ZC3H12D-hsa-miR-4443-ENST00000630242 axis in the prognosis of LUAD, especially the role of ZC3H12D concerning immunomics.
Patients and clinical samples
Ten paired LUAD and paracancerous tissues were collected between January 1, 2019 and May 31, 2019 at Fujian Medical University Second Affiliated Hospital. All samples were obtained from LUAD patients who only received primary surgical treatment. Patient characteristics are provided in Table 1. All the Surgically-resected samples were immediately snap-frozen by liquid nitrogen, and then transferred to a − 80 °C refrigerator for RNA extraction. According to the World Health Organization (WHO) guidelines (2015), we had two well-experienced pathologists to confirm the clinicopathological diagnosis. The study was approved by the bioethical committees at The Second Affiliated Hospital of Fujian Medical University, China (2020-206). And, all participating patients provided written informed consent.
RNA extraction and sequencing
Using the RNeasy Mini Kit (Qiagen, Germany) to isolate the total RNA of LUAD tissues and paracancerous tissues from collected frozen tissues following the standard manufacturer’s instructions. Qubit 4.0 (Thermo Fisher Scientific, Wilmington, DE, USA) was used to evaluate the RNA concentration, and agarose gel electrophoresis was applied to assess the RNA quality.
Then, the maximum residual non-coding RNA (ncRNA) was retained after the ribosomal RNA was removed from the total RNA. Using the TruSeq RNA Sample Prep Kit (Illumina, San Diego, CA, USA) to perform the cDNA library construction after fragments of rRNA-depleted RNA. Following the standard manufacturer’s instructions, the VAHTS total RNA-seq Library Prep kit for Illumina (Vazyme NR603, China) was used for generating lncRNA/mRNA sequencing libraries. The 150-bp paired-end reads’ cDNA fragments were generated for RNA sequencing. Then, establish the miRNA library for samples using the NEBNext® Multiplex Small RNA Library Prep Set for Illumina® (NEB). With a single lane of Illumina HiSeq Xten sequencing platform, 12 libraries were pooled and sequenced. And establish the miRNA library with 50-bp paired-end reads using the Illumina’s TruSeq small RNA library preparation kit. Using the Illumina HiSeq Xten platform to carry out the sequencing for both lncRNA/mRNA and miRNA after library construction.
Identifications of differentially expressed mRNAs, miRNAs and lncRNAs
Mirdeep2 (v184.108.40.206) was applied to predict new miRNA, whose expression was calculated and standardized using counts per million (CPM) read . While, LncRNAs were annotated by these databases, including CNCI (https://github.com/www-bioinfo-org/CNCI) , CPC2 (http://cpc2.cbi.pku.edu.cn/) , CPAT (https://sourceforge.net/projects/rna-cpat/)  and PLEX (https://sourceforge.net/projects/plek/) . The intersection was retained for further analysis.
Screen the differentially expressed mRNAs (DEmRNAs), differentially expressed miRNAs (DEmiRNAs) and differentially expressed lncRNAs (DElncRNAs) between LUAD and normal tissues using the DESeq2Rpackage (https://bioconductor.org/packages/release/bioc/html/DESeq2.html) in the Bioconductor project. The cut-off criteria was |Log2(fold change)|(|log2FC|) ≥1 and statistical P < 0.05. Subsequently, unsupervised hierarchical clustering was performed for DE-RNAs https://cran.r-project.org/web/packages/pheatmap/index.html.
Analysis of the DElncRNAs enrichment pathway
Using the gene ontology (GO, http://www.bioconductor.org/packages/release/bioc/html/topGO.html) function analysis to screen enrichment of targeted genes to annotate the biological functions regulated by DElncRNAs, including molecular function (MF), biological processe (BP), and cellular component (CC). In addition, using the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/keggbin/show_organism?menu_type=pathway_maps&org=hsa) analysis to determine the vital signaling pathways associated with DElncRNAs. Both GO and KEGG enrichment analysis set gene count ≥2 and P value < 0.05 as the threshold for statistical significance.
Predication of miRNA regulation relationship
The miRWalk 2.0 (http://zmf.umm.uni-heidelberg.de/pps/zmf/mirwalk2/) was applied to perform the prediction of miRNA-gene analysis of DEmiRNA . Additionally, using the miRWalk, miRanda, miRDB, miRMap and TargetScan databases to predict the potential DEmiRNA-DEmRNA regulatory relationships. Subsequently, the StarBase (http://starbase.sysu.edu.cn/) database was used to predict the potential miRNA-lncRNA regulatory relationships by DEmiRNAs . Then, based on shared DEmiRNAs with which DEmRNAs and DElncRNAs interact, the DEmRNA-DEmiRNA and DEmiRNA-DElncRNA regulatory relationships were successfully constructed, and were visualized by Cytoscape software.
Construction of lncRNA-miRNA-mRNA ceRNA network
Based on the hypothesis of miRNA sponge, we focused on the positive correlation expression of DElncRNAs-DEmRNAs, and obtained the co-expressed relationships between DEmRNAs and DElncRNAs simultaneously regulated by DEmiRNAs. Subsequently, based on shared miRNAs, we constructed the competing endogenous RNA (ceRNA) network. Furthermore, we applied a hypergeometric cumulative distribution function test to predict the possible ceRNA pairs, and only the pairs with correlation coefficient>0.5 and P value <0.05 were selected. The Sankey diagram was built based on the R software package ‘ggalluval’ (https://github.com/corybrunson/ggalluvial).
Public data source
Raw counts of RNA-sequencing data and corresponding clinical information were obtained from The Cancer Genome Atlas (TCGA) dataset (https://portal.gdc.cancer.gov/) in January 2020, and the method of acquisition and application for data were complied with the guidelines and policies. LUAD-related RNA-sequencing data and corresponding clinical information were obtained from the TCGA database as following criteria: 1) histologically diagnosed as LUAD; 2) data with clinical information. Finally, a total of 510 paired LUAD tissues was included for further analysis.
MiRNAs data was downloaded from KM-plotter (https://kmplot.com/analysis/) , and LncRNAs data was obtained from TCGA database. The single-cell sequencing data of normal lung tissue from 8 mice came from Tabula Muris (https://tabula-muris.ds.czbiohub.org/)  https://www.synapse.org/. While, the single-cell sequencing data of NSCLC, including LUAD, resulted from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) and European Molecular Biology Laboratory (EMBL, https://www.ebi.ac.uk/). The databases were included only if they contained single cell RNA sequencing of NSCLC. Finally, we founded 6 valuable datasets. The ScRNA-seq data of EMTAB6149 was prepared from 5 dissected lung tumors , while GSE117570 contained scRNA-seq data of tumor-infiltrating immune cells from 4 untreated NSCLC patients . Besides, GSE127465 contained scRNA-seq data of red blood cell (RBC)-depleted cells from NSCLC tumor and from blood of 7 patients, as well as CD45-positive cells from lungs from two healthy mice and two tumor-bearing mice . GSE127471 data were collected from cryopreserved peripheral blood mononuclear cells of NSCLC and GSE131907 contained data of 208,506 cells derived from 44 patients with LUAD, taken from normal lung tissues, LUAD tissues, normal lymph nodes, invaded lymph nodes, pleural effusion, and brain metastases [32, 33]. Additionally, GSE139555 contained data of pretreatment samples from 14 NSCLC patients, which covered normal tissues, primary tumors and peripheral blood . The present study meets the criteria of data usage and publishing of the National Cancer Institute of National Institutes of Health.
Survival prognosis and clinical factors
The Kaplan–Meier survival analysis was used to compare the survival difference between the above two groups using ‘survival’ (https://cran.r-project.org/web/packages/survival/index.html) and ‘survminer’ (https://cran.r-project.org/web/packages/survminer/index.html) R packages. When processing mRNA or lncRNA expression data, we splited patients by median. However, due to zero miRNA expression data in some LUAD samples, we decided to use auto select best cutoff to group patients better. TimeROC analysis was performed to compare the predictive accuracy of each gene and risk score. For Kaplan–Meier curves, p-values and hazard ratio (HR) with 95% confidence interval (CI) were generated by log-rank tests. Univariate and multivariate cox regression analysis were performed to identify the proper terms to build the nomogram. The forest plot was used to display P values, HR and 95% CI of each variable using ‘forestplot’ R package (https://CRAN.R-project.org/package=forestplot). Nomograms were developed to predict LUAD patients' overall survival (OS) at 1, 3, 5 years, respectively, based on the results of multivariate Cox proportional hazards analysis. The nomograms quantified the risk so that we could evaluate the prognosis of LUAD patients by the points associated with each risk factor through ‘rms’ R package (https://cran.r-project.org/web/packages/rms/index.html) . Both the above analysis methods and R packages were implemented by R foundation for statistical computing (2020) version 4.0.3 and ggplot2 (v3.3.2). P value< 0.05 was considered statistically significant.
Gene expression atlas based on single-cell RNA sequencing data
We chose the t-distributed Stochastic Neighbor Embedding (t-SNE) algorithm [36, 37] or Uniform Manifold Approximation and Projection for Dimension Reduction (UMAP)  to reduce the dimensionality of the quality-controlled scRNA-seq data of normal lung cell from mouse . Moreover, scRNA-seq data from lung cell atlas of human  and lung cancer brain metastases  of human were re-analyzed through the UCSC cell browser (https://cells.ucsc.edu/) . To understand the expression of ZC3H12D in different cell types across selected datasets, we applied TISCH (http://tisch.comp-genomics.org/) to obtain the ZC3H12D average expression data from multiple databases . Furthermore, we explored the expression of ZC3H12D in LUAD and other well-characterized NSCLC at single-cell level and identified the distribution of expression of ZC3H12D in different crucial cell-types across datasets. The expression of ZC3H12D was collapsed by mean value. The gene expression level displayed using UMAP and violin plots was quantified by the normalized values.
Screening of co-expressed genes and enrichment analysis
We applied Spearman’s correlation analysis to identify genes that were related to the expression level of ZC3H12D. Those genes with correlation coefficient>0.5, P < 0.01 and FDR < 0.01 were considered as co-expressed genes for GO functional enrichment analysis and KEGG pathway enrichment analysis. The top 10 categories by GO analysis and top 10 pathways enriched by KEGG analysis were displayed, respectively. The volcano plot and heat map were drawn by LinkedOmics (http://www.linkedomics.org/login.php)  and Metascape (http://metascape.org/) , and the bubble diagram was drawn by ‘ggplot2’ R package (https://cran.r-project.org/web/packages/ggplot2/index.html). The sangerbox tool also provided assistance with image adjustments during drawing, a free online data analysis platform (http://www.sangerbox.com/tool).
The immune cell infiltration level and immune molecules expression
To obtain reliable immune infiltration estimations, we utilized the ‘immunedeconv’(https://www.rdocumentation.org/packages/immunedeconv/versions/2.0.0) R package that integrated six state-of-the-art algorithms , including Estimating the Proportion of Immune and Cancer cells (EPIC)  and Tumor Immune Estimation Resource (TIMER, http://timer.cistrome.org/) . Besides, immunostimulators and immunoinhibitors were selected and the expression levels (transcripts per kilobase of exonmodel per million mapped reads) of these genes were extracted. The two-gene, and the multi-gene correlation map was displayed using ‘ggstatsplot’ (https://cran.r-project.org/web/packages/ggstatsplot/index.html) and ‘pheatmap’ R package (https://cran.r-project.org/web/packages/pheatmap/index.html), respectively. Additionally, we used Spearman’s correlation analysis to describe the correlations among gene expression levels. And P < 0.05 was considered statistically significant.
Construction of ceRNA regulatory network in LUAD
Pearson correlation coefficient and significant p value between mRNAs and lncRNAs expression were calculated for differentially expressed miRNAs and their differentially expressed target mRNAs and lncRNAs, and for all known miRNAs and their differentially expressed target mRNAs and lncRNAs. In the present study, a total of 1955 differentially expressed mRNAs, 165 differentially expressed miRNAs and 1107 differentially expressed lncRNAs were identified by the RNA sequencing profiles of 10 stage I LUAD patients. The top 50 miRNAs and their target genes were ranked according to the respective network degree of both mRNAs as well as lncRNAs for display. Totally, 49 mRNAs, 99 miRNAs and 50 lncRNA were obtained to construct ceRNA regulatory networks (Supplement 1a). To identify a powerful pathway from this complicated ceRNA regulatory network, we focused on the function of lncRNA-miRNA-mRNA axises that could predict OS of LUAD. Verified by the Kaplan–Meier survival analysis, five mRNAs including ZC3H12D were identified to be significant associated with OS, including ZC3H12D, GNAO1, KSR2, SBK1 and SLIT3 (Fig. 1a-b and Supplement 2). Compared with TCGA normal data, LUAD samples had significantly higher expression of SBK1 and ZC3H12D, while lower SLIT3 expression was observed in LUAD. However, there was no significant difference in the expression of GNAO1 and KSR2 between tumor and paracancerous tissue samples. Meanwhile, both SBK1 and SLIT3 were found to play a vital role in lung cancer by miRNA or lncRNA [49,50,51]. Thereinto, we focused on ZC3H12D to construct a ceRNA regulatory network, and screened for potentially functional miRNAs and lncRNAs (Fig. 1c). Then, high expression of hsa-miR-4443 was not conducive to the prognosis of patients with LUAD (Fig. 1d), but high expression of ENST00000630242 was beneficial to the prognosis of patients with LUAD (Fig. 1e). Furthermore, ZC3H12D-hsa-miR-4443-ENST00000630242 axis was established, which might play a crucial role in the prognosis of LUAD patients.
Establishment of cox prognostic model
To explore the predictive role of ZC3H12D in LUAD patients, we constructed Cox regression models. Fortunately, ZC3H12D was beneficial to the overall survival of patients in both univariate and multivariate Cox regression models. (Fig. 2a and b). Based on the identifications from the multivariate cox proportional hazards analysis, a prediction nomogram was developed to calculate the risk by the points associated with the three risk factors (ZC3H12D, N stage, T stage). The total score ranged from 0 to 160, and was calculated by summing all the scores in each variable. According to the score calculated, we could roughly predict one-year survival rate, three-year survival rate and five-year survival rate of patients. A higher score indicated a higher risk of death (Fig. 2c). For the purpose of internal validation, the calibration plot showed that the nomogram predicting 1-year OS was relatively more accurate than the nomogram predicting 3-year OS and 5-year OS of LUAD patients (Fig. 2d).
Relationship between vital clinical factors and ZC3H12D expression level
As identified from the prognostic model, ZC3H12D, N stage and T stage could predict OS of LUAD patients. To explore the relationship between ZC3H12D expression and clinical factors associated with prognosis, we analyzed the significance P values by chi-square test. The results showed that T and N stages displayed significant differences between in group C1 and in group C2. In group C1, the expression of ZC3H12D was significantly higher than the median, while in group C2 was significantly lower than the median (Fig. 3a and b). From the correspondent bar chart, it exhibited that the proportion of ZC3H12D high expression in T1 and N0 was significantly higher than that of ZC3H12D low expression (Fig. 3c and d). To further verify the results, we combined N1, N2 and N3 group into one group, which was compared with N0 group to observe the difference of expression between the two groups. In N0 group, the ZC3H12D expression was significantly higher than another group (Wilcoxon rank sum test, P = 8.3e-05) (Fig. 3e). In addition, eliminating unknown samples, we took T2, T3 and T4 as a whole. Compared with T1, we found that there was a statistical significance between them and ZC3H1 2D expression was relatively higher in T1 phase (Fig. 3f).
ScRNA sequencing reveals ZC3H12D expression in normal lung tissue and LUAD
Different from conventional bulk RNA-sequencing, we applied scRNA-seq database to explore the expression level of ZC3H12D in normal lung tissue and LUAD. We found that ZC3H12D was mainly expressed in immune cells, which can be found in both human and mouse lung specimens (Fig. 4a-c). In parallel, we explored data from single-cell RNA sequencing of non-small cell lung cancer in public databases. Using scRNA-seq data from the GEO database, we found that ZC3H12D expression transformed by log (TPM/10 + 1) displayed heterogeneity in different clusters of cells in different NSCLC datasets (Fig. 4d). ZC3H12D was more abundantly expressed in CD4 T cells using GSE127465 and in Treg cells using GSE99254. Then, to eliminate the effect of lung squamous cell carcinoma we used another scRNA-seq dataset from GSE131907 to unravel the expression level in LUAD, indicating that it was relatively higher expressed at CD8 Tex(0.11), CD8 T cells(0.08), B cells(0.08) and CD4 T cells(0.07) over other cells(≤0.05) (Fig. 4e-h). Also, it was expressed in some plasma cells, DC cells, monocytes or macrophages, while ZC3H12D was not abundantly expressed in fibroblasts, mast and endothelial cells. However, when lung cancer metastasized to the brain, the expression levels of various types of cells changed. The ZC3H12D expression in monocytes was significantly increased, even higher than that in CD4 and CD8 T cells sometimes (Fig. 4i). And, the heterogeneity among individuals was also very significant. In another patient with leptomeningeal metastasis of lung cancer, the ZC3H12D expression was relatively higher in CD8 T cells and CD4 T cells (Fig. 4j).
Screening for important co-expressed genes
Base on Spearman correlation test, we identified 19,987 genes from the TCGA LUAD data. A total of 12,201 genes were positively correlated with ZC3H12D expression, while 7786 genes were negatively correlated with ZC3H12D expression, which displayed in the volcanic map (Fig. 5a). The top 10 positively correlated significant genes and the top 10 negatively correlated significant genes were screened and plotted, respectively (Fig. 5b and c).
GO functional enrichment analysis and KEGG pathway enrichment analysis
After further filtration of 12,201 positive related genes, 345 genes were obtained using correlation coefficient > 0.5, p-value < 0.01 and FDR < 0.01 as criteria for screening Gene ontology (GO) classification functional enrichment and Kyoto encyclopedia of genes and genomes (KEGG) classification pathway enrichment of 345 genes were performed. The results demonstrated that these genes were mainly enriched in the immune function of lymphocyte activation, immune response, B cell activation and alpha-beta T cell activation (Fig. 5d and Supplement 3). Meanwhile, co-expressed genes also participated in many pathways, such as cytokine-cytokine receptor interaction, cell adhesion molecules and Th17 cell differentiation (Fig. 5e).
Correlation between gene expression and immune cell infiltration
In the face of so much evidence showing the potential relationship between ZC3H12D and immunity, we decided to explore the correlation between expression of ZC3H12D and various immune cell infiltration levels in LUAD. Divided by the median, G1 group represents high expression of ZC3H12D, while G2 group represents low expression of ZC3H12D. Based on the experimental peritoneal cancer index (EPCI) algorithm, six major immune cells showed a trend: the scores of the G1 group with high expression were higher than the scores of the G2 group with low gene expression, interestingly, in the undefined cells, the scores of the G1 group with high gene expression were relatively lower than those of the G2 group with low gene expression (Fig. 6a).
Correlation between ZC3H12D and important immune molecules
Five hundred and ten LUAD samples were used for co-expression between ZC3H12D and other genes. Among the co-expressed genes screened, ZC3H12D expression was linked to many immune-related genes. To further explore its regulatory effect, we analyzed the correlation between the expression of ZC3H12D and the common immunostimulators. Ultimately, we identified 10 immunostimulators (CD27, CD28, CD40LG, CD48, CXCR4, ICOS, KLRK1, LTA, TNFRSF13B, TNFRSF13C) with statistical significance (Fig. 6b). Interestingly, we also found that the expression of ZC3H12D was related to immunoinhibitors. 10 statistically significant immunoinhibitors (BTLA, CD274, CD96, CTLA4, HAVCR2, LAG3, PDCD1, PDCD1LG2, TIGIT, SIGLEC15) were also screened, including several immune-checkpoint–relevant transcripts, such as CTLA4, CD96, TIGIT. However, the correlation between the immune checkpoint SIGLEC15 and ZC3H12D expression was relatively weak (Fig. 6c).
Screened from the ceRNA regulatory network of LUAD, ZC3H12D-hsa-miR-4443-ENST00000630242 pathway was found, which was closely related to the survival of patients. High expression of hsa-miR-4443 was not conducive to the prognosis of patients with LUAD, but high expression of ENST00000630242 was beneficial to the prognosis of patients with LUAD. Besides, high ZC3H12D expression was linked to various clinical characteristics and better prognosis. Enriched in immune cells, ZC3H12D was associated with various immune cell infiltration levels and immune molecules. The functional enrichment analysis also showed that the co-expressed genes mainly played a role in lymphocyte activation and cytokine-cytokine receptor interaction.
ZC3H12D is also called MCPIP4, C6orf95, and dJ281H8.1. It is a tumor suppressor gene, which plays a critical role in many cancers, including tongue cancer , osteosarcoma  and lung cancer. As for lung cancer, most of previous studies focused on the effect of genetic polymorphisms of ZC3H12D, and little was known about its potential regulatory mechanisms on lung cancer . Previous studies have demonstrated that ZC3H12D was associated with memory T lymphocytes and macrophages, participating in the regulation of inflammation [54, 55]. Therefore, we hypothesized that this gene might have an effect on LUAD through immune regulatory mechanism. In our study, we found that many immune-related genes were positively correlated with the expression of ZC3H12D, such as ZNF831, SLAMF1 and IL-16. Both ZNF831 and ZC3H12D were linked to Zinc Finger Family, and it has been reported that ZNF831 was specifically significant in the high immunity subtype of triple-negative breast cancer, which was characterized by anti-tumor immune activities, better immune cell infiltration and greater probability of OS . And SLAMF1 could both inhibit proliferation and impair responses to B cell receptor ligation in IGHV mutated chronic lymphocytic leukemia, which was similar to the anti-tumor effect of ZC3H12D . In the previous literature, ZC3H12D could regulate IL-6, which was a member of the interleukin family, and we further found that IL-16, another member of the interleukin family, was also closely related to ZC3H12D [12, 58]. As a pro-inflammatory cytokine, IL-16 was associated with high grade immune related adverse events in advanced NSCLC treated with immune checkpoint inhibitors .
In addition to the correlation of co-expressed genes, it is more convincing to analyze the function and pathway of positively related gene sets in LUAD samples. In our study, a lot of genes participated in the functions and pathways associated with immunity. The function enrichment results showed that many co-expressed genes were involved in lymphocyte activation and immune response, which were consistent with the results of previous studies [60, 61]. In addition to the immune-related pathways, we also found that the co-expressed genes of ZC3H12D expression were enriched in well-known cancer-related pathways, such as NF-kB signaling pathway, which needed to validate its function in cancer in the further.
Besides, the scRNA-seq technology provided us an innovative method to reveal the gene expression level in immune cells under different conditions [62, 63]. Based on the ScRNA-seq in our study, it showed that ZC3H12D was not homogeneous among different clusters in tumor, but it selectively highly expressed in immune cells. Compared with the samples of normal lung tissues, LUAD tissues, non-small cell lung cancer and brain metastasis of lung cancer, we found that the expression of ZC3H12D in different immune cell types, including conventional CD4 T cells, regulatory T cells, monocytes or macrophages, would change under different conditions, which reflected the plasticity of ZC3H12D to a certain extent. Also, when it was highly expressed, the EPIC scores of the main immune cells were correspondingly higher. Meanwhile, the high expression of ZC3H12D was also more common in T1 and N0, and was related to some immune molecules, so we speculated that it may play an anti-cancer effect by regulating immunity in the early stage of LUAD. Future studies are required to compare the immune changes caused by ZC3H12D between early-stage and advanced-stage, so as to further reveal the function of ZC3H12D on the dynamic heterogeneity of LUAD.
Apart from immune mechanism, ceRNA regulatory network may also involve in the prognosis of LUAD. In the present study, lncRNA ENST00000630242, acting as a ceRNA, could “sponge” hsa-miR-4443 to regulate the expression of target ZC3H12D. Overexpression of ZC3H12D was beneficial to prolong the survival time of LUAD patients, while hsa-miR-4443 was not conducive to the prognosis of patients, which were conformity with the pertinent literature [14, 64]. Although there were few studies on ENST00000630242, the data showed that it was beneficial to the prognosis of patients with LUAD. Therefore, ZC3H12D-hsa-miR-4443-ENST00000630242 axis could be served as a novel potential target for LUAD treatment.
The major limitation of this study is that we have not confirmed the ENST00000630242-hsa-miR-4443-ZC3H12D axis by experiments, though the correlations have been primarily uncovered through RNA-seq from clinical samples. Although the function of ZC3H12D was revealed in the present study, the mechanisms of ENST00000630242 and hsa-miR-4443 in LUAD are still unclear. Therefore, further study is needed to explore the role of ENST00000630242 and hsa-miR-4443 in LUAD. Additionally, due to the small number of samples performed by RNA-seq, it might induce bias, so we used the public database to validate our study. In the future, we will also make an effort to enlarge samples to reduce selection bias.
In summary, we found that ENST00000630242-hsa-miR-4443- ZC3H12D axis might be involved in the OS of LUAD patients. lncRNA ENST00000630242 could “sponge” hsa-miR-4443 to regulate the expression of ZC3H12D. ZC3H12D or ENST00000630242 could be beneficial to improve the prognosis of LUAD patients, while hsa-miR-4443 was not conducive to the overall survival time of patients. ZC3H12D, the core part of this pathway, was combined with some clinical factors to establish a cox model together. Furthermore, ZC3H12D expression at single level was unraveled in both normal lung tissues and lung tumors. Also, we found ZC3H12D expression was associated with some clinical features, important functions and pathways. Meanwhile, we explored the correlation between ZC3H12D and immune mechanisms to understand LUAD better. Therefore, it may serve as a vital predictive marker and could be regarded as a potential therapeutic target for LUAD in the future.
Availability of data and materials
The datasets generated and analysed during the current study are available, as follows: TCGA database (https://portal.gdc.cancer.gov/), GEO database (https://www.ncbi.nlm.nih.gov/geo/), Tumor Immune Estimation Resource (TIMER, http://timer.cistrome.org/), GO (http://www.bioconductor.org/packages/release/bioc/ html/topGO.html), KEGG (http://www.genome.jp/keggbin/show_organism?menu_type=pathway_maps&org=hsa), miRWalk 2.0 (http://zmf.umm.uni-heidelberg.de/ pps/zmf/mirwalk2/), StarBase (http://starbase.sysu.edu.cn/), KM-plotter (https://kmplot.com/ analysis/), Tabula Muris (https://tabula-muris.ds.czbiohub.org/), European Molecular Biology Laboratory (EMBL, https://www.ebi.ac.uk/), UCSC cell browser (https://cells.ucsc.edu/), TISCH (http://tisch.comp-genomics.org/), LinkedOmics (http://www.linkedomics.org/login.php), Metascape (http://metascape.org/) CNCI (https://github.com/www-bioinfo-org/CNCI), CPC2 (http://cpc2.cbi.pku.edu.cn/), CPAT (https://sourceforge.net/projects/rna-cpat/) and PLEX (https://sourceforge.net/projects/plek/). A dataset generated during the current study is not publicly available due to concerns regarding patient confidentiality and proprietary information but is available upon reasonable request from the corresponding author. We provided a reviewer link of unpublished BioProject. Use the following URL: https://dataview.ncbi.nlm.nih.gov/object/PRJNA732584?reviewer=jo5ph00rjrkr7u19thhft8eho.
Zinc finger CCCH-type containing 12D
Competing endogenous ribonucleic acid
Cytotoxic T-Lymphocyte Associated Protein 4
T-cell Immunoreceptor with Ig and ITIM domains
Messenger ribonucleic acid
Long non-coding RNA
Non-small cell lung cancer
Centromere proteins F
Single-cell RNA sequencing
Immediate early response 3
Tumor necrosis factor
Nuclear factor kappa-light-chain-enhancer of activated B cells
Family With Sequence Similarity 30 Member A
World health organization
Counts per million
Computerized progressive attentional training
Coding potential calculator 2
Predictor of long non-coding RNAs and mRNAs based on k-mer scheme
Differentially expressed mRNAs
Differentially expressed miRNAs
Differentially expressed lncRNAs
Kyoto encyclopedia of genes and genomes
The cancer genome atlas
European molecular biology laboratory
Gene expression omnibus
Red blood cell
Receiver operation characteristic
t-distributed stochastic neighbor embedding
Uniform Manifold Approximation and Projection for Dimension Reduction
University of California Santa Cruz
Tumor immune single-cell hub
False discovery rate
Tumor immune estimation resource
Guanine nucleotide-binding protein G(o) subunit alpha
Kachin special region II
SH3 domain binding kinase 1
Slit guidance ligand 3
Experimental peritoneal cancer index
C-X-C motif chemokine receptor 4
Inducible T cell costimulator
Killer cell lectin like receptor K1
Tumor necrosis factor receptor superfamily member
MCP Induced Protein 4
Chromosome 6 open reading frame 95
Zinc finger protein 831
Signaling lymphocytic activation molecule family member 1
Immunoglobulin heavy chain variable region
Janus kinase-signal transducer and activator of tran-ions
Bray F, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68(6):394–424.
Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin. 2019;69(1):7–34.
Guo L, et al. Construction and investigation of a combined hypoxia and stemness index lncRNA-associated ceRNA regulatory network in lung adenocarcinoma. BMC Med Genet. 2020;13(1):166.
Zhang Y, et al. The long noncoding RNA Linc01833 enhances lung adenocarcinoma progression via MiR-519e-3p/S100A4 axis. Cancer Manag Res. 2020;12:11157–67.
Zeng H, et al. Stemness related genes revealed by network analysis associated with tumor immune microenvironment and the clinical outcome in lung adenocarcinoma. Front Genet. 2020;11:549213.
He D, et al. Single-cell RNA sequencing reveals heterogeneous tumor and immune cell populations in early-stage lung adenocarcinomas harboring EGFR mutations. Oncogene. 2021;40(2):355–68.
Liu Y, et al. Single-cell transcriptome analysis demonstrates inter-patient and intra-tumor heterogeneity in primary and metastatic lung adenocarcinoma. Aging. 2020;12(21):21559–81.
Ruan H, et al. Circulating tumor cell characterization of lung cancer brain metastases in the cerebrospinal fluid through single-cell transcriptome analysis. Clin Transl Med. 2020;10(8):e246.
Guo L, et al. TOX correlates with prognosis, immune infiltration, and T cells exhaustion in lung adenocarcinoma. Cancer Med. 2020;9(18):6694–709.
Liang J, et al. A novel CCCH-zinc finger protein family regulates proinflammatory activation of macrophages. J Biol Chem. 2008;283(10):6337–46.
Wang M, et al. Identification of a novel tumor suppressor gene p34 on human chromosome 6q25.1. Cancer Res. 2007;67(1):93–9.
Wawro M, et al. Intact NYN/PIN-like domain is crucial for the degradation of inflammation-related transcripts by ZC3H12D. J Cell Biochem. 2017;118(3):487–98.
Zhu M, et al. miR-128-3p serves as an oncogenic microRNA in osteosarcoma cells by downregulating ZC3H12D. Oncol Lett. 2021;21(2):152.
Zhang W, Qiao B, Fan J. Overexpression of miR-4443 promotes the resistance of non-small cell lung cancer cells to epirubicin by targeting INPP4A and regulating the activation of JAK2/STAT3 pathway. Die Pharmazie. 2018;73(7):386–92.
Gong J, et al. lncRNA FEZF1-AS1 contributes to cell proliferation, migration and invasion by sponging miR-4443 in hepatocellular carcinoma. Mol Med Rep. 2018;18(6):5614–20.
Gao Y, et al. lncRNA MNX1-AS1 promotes glioblastoma progression through inhibition of miR-4443. Oncol Res. 2019;27(3):341–7.
Meerson A, Yehuda H. Leptin and insulin up-regulate miR-4443 to suppress NCOA1 and TRAF4, and decrease the invasiveness of human colon cancer cells. BMC Cancer. 2016;16(1):882.
Wu D, et al. Expression profiling and cell type classification analysis in periodontitis reveal dysregulation of multiple lncRNAs in plasma cells. Front Genet. 2020;11:382.
de Lima DS, et al. Long noncoding RNAs are involved in multiple immunological pathways in response to vaccination. Proc Natl Acad Sci U S A. 2019;116(34):17121–6.
Kang Y, et al. CPC2: a fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017;45(W1):W12-6.
Li G, et al. CPAT: Coding-Potential Assessment Tool using an alignment-free logistic regression model. Nucleic Acids Res. 2013;41(6):e74.
Li A, et al. PLEK: a tool for predicting long non-coding RNAs and messenger RNAs based on an improved k-mer scheme. BMC Bioinformatics. 2014;15(1):311.
Harsh Dweep, et al. miRWalk – Database: Prediction of possible miRNA binding sites by “walking” the genes of three genomes. J Biomed Inform. 2011;44(5):839–47.
Li J, et al. starBase v2.0: decoding miRNA-ceRNA miRNA-ncRNA and protein–RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 2014;42(D1):D92–7.
Der Maaten LV, Hinton GE. Visualizing data using t-SNE. J Mach Learn Res. 2008;9 (Nov):2579-605.
George C, et al. Fast interpolation-based t-SNE for improved visualization of single-cell RNA-seq data. Nature Methods. 2019;16(3):243–5.
Nagy Á, Munkácsy G, Győrffy B. Pancancer survival analysis of cancer hallmark genes. Sci Rep. 2021;11(1):6047.
Tabula MC, et al. Single-cell transcriptomics of 20 mouse organs creates a Tabula Muris. Nature. 2018;562(7727):367–72.
Lambrechts D, et al. Phenotype molding of stromal cells in the lung tumor microenvironment. Nat Med. 2018;24(8):1277–89.
Song Q, et al. Dissecting intratumoral myeloid cell plasticity by single cell RNA-seq. Cancer Med. 2019;8(6):3072–85.
Zilionis R, et al. Single-cell transcriptomics of human and mouse lung cancers reveals conserved myeloid populations across individuals and species. Immunity. 2019;50(5):1317–1334.e10.
Newman AM, et al. Determining cell type abundance and expression from bulk tissues with digital cytometry. Nat Biotechnol. 2019;37(7):773–82.
Kim N, et al. Single-cell RNA sequencing demonstrates the molecular and cellular reprogramming of metastatic lung adenocarcinoma. Nat Commun. 2020;11(1):2285.
Wu TD, et al. Peripheral T cell expansion predicts tumour infiltration and clinical response. Nature. 2020;579(7798):274–8.
Jeong SH, et al. Nomogram for predicting gastric cancer recurrence using biomarker gene expression. Eur J Surg Oncol. 2020;46(1):195–201.
McInnes L, Healy J, Melville J. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv preprint. 2018. arXiv:1802.03426
Li X, Li X, Ding L. Comprehensive Analysis to Identify Enhancer-Regulated Inflammation-Associated Genes in Lung Adenocarcinoma. Cancer Management and Research Volume. 2021;13:7115–29.
Jin X, et al. Microarray data analysis on gene and miRNA expression to identify biomarkers in non-small cell lung cancer. BMC Cancer. 2020;20:329.
Consortium TM. A single-cell transcriptomic atlas characterizes ageing tissues in the mouse. Nature. 2020;583(7817):590–5.
Travaglini KJ, et al. A molecular cell atlas of the human lung from single-cell RNA sequencing. Nature. 2020;587(7835):619–25.
Chi Y, et al. Cancer cells deploy lipocalin-2 to collect limiting iron in leptomeningeal metastasis. Science. 2020;369(6501):276–82.
Lee CM, et al. UCSC genome browser enters 20th year. Nucleic Acids Res. 2020;48(D1):D756–61.
Sun D, et al. TISCH: a comprehensive web resource enabling interactive single-cell transcriptome visualization of tumor microenvironment. Nucleic Acids Res. 2021;49:D1420–30.
Vasaikar SV, et al. LinkedOmics: analyzing multi-omics data within and across 32 cancer types. Nucleic Acids Res. 2018;46:D956–63.
Zhou Y, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10(1):1523.
Sun L, et al. Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 2013;41(17):e166–e166.
Friedländer MR, et al. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40(1):37–52.
Li T, et al. TIMER: a web server for comprehensive analysis of tumor-infiltrating immune cells. Cancer Res. 2017;77(21):e108–10.
Li Z, et al. LncRNA‐ENST00000501520 promotes the proliferation of malignant‐transformed BEAS‐2B cells induced with coal tar pitch mediated by target genes. Environ Toxicol. 2019;34(7):869–77.
Racle Julien, et al. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife. 2017;06. Available from: https://doi.org/10.7554/eLife.26476.
Etienne B, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biology. 2016;17(1):218.
Shen T, Wang M, Wang X. Identification of prognosis-related hub RNA binding proteins function through regulating metabolic processes in tongue cancer. J Cancer. 2021;12(8):2230–42.
Vizoso M, et al. Aberrant DNA methylation in non-small cell lung cancer-associated fibroblasts. Carcinogenesis. 2015;36(12):1453–63.
Huang S, et al. The putative tumor suppressor Zc3h12d modulates toll-like receptor signaling in macrophages. Cell Signal. 2012;24(2):569–76.
Emming S, et al. A molecular network regulating the proinflammatory phenotype of human memory T lymphocytes. Nat Immunol. 2020;21(4):388–99.
He Y, et al. Classification of triple-negative breast cancers based on Immunogenomic profiling. J Exp Clin Cancer Res. 2018;37(1):327.
von Wenserski L, et al. SLAMF receptors negatively regulate B cell receptor signaling in chronic lymphocytic leukemia via recruitment of prohibitin-2. Leukemia. 2021;35(4):1073–86.
Wawro M, et al. ZC3H12B/MCPIP2, a new active member of the ZC3H12 family. RNA. 2019;25(7):840–56.
Costantini A, et al. Plasma biomarkers screening by multiplex ELISA assay in patients with advanced non-small cell lung cancer treated with immune checkpoint inhibitors. Cancers. 2020;13(1):97.
Skundric DS, et al. Emerging role of IL-16 in cytokine-mediated regulation of multiple sclerosis. Cytokine. 2015;75(2):234–48.
Gordiienko I, et al. SLAMF1/CD150 in hematologic malignancies: silent marker or active player? Clin Immunol. 2019;204:14–22.
Mantsoki A, Devailly G, Joshi A. Gene expression variability in mammalian embryonic stem cells using single cell RNA-seq data. Comput Biol Chem. 2016;63:52–61.
Wang L, et al. Myeloid cell-associated resistance to PD-1/PD-L1 blockade in urothelial cancer revealed through bulk and single-cell RNA sequencing. Clin Cancer Res. 2021;27(15):4287–300.
Yang L, et al. Development and validation of a prediction model for lung adenocarcinoma based on RNA-binding protein. Ann Transl Med. 2021;9(6):474.
We acknowledge funding from funding project of Fujian Medical University College Student Innovation and entrepreneurship training program (Grant No. C19046) and funding project of Fujian Medical University College Student Innovation and entrepreneurship training program (Grant No. C19071). The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Ethics approval and consent to participate
The study was approved by the bioethical committees at The Second Affiliated Hospital of Fujian Medical University, China (2020-206).
Consent for publication
All participating patients provided written informed consent.
The authors declare that they have no competing interests in this section.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
a. A ceRNA regulatory network based on 49 mRNAs, 99 miRNAs and 50 lncRNAs. b. Association between ZC3H12D expression and LUAD patient survival time.
The survival curves of GNAO1, KSR2, SBK1 and SLIT3.
The genes were mainly enriched in the immune function of lymphocyte activation, immune response, B cell activation and alpha-beta T cell activation.
lncRNA and mRNA ceRNA network list top.
About this article
Cite this article
Chen, W., Guo, Z., Wu, J. et al. Identification of a ZC3H12D-regulated competing endogenous RNA network for prognosis of lung adenocarcinoma at single-cell level. BMC Cancer 22, 115 (2022). https://doi.org/10.1186/s12885-021-08992-1
- Zinc finger CCCH-type containing 12D (ZC3H12D)
- Competitive endogenous RNA
- Single-cell RNA sequencing
- Lung adenocarcinoma