In this study, six pairs of cancerous tissues and para-cancerous tissues were collected from male patients. The inclusion criteria for patients were as follows: male patients pathologically diagnosed with HCC, aged between 45 and 65 years, classified as Barcelona Clinic Liver Cancer (BCLC) Stage A [10], with no history of smoking or alcohol consumption, and who had not received preoperative radiotherapy or chemotherapy. The exclusion criteria included the presence of other malignant tumors or severe extrahepatic diseases. Two seasoned pathologists diagnosed these samples by relying on the rapid freezing pathology results obtained during the surgical procedure. A sample had a volume of 1 mm³ ($>$200 mg) to meet the criteria necessary for transcriptomics analysis and snRNA sequence. The para-cancerous tissue was collected by cutting at least 2 cm from the margin of the corresponding tumor. This study received approval from the Ethics Committee of Shandong University affiliated Shandong Provincial Third Hospital (No. KYLL-2024044; https://www.medicalresearch.org.cn/), and informed consent was acquired from the patients. They have been performed according to the Declaration of Helsinki.

The MHCC97-H cell line was purchased from Zhongqiaoxinzhou Company of China (ZQ0020, Shanghai, China) and cultured in dulbecco’s modified eagle medium (DMEM, ZQ-122, Zqxzh-Bio, Suzhou, China) at 37 °C in a 5% CO₂ incubator. All cell lines were validated by short tandem repeat (STR) profiling and tested negative for Mycoplasma (MycoAlert™ kit; LT07-318; Merck KGaA, Darmstadt, German). The primary TCs were carefully dissected from fresh liver tumor tissues under sterile conditions. Minimized the inclusion of fat and cut tissues into tiny pieces (~1 mm³) by surgical scissors, then digested them by dissolving collagenase (0.5 mg/mL, 2593923, BioMart, Shanghai, China), dispase II (2.4 U/mL, 40104ES60, BioMart, Shanghai, China), and DNase I (20 µg/mL, ALH044-DLW, Beijing Baiaoleibo Technology Co., Ltd., Beijing, China) in sterile PBS for one minute and filtered the suspension through a 55 µm cell strainer to remove undigested debris. Centrifuged the cell suspension at 300 $\times$g for 5 minutes to pellet the cells and cultured successfully by the complete medium (ZQ-121, Zhongqiaoxinzhou Company, Suzhou, China) for 48 hours. Collect the non-adherent cells which are more enriched in TCs into a new dish and incubated in a humidified incubator with 5% CO₂ at 37 °C. Then, CD34 (No130-097-418, Guiechem, Xian, China) and CD31 immunomagnetic beads (130-097-047, Nuowei company, Beijing, China) were used to separate the TCs from perithelial and epithelial cells, and TCs were selected successfully [7, 11, 12, 13]. The identification of TCs was performed by transmission electron microscope assay and immunofluorescent staining with CD34, PDGFRA, and CD31 according to our previous studies [7, 14, 15]. The primary TCs were validated by STR profiling and tested negative for Mycoplasma. Moreover, these two cell lines were rechecked out of “Error cells directory” (http://www.cellresource.cn; Procell Life Science & Technology Co., Ltd.). Modified version: Experimental groups were established by co-culturing TCs (1 $\times$ 10⁴ cells) alongside HCC cell lines (3 $\times$ 10⁵ cells) within double-deck petri dishes. In contrast, the control group involved only TCs that were cultured for 48 hours at 37 °C in a 5% CO₂ incubator using the same ZQ-121 medium. Following this, nucleic acid data from both TCs groups were analyzed through high-throughput gene sequencing technology, and the resulting transcriptomic data were employed for further analysis.

Samples undergo treatment with a particular enzyme to produce a suspension of cell nuclei. This cell suspension goes through several processing steps: it is initially filtered with a 60 µm filter screen, and then a second filtration is performed using a 40 µm filter screen to improve the removal of cell clusters. Once the reverse transcription step is completed, the emulsions are disrupted, and the barcoded cDNA is extracted, followed by RNA amplification. The amplified cDNA obtained is employed for the creation of a 3^′ gene expression library. Within the construction of the library, 50 ng of the amplified cDNA is fragmented, end-repaired, subjected to size selection using SPRIselect beads, and then sequenced on the NovaSeq platform (Berry Genomics Company, Beijing, China) 16 to yield 150 bp paired-end reads. The full-length cDNA is subsequently amplified via PCR to ensure an adequate quantity for library preparation. During the incubation of gel bead in emulsion (GEMs), TruSeq Read 1 (a specific primer sequence for Read 1) is incorporated into the molecules. The following processes of End Repair, A-tailing, Adaptor Ligation, and PCR incorporate P5, P7, a sample index, and TruSeq Read 2 (a specific primer sequence for Read 2) into the mixture. Quality control is performed using Qubit HS quantification and Agilent 2100 Bioanalyzer/Fragment Analyzer 5300, resulting in a final library size of approximately 450 bp. After sequencing and quality assessment, raw reads are demultiplexed and aligned to the reference genome via the 10X Genomics Cell Ranger pipeline [16], utilizing default settings. The Seurat package was used [17, 18]. The CellMarker website (http://xteam.xbio.top/CellMarker/) was utilized to guide the classification of cell-types [19]. All raw data was produced by Berry Genomics Corporation (http://www.berrygenomics.com).

After a 48-hour co-culture with HCC cell lines, the TCs cells and those from the control group were gathered. Total RNA extraction from both groups was performed utilizing the R1200 kit from Solarbio, located in Shanghai, China. The RNA underwent qualification and quantification through the following procedures: (1) RNA purity and concentration assessments were carried out using a NanoDrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA); (2) RNA integrity and amount were evaluated with the Agilent 2100/4200 system (Agilent Technologies, Santa Clara, CA, USA). Following the construction of the library, its concentration was assessed using the Qubit® fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). The precise concentration of the cDNA library was further evaluated through qPCR. The library’s size distribution was examined using agarose gel electrophoresis. Initial raw data in FASTQ format was processed with custom Perl scripts. Feature-Count was applied to quantify the reads. Differential expression analysis was performed using EdgeR (Walter and Eliza Hall Institute of Medical Research, Parkville, VIC, Australia), with p-values adjusted via the Benjamini and Hochberg method to control for the false discovery rate. Genes exhibiting $|$log2 (FoldChange)$|$ $>$1 and q value 0.05 were classified as differentially expressed [20, 21].

The HCC tissue sample is collected and immediately fixed using a chemical solution (H910607, Haisen Biology company, Hefei, China), to preserve its cellular structures and prevent degradation. Following fixation (Glutaraldehyde, H910607, Hisen company, Shanghai, China), the sample undergoes a series of dehydration steps in graded concentrations of ethanol or acetone to remove all water content. The dehydrated sample is then subjected to critical point drying, a technique that replaces the dehydrant with a liquid that can be transitioned to gas without surface tension effects, preserving the tissue’s structure. Next, the dried sample is mounted on a suitable stub to make it conductive and improve electron signal quality. Finally, the prepared sample is placed inside the scanning electron microscope (SEM) (HITACHI, SU9000, Sinaida, Shanghai, China) chamber, and the microscope is adjusted to focus on the area of interest. High-energy electrons are scanned over the surface, and the secondary or backscattered electrons emitted from the sample surface are collected to form a detailed image of tissue’s surface morphology, revealing the complex cellular and extracellular structures at high magnification.

The Cancer Genome Atlas (TCGA)-LIHC database (https://portal.gdc.cancer.gov/projects/TCGA-LIHC), including 424 RNA-seq samples, was used to assess the prognostic values of overall survival (OS), different expression levels with TNM stages, and receiver operating characteristic (ROC) curves of which differential expression genes (DEGs) were screened from snRNA-seq assay. Xiantaoxueshu website was a data visualization software and purchased from Jieluoxuan company of China (https://www.xiantaozi.com/) which could construct heatmap, OS curves, ROC curves, column diagrams and chord diagrams of relationships between DEGs in our study. We utilized the Gene Ontology (GO) database (http://geneontology.org) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (https://www.genome.jp/kegg/) [22, 23], accessed via Bioconductor (http://www.bioconductor.org/) to explore the signaling pathways associated with the DEGs. For identifying protein-protein interactions (PPI), we used the STRING database (http://string-db.org/), applying an interaction score threshold of $>$0.4. A confidence score of $\ge$700 was set to determine significance. After downloading PPI data, Cytoscape software (3.9.1 version, https://cytoscape.org) with the CytoHubba plug-in was then used to search and visualize the hub genes with the highest degree of connectivity via MCC algorithms. Correlation analysis, Venn DEGs screening, and heatmap of transcriptomics and snRNA-seq results were screened by limma R package (fgsea and cluster-Profiler plug-in; Walter and Eliza Hall Institute of Medical Research, Parkville, VIC, Australia; https://bioconductor.org/packages/release/bioc/html/limma.html) [24, 25], and infiltrated DEGs were next uploaded into Gene Set Enrichment Analysis (GSEA) system to verify enriched pathways shown in our study. The iTALK software (https://github.com/Coolgenome/iTALK) is utilized in the analysis of single-nucleus transcriptomes to evaluate the variances in gene expression of ligands and receptors across different cell types by comparing various samples. Meanwhile, PySCENIC (version 0.12.1; https://scenic.aertslab.org/) identifies modules of co-expression among transcription factors and target genes. It constructs gene mediated networks specific to distinct cell subsets and uncovers critical regulators that uphold cellular characteristics through the area under the curve (AUC)-cell algorithm, which incorporates adjusted values of the Regulon Specificity Score (RSS) [26].

Comparisons among several groups were performed utilizing one-way analysis of variance (ANOVA), and subsequent post hoc analyses were carried out using Tukey’s method. The findings are presented as means $\pm$ standard deviation (SD), with a significance threshold established at p 0.05. Additionally, univariate and multivariate analyses were performed, including Cox proportional hazards regression, to assess the impact of various clinical and pathological factors on overall survival. We also employed ROC curve analysis to determine the diagnostic performance of biomarkers and identify optimal cut-off values. Overall survival was calculated using the Kaplan-Meier method, with differences in survival curves analyzed through the log-rank test.

The snRNA-seq technique was employed to create an extensive representation of the cell clusters present in the tumor microenvironment (Fig. 1). The Loupe browser software, along with Cell Ranger data, was then employed to define the machine clustering outcomes, which indicated the identification of 13 cell clusters within the cancer tissues and 20 clusters in the para-cancerous tissues. These clusters were additionally classified into seven unique subsets for both HCC and the para-cancerous tissues: activated T cells, cancer cells, dendritic cells, Kupfer cells, TCs, vascular endothelial cells, and hepatocytes (Fig. 2A,B,H,I; Supplementary Fig. 1A). Remarkably, TCs were identified by employing 25 biomarkers derived from prior study [27] (Supplementary Table 1), and the mean frequency of TCs were initially quantified as 0.35% in the HCC tissues (Fig. 2C). The proportion of TCs and Kupfer cells was discerned to be lower in the cancer tissues compared with the para-cancerous tissues (0.35% vs 8.19%; 13.64 vs 18.89%) (Fig. 2C,J). Subsequently, the diverse heatmaps of DEGs in sub-cell populations underscored their heterogeneity in HCC (Fig. 2D–G,K–N; Supplementary Tables 2,3,4), and Uniform Manifold Approximation and Projection (UMAP) diagrams along with violin plots of the top 14 DEGs were created (Supplementary Fig. 1B,C). To elucidate the functional implications of DEGs, gene ontology and pathway enrichment analyses were conducted using the GO/KEGG browser (Supplementary Fig. 2). These analyses revealed that the DEGs were implicated in several signaling pathways including Rap1 (Rap1 orthologs from arthropods), PI3K/AKT, HIF-1 (Hypoxia inducible factor 1), ECM (Multimerin), MAPK, VEGF, and JAK/STAT. These pathways are instrumental in regulating cellular processes such as proliferation, migration, morphogenesis, and angiogenesis. The DEGs in tumor cells were linked to the cell cycle, resistance to platinum-based drugs, the p53 signaling pathway, and the Wnt signaling pathway. These pathways play roles in chromosome segregation, the regulation of the mitotic spindle, and the activity of ATP-dependent microtubule motors.

Fig. 1.

The flow chart of materials and methods. SnRNA-seq, single nucleus RNA sequencing; HCC, hepatocellular carcinoma; TCs, telocytes; MH97, MHCC97 cells; GO/KEGG, Gene Ontology/Kyoto Encyclopedia of Genes and Genomes; PPI, protein-protein interaction; ROC, receiver operation characteristic curve; OS, overall survival.

Fig. 2.

Analyses of single nucleus RNA sequencing of HCC tissues. (A,B) The UMAP and t-SNE graphs showed sub-cell clusters in HCC tissues. (C) The proportions of sub-cells in HCC tissues. (D,E) The heatmap of DEGs in different sub-cells (Top 35) with Z-scores. (F) The volcano plot of DEGs was screened. (G) The bubble diagram of parts of DEGs with average levels. (H,I) The UMAP and t-SEN graphs showed sub-cell clusters in para-cancerous tissues. (J) The proportions of sub-cells in para-cancerous tissues. (K,L) The heatmap of DEGs in different sub-cells (Top 35) with Z-scores. (M) The volcano plot of DEGs was screened. (N) The bubble diagram of parts of DEGs with average levels. DEGs, differential expressed genes.

To investigate the intricate interactions between TATCs, immune cells and cancer cells in HCC, we showcased their cluster UMAP graphs for them using the Loupe browser. These graphs indicated that activated T cells were predominant in terms of quantity (Fig. 3A). A comprehensive GO/KEGG analysis showed that the TATCs were closely related with CD8⁺ T cell differentiation, but no clustering information on dendritic cells’ behavior with TATCs was obtained (Fig. 3B). Moreover, TOP2A (DNA topoisomerase Ⅱ alpha), BUB1B (BUB1 mitotic checkpoint serine/threonine kinase B), KIF11 (kinesin family member 11), and CENPF (centromere protein F) were identified as hub genes of TCs through the use of the String website and Cytoscape software (Spearman’s coefficient $>$0.8; p 0.05), while CD74, the PRKACB (protein kinase cAMP-activated catalytic subunit beta) of T cells, IRF8 (interferon regulatory factor 8), and the SPI1 (Spi-1 proto-oncogene) of dendritic cells were also highlighted (Fig. 3C–E). The ring heatmap depicted the gene expressions in divergent cell clusters, confirming that TATCs possess unique heterogeneity in their molecular profile compared with immune cells (Fig. 3F). To delve deeper into the cell-to-cell interaction signals that are mediated by these DEGs, we utilized iTALK and PySCENIC software (University of California, San Francisco, CA, USA), renowned for their proficiency in bioinformatics analysis and visualization [28]. Notably, the scatter diagram of regulon activity demonstrated that E2F8 (E2F transcription factor 8) in TATCs achieved a high Z-score, signifying a crucial role in cell proliferation (Fig. 3G; Supplementary Table 5). The visualization analysis demonstrated that the ITGB1 in TCs acts as a crucial receptor for PLG (plasminogen), FN1 (fibronectin 1), and RELN (reelin) within hepatocytes, cancer cells, and vascular endothelial cells. Furthermore, the receptor-ligand interactions revealed molecular communication among cancer cells, vascular endothelial cells, dendritic cells, and TATCs (Fig. 3H). RRM2 (ribonucleotide reductase regulatory subunit M2) and KNL1 (kinetochore scaffold 1) emerged as the top two genes in TATCs with Log2FC $>$2, and the lollipop graph implied a negative correlation of TATCs with dendritic cells, neutrophils, CD8⁺ T cells, and NK cells, while showing a positive correlation with Th2 cells and T helper cells (Fig. 3I,J). Collectively, the snRNA-seq technique facilitates the distinction of TCs from other cells in the tumor microenvironment (TME) and delineates their cross-talk with T cells and dendritic cells via regulons and hub genes.

Fig. 3.

Bioinformatic analysis of TATCs and immune cells in HCC. (A) The UMAP of T cells, cancer cell, dendritic cell and TCs. (B) GO/KEGG analyses of dendritic cell, activated T cell, and TCs (https://www.kegg.jp/kegg/kegg1.html). (C–E) The protein-protein interactions and hub gene of dendritic cell, activated T cell, and TCs by Cytoscape software. (F) The ring heatmap of DEGs in dendritic cell, activated T cell, Kupfer cell and TCs. (G) The bubble graph of regulons in seven sub-cell clusters with Z-scores. (H) The correlations among sub-cell clusters according to snRNA-seq data by iTALK and PySCENIC software. (I,J) The correlation between TCs and immune cells with top two expression in TCs. p 0.05 as a significant value. *p 0.05, **p 0.01, ***p 0.001, ns $>$ 0.05.

The scanning electron microscope assay, acknowledged as the optimal method for visualizing TCs in cancer tissues [14], enables the identification of the morphology of TCs and HCC cells, thereby substantiating their associative relationship (Fig. 4A). To mimic the presence of TCs in HCC, we executed in vitro experiments employing TCs and the HCC cell line. Following this, we examined the transcriptomics (mRNA) of the TCs and TATCs, uncovering that 622 genes were upregulated while 758 were downregulated in TATCs contrasted to TCs (Fig. 4B; Supplementary Table 6). Furthermore, these DEGs were subjected to a GO/KEGG network analysis, demonstrating their involvement in the mechanisms of channel activity, transmembrane transporter activity, and the sensory perception of the channel complex (Fig. 4C; Supplementary Table 7). This finding illuminates the direct influence of HCC cancer cells on the genetic modulation of TCs.

Fig. 4.

Transcriptomic changes and clinical significance of TCs after co-cultured with the MH97 cell line. (A) Graphs of HCC cancer cells and TCs under the scanning electron microscope. Scale bars = 10 µm, and 20 µm. (B) The heatmap of DEGs between TCs and TCs co-cultured with MHCC97-H cells. (C) Network chart of DEGs by the GO/KEGG enrichment. (D) The Venn graph of TCs in HCC tissues and TCs co-cultured MH97 cells. (E) Nine quadrants graph of TCs’ DEGs from transcriptomic testing and snRNA-seq technology with six co-expression genes. (F) The correlation of six co-expression genes. (G) Network chart of co-expression genes. (H) Expression levels of co-expression genes in HCC tissues and adjacent cancer tissues by TCGA database. (I–N) ROC curves and OS curves of co-expression genes in HCC. **p 0.01, ***p 0.001.

To identify biomarkers of TATCs in HCC, we combined snRNA-seq with transcriptomics data, focusing on the analysis of shared genes. The Venn diagram displayed six co-clustered genes derived from the raw sequencing data: TMC5 (transmembrane channel-like 5), KCNQ5 (potassium voltage-gated channel subfamily Q member 5), ZNF385D (zinc finger protein 385D), SNAP91 (synaptosome-associated protein 91), CRP (C-reactive protein), and SLC35F3 (solute carrier family 35 member F3) (Fig. 4D,E). It is worth noting that TMC5 showed a strong correlation with CRP, which is potentially involved in opsonin binding, phosphatidylinositol binding, and regulation of macrophage derived foam cell differentiation (Fig. 4F,G). Further analysis revealed that five out of these six biomarkers demonstrated significant alterations between HCC tissues and adjacent cancer tissues according to the TCGA-LIHC database: CRP was expressed at lower levels in tumors, whereas TMC5 and ZNF385D exhibited higher expression levels in tumors, while KCNQ5, SLC35F3, and SNAP91 had no divergency in either tumors or adjacent cancer tissues based on the TCGA database (Fig. 4H). However, these significant biomarkers might be deemed insufficient as diagnostic predictors due to their inadequate an AUC values (ZNF385D = 0.762; TMC5 = 0.647; Fig. 4I). Furthermore, the OS curves showed that high levels of SLC35F3 and TMC5 in HCC patients were associated with a poor prognosis, while CRP and ZNF385D had no statistical significance based on the TCGA database (Fig. 4J–N. The univariate Cox regression analyses of TMC5, CRP, SLC35F3, and ZNF385D based on detailed information showed that high TMC5 expression is significantly associated with a poor prognosis (HR = 1.499; 95% CI: 1.059–2.122; p = 0.023), as is high SLC35F3 expression (HR = 1.562; 95% CI: 1.102–2.214; p = 0.012) (Table 1). Taking together, although TMC5 and SLC35F3 were corelated with OS and considered markers of TCs, TMC5 is more suitable as a biomarker due to its bigger AUC of the ROC curve.