Huh-7 (KCLB60104; Korean Cell Line Bank, Seoul, Korea), was cultured at 37℃ with 5% CO₂ in Dulbecco’s modified eagle medium (Gibco; Thermo Fisher Scientific, Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS) (Gibco; Thermo Fisher Scientific), 4.5 g/L glucose, L-glutamine, and 1% penicillin/streptomycin.

Cells were harvested from Huh7 cell lines, and suspended in phosphate-buffered saline (PBS) supplemented with 2% FBS at a concentration of 1×10⁷ cells/mL. The cells were then incubated with an anti-CD133 (W6B3C1; BioLegend, San Diego, CA, USA) antibody conjugated to fluorescein isothiocyanate (FITC) at the dilution recommended by the manufacturer for 30 min on ice in the dark. Following incubation, the cells were washed twice with PBS and 2% FBS to remove excess antibodies, and centrifuged at 300× g for 5 minutes. The resulting pellets were carefully resuspended in a small volume of PBS containing 2% FBS. The MACS system (Miltenyi Biotec, Bergisch Gladbach, Germany) was used according to the manufacturer’s protocol as described.13 A fluorescence-activated cell sorter (FACS) equipped with appropriate lasers and detectors was utilized for the selected fluorophores. The FACS was set to identify and sort CD133+ (positive) and CD133- (negative) populations based on fluorescence intensity. CD133+ cells were collected in one tube and CD133- cells in another. To confirm the purity of isolated populations, additional staining of sorted cells with anti-CD133 antibody was performed following FACS analysis. The percentage of CD133+ and CD133- cells in both sorted populations was assessed, with an expected purity greater than 90%. FACS analysis software (e.g., FlowJo) was used to determine the percentage and absolute count of CD133+ and CD133- cells in the populations (data not shown).

For RT-qPCR studies, cells were seeded in 6-well plates, allowed to adhere until reaching 85-90% confluence, and then incubated overnight in serum-containing conditions. The mRNA levels of human TERT, HKR3, and cell cycle factors were determined using SYBR green-based semi-quantitative PCR. RNA was extracted using a QIAGEN RNeasy Mini Kit (QIAGEN, Hilden, Germany) according to the manufacturer’s instructions and reverse transcribed (RT) using a Clontech RT Kit (Clontech, Mountain View, CA, USA). Following reverse transcription, RT-qPCR analysis was performed using Applied Biosystems (Thermo Fisher Scientific) SYBR Green Master Mix and specific PCR primers (Supplementary Table 1). Amplification efficiencies were calculated for all primers using serial dilutions of the pooled cDNA samples. Data were analyzed using the comparative ΔΔCt method and normalized to beta-actin, the housekeeping gene. Results are presented as mean± standard error of the mean (SEM) (error bars) from at least three independent experiments and represented as fold-change relative to controls. Melting curves were generated to ensure the presence of a single gene-specific peak. No-template controls were included for each run, and each set of primers was used as a control for non-specific amplification.

Total RNA was extracted using TRIzol^® RNA Isolation Reagents (Life Technologies, Carlsbad, CA, USA). The quantity and quality of total RNA were evaluated using an Agilent 2100 Bioanalyzer RNA Kit (Agilent, Santa Clara, CA, USA). The isolated total RNA was processed to prepare an mRNA-sequencing library using the Illumina TruSeq Stranded mRNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s protocol. Library quality and fragment size were assessed using an Agilent 2100 Bioanalyzer DNA Kit (Agilent). All libraries were quantified by qPCR using the CFX96 Real-Time System (Bio-Rad, Hercules, CA, USA) and sequenced on a NextSeq500 sequencer (Illumina) with a paired-end 75 bp read and a single 8 bp index read run.

Potential sequencing adapters and law-quality bases in the raw reads were trimmed using Skewer software. The options -x AGATCGGAAGAGCACACGTCTGAACTCCAGTCA and -y AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT were applied for removing common TruSeq adapter sequences (Illumina). Additionally, the parameters -q 0 -l 25 -k 3 -r 0.1 -d 0.1 were used to trim low-quality bases from the 5′ and 3′ ends of the raw reads. After trimming, the cleaned high-quality reads were mapped to the reference genome using STAR software. Because the sequencing libraries were prepared in a strand-specific manner using strand-specific library preparation kit (Illumina), the strand-specific library option, library-type=fr-firststrand, was applied during the mapping process.

To quantify the mapped reads on the reference genome for gene expression analysis, Cufflinks with the strand-specific library option (library-type=fr-firststrand), and other default settings were used. The hg19 gene annotation from the UCSC genome browser (https://genome.ucsc.edu) in GTF format was used as the gene model, and expression values were calculated in fragments per kilobase of transcript per million fragments mapped units. The differentially expressed genes between the two selected biological conditions were analyzed using Cuffdiff software in the Cufflinks package5 with the strand-specific library option, (library-type=fr-firststrand) and other default settings. To compare the expression profiles among samples, the normalized expression values of a selected subset of a few hundred differentially expressed genes were clustered using unsupervised hierarchical clustering in in-house R scripts. Additionally, scatter plots for gene expression values and volcano plots for expression-fold changes and P-values were generated using in-house R scripts.

To gain insight into the biological functions associated with the differential gene expression between the compared biological conditions, a gene set overlapping test was performed between the analyzed differentially expressed genes and functionally categorized genes, including biological processes of gene ontology, KEGG pathways, and other functional gene sets using g:Profiler.

For the proteomic analysis, each group of samples was cultivated in triplicate under the conditions described above. The samples were dissolved using RIPA buffer. The protein concentration was determined by bicinchoninic acid (BCA) assay, and the samples were stored at -70℃ for further study. A total of 15 µg of protein samples was separated using 12% sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS-PAGE). The gels were fractionated into seven sections based on their molecular weights. In-gel digestion was performed as previously described.14 The resulting tryptic peptides were dissolved in 0.5% trifluoroacetic acid for liquid chromatography with tandem mass spectrometry (LC-MS/MS) analysis.

Peptides were separated using an Ultimate 3000 UPLC system (Dionex, Sunnyvale, CA, USA) connected to a Q-Exactive Plus mass spectrometer (Thermo Fisher Scientific) equipped with a nanoelectrospray ion source (Dionex). Peptides were first concentrated on a trapping column (75 μm inner diameter) packed with 5 μm C18 particles (Acclaim PepMap; Thermo Fisher Scientific) and then separated using a 15 cm analytical column packed with RSLC C18 reversed-phase column (Thermo Fisher Scientific) at a flow rate of 300 nL/min. The peptides were eluted using a 0-95% acetonitrile gradient in 0.1% formic acid over 120 minutes. All MS and MS/MS spectra were acquired in the data-dependent top-10 mode, with automatic switching between full-scan MS and MS/MS acquisition. Survey full-scan MS spectra (m/z 150-2,000) were acquired using the Orbitrap at a resolution of 70,000 (m/z 200) after ion accumulation to a 1×106 target value based on predictive automatic gain control (AGC) from the previous full scan. The MS/MS analysis was performed three times for each sample. The raw MS/MS spectral files were converted to mgf files using the Proteome Discoverer Daemon ver.1.4 (Thermo Fisher Scientific). The converted mgf files were used for protein identification using MASCOT ver.2.4.1 (Matrix Science, Boston, MA, USA; www.matrixscience.com). Protein quantification was performed by calculating the exponentially modified protein abundance index (emPAI)15 and the mol %, a percentile normalized value, based on the emPAI values. The search parameters applied in the database searches were as follows: enzyme specificity-trypsin/P, maximum missed cleavage-1, carboxymethyl as a static modification, oxidation and N-terminal acetylation as dynamic modifications, precursor mass tolerance of 10 ppm, and MS/MS mass tolerance of 0.8 Da. The MS/MS data were filtered according to a false discovery rate (FDR) of <1% to calculate the ratio of false positive matches in the decoy database based on the number of matches in the database. The UniProt human database (2018.05, www.uniprot.org), which has a total 93,614 sequences and 37,039,836 residues, was combined with cRAP database (2009.05.01, ftp://ftp.thegpm.org/fasta/cRAP/crap.fasta) which includes 116 sequences and 38,459 residues, and the contaminants database (2016.01.30, http://www.matrixscience.com/downloads/contaminants.zip), which consists of 247 sequences and 128,130 residues. These databases of common contaminants were used for protein identification. Relative quantitation and significance were provided for all the identified variables. Raw data are available via ProteomeXchange under the identifier PXD#331223 (https://www.ebi.ac.uk/pride/archive/).

Functional analysis of the dataset was performed using protein network analysis (https://string-db.org/). The functional analysis identified the most significant biological functions and/or diseases associated with the dataset. The analysis was performed against the string test to calculate a P-value, determining the probability that each assigned biological function and/or disease was due to chance alone. The focus genes from the two shortlists were overlaid onto a global molecular network developed using the information contained in the STRING database. This generated network based on the connectivity of individual proteins. Simultaneously, protein clusters and assignments of differentially expressed genes were generated. To gain insight into the biological functional role of the differential gene expression between the compared biological conditions, a gene set overlapping test was performed. This test compared the differentially expressed genes with functionally categorized gene sets, including biological processes of Gene Ontology, KEGG pathways, and other functional gene sets using g:Profiler.

Results are expressed as means±SEM or frequency (%). An independent t-test was performed to compare the differences in means between the control and experimental groups. All statistical analyses were performed using SPSS version 12.0 (IMB, Armonk, NY, USA). Statistical significance was set at P<0.05.

To analyze the propensity of CSCs specific to liver cancer, cell lines with CSC characteristics were isolated from HCC cell lines and confirmed to exhibit CSC properties. Based on CD133, a representative CSC marker, a CD133+ cell line was isolated from the HCC cell line Huh7, and a CD133- cell line lacking CSC characteristics was also isolated. To confirm whether the isolated CD133+ Huh7 cell line exhibited CSC propensity, the expression of the representative CSC markers CD44, CD24, epithelial cell adhesion molecule (EpCAM), and CD90 was analyzed (Fig. 1). Compared to wild-type Huh7 cells before isolation, EpCAM, CD44, and CD133 showed significantly higher expression patterns in the CD133+ Huh7 cell line, and although not significant, CD24 also showed high expression patterns. In contrast, CD90, which is recognized as a novel CSC marker in HCC, showed high expression overall but lower expression, although not significant, compared to the wild-type Huh7 cells. However, this was still higher than the CD133- cell line.

The proteomic profiling analysis of CD133+ Huh7, which expresses four out of five CSC markers, was compared with that of wild-type and CD133- cells. CD133+ HCC cell lines exhibited distinct expression patterns compared to wild-type or CD133-cells, while wild-type and CD133- cells showed similar expression patterns (Fig. 2A). In CD133+ cells, the expression of genes corresponding to biological regulation was the highest, accounting for approximately 60.6%. Genes related to cellular components and metabolic processes were expressed at approximately 18.0% and 17.9%, respectively. Additionally, genes related to the immune response and rhythmic processes were identified with expression levels of 6.9% (Fig. 2B). In the metabolic process group, the amine metabolic process showed the highest expression at approximately 39%, followed by the alcohol metabolic process, a key function of hepatocytes, at approximately 30%. In the biological regulation category, various expression patterns were observed, ranging from cell proliferation to response to stimulation. Notably, the highest expression difference was approximately 67% in the biological regulation part of the cell. Among the cellular components, which showed the third highest expression pattern, a high expression level of approximately 57% was observed in the structural-developmental part (Fig. 2C).

Analysis of the enrichment protein patterns among the biological regulation category, which showed high expression changes, revealed that the expression z-score of calcium ion binding-related proteins was higher than that of the CD133- type (Fig. 3A; Table 1). Among the 19 calcium ion binding proteins, 13 were observed in the functional annotation enrichment analysis. Correlation analysis of these 13 proteins confirmed that 10 of them have a close relationship with each other on the gene-to-gene network (https://string-db.org/) (Fig. 3B). Nearly all proteins forming these 13 networks are associated with tumor proliferation, metastasis, and suppression of the immune response and are closely linked with CSC survival, invasion, metastasis, and resistance to therapeutic drugs. In particular, the S100 family members MYL12B and TPT1, which are highly correlated with CD133+ status, were found to be closely associated with CSC infiltration and multidisciplinary resistance.

To confirm the effect of proteomic profiling on actual expression, the expression patterns of ten closely related calcium ionbinding proteins in the CD133+ and CD133- groups were compared relative to the wild-type cells (Fig. 4). Among the calcium-binding proteins, the S100 family was highly expressed in CD133+ cells, with MYL12B and TPT1 also showing significantly higher expression (P<0.05, P<0.01). Although not significant, LCP1, PLS3, and PDCD6 also exhibited higher expression levels compared to the wild-type cells, displaying an expression pattern opposite to that of CD133-. Conversely, ANXA1 expression was significantly suppressed in CD133+ cells compared to both wild-type and CD133- cells (P<0.01), while ACTN1 did not show a significant difference among the three groups.