The expression arrays and associated clinical information of patients with HCC from The Cancer Genome Atlas (TCGA) were obtained via the University of California Santa Cruz Xena tool (http://xena.ucsc.edu/) (18) for subsequent analysis. After excluding samples with no clear clinical information and OS <30 days, a total of 362 patients participated. The Human Protein Atlas portal (https://www.proteinatlas.org/) was searched to acquire images of HCC and adjacent normal tissue immunohistochemistry (IHC) to show dysregulation of a key gene at the protein level. The gene set for FA metabolism was obtained from the Molecular Signatures Database (MSigDB) (19,20). After entering ‘fatty’ in the keywords box with the filter of Homo sapiens, a total of 187 gene sets were obtained. The target gene set comprised genes from the FA metabolism gene sets that overlapped with gene sets associated with metastasis and worse OS risk. Kaplan-Meier survival curve analysis was conducted using the Gene Expression Profiling Interaction Analysis 2 web application (http://gepia2.cancer-pku.cn/#index) to explore the association of specific genes with OS in HCC. The University of Alabama at Birmingham Cancer Data Analysis (UALCAN) portal (http://ualcan.path.uab.edu/) was utilized to analyze the correlation between promoter methylation level and protein level of the key gene mitochondrial ribosomal protein L35 (MRPL35) and clinical status. UALCAN (ualcan.path.uab.edu/analysis-prot.html) was used to analyze the association between protein level and clinical characteristics from the Clinical Proteomic Tumor Analysis Consortium (CPTAC). TNMplot (https://tnmplot.com/analysis/) (21) public database was used to explore the relationship between MRPL35 and metastatic status in HCC.

The workflow of the present study is shown in Fig. 1. Following the acquisition of the RNA-sequencing (RNA-seq) array and clinical information of patients with HCC, a weighted correlation network analysis (WGCNA) (22) was conducted using R (Vision: 3.6.3; r-project.org/). WGCNA was performed to investigate the association between gene modules and clinical variables (age, T stage, N stage and M stage). In this analysis, gene significance (GS) was determined by taking the logarithm (base 10) of the P-value resulting from the linear regression analysis conducted on gene expression and clinical features. Module significance (MS) was defined as the mean GS across all genes within a module. The module exhibiting the highest absolute MS among all the modules selected for a specific clinical feature was considered to be associated with that feature. To condense the expression patterns of all genes within a module, a singular representative expression value known as the module eigengene (ME) was calculated, which served as the primary component in the principal component analysis. By evaluating the associations between MEs and clinical features, the most pertinent gene module for each clinical feature of interest was identified. Further analysis was conducted on the module displaying the strongest association with the desired clinical features.

A combined WGCNA and FA (WF) score was calculated by multiplying the expression value of each gene by the sum of its respective coefficients. The coefficients calculated by Least Absolute Shrinkage and Selection Operator (LASSO) (23) were as follows: CYP2S1, 0.043; MRPL35, 0.178; and PLCD3, 0.061. The HCC cohort was divided into high and low WF score groups according to the median WF score. Gene Set Enrichment Analysis (GSEA; gsea-msigdb.org/gsea/index.jsp) (19) was utilized to investigate the biological functions that differed between the two groups. The GSEA employed a cut-off threshold of nominal P<0.05 and a false discovery rate (FDR) <0.25. The LinkedOmics (https://www.linkedomics.org/admin.php) (24) online tool was used to analyze the potential function of MRPL35 in HCC.

Cox regression modeling using LASSO was performed with the LASSO-Cox regression tool available in Sangerbox (25) to select the most relevant genes. LASSO regression is characterized by variable selection and complexity regularization while fitting a generalized linear model. Therefore, regardless of whether the target dependent variable is continuous or discrete, LASSO regression is used to model and predict. COX survival analysis was conducted by glmnet, survival and survminer packages of R based on clinical information.

The differential infiltration of immune cells was assessed using the CIBERSORT algorithm (https://cibersortx.stanford.edu) (26). In this analysis, the default signature matrix comprising 1,000 permutations was employed. To enhance the reliability of the findings, only data with deconvolution P<0.05 were retained. Wilcoxon rank-sum test and signed-rank test were utilized appropriately to analyze the differences between clusters and in immune cell infiltration between the high and low WF score groups.

To assess the association between MRPL35 expression and immune infiltration, the TIMER2.0 (timer.comp-genomics.org/) (27) systematic platform was used to examine immunological attributes of cancer in TCGA data.

Tumor Immune Single-Cell Hub (tisch.comp-genomics.org/home/) (TISCH) is a public source of scRNA-seq data on the tumor microenvironment (TME). TISCH collected a total of 76 sets of high-quality tumor single-cell transcriptome data from GEO and ArrayExpress and corresponding patient information, covering 27 cancer types. It provides specific annotation of cell type at a single-cell level and allows gene expression in the TME to be investigated (28).

The HepG2 liver cancer cell line (identified by Short Tandem repeat) was purchased from The Cell Bank of Type Culture Collection of The Chinese Academy of Sciences. The cells were cultured at 37°C in a humid atmosphere with 5% CO₂ and maintained in Dulbecco's modified Eagle's medium (DMEM; Gibco; Thermo Fisher Scientific, Inc.) supplemented with 10% fetal bovine serum (FBS; Gibco; Thermo Fisher Scientific, Inc.), 100 U/ml penicillin and 100 µg/ml streptomycin (Beijing Solarbio Science & Technology Co., Ltd.).

Two siRNAs targeting MRPL35 (siMRPL35-1 and siMRPL35-2) and a control siRNA (siCtrl) were synthesized and designed by Sangon Biotech Co., Ltd. to suppress the expression of MRPL35 in HepG2 cells. To initiate transfection, cells were plated in 6-well plates at a density of 1×10⁵ cells/well and incubated at 37°C until 70% confluency was reached. Subsequently, Lipofectamine^® 3000 (Invitrogen; Thermo Fisher Scientific, Inc.) was used to transfect the cells with siRNAs or siCtrl (<60 nM, following the protocols provided by the manufacturer at room temperature for 20 min. The antisense sequences of the three siRNAs were as follows: siMRPL35-1, TAGTCATATTCAGACACCAGTTGTT; siMRPL35-2, CAACTGTGTCAAGAATGCCTCTCTT; and siCtrl, TAGATACTTAGACACGACTTTCGTT.

To confirm the efficacy of siRNA transfection, RT-qPCR was utilized to detect the level of MRPL35 after transfection. Following 36 h incubation under 37°C, total RNA was isolated from cells using TRIzol^® (Invitrogen; Thermo Fisher Scientific, Inc.). Subsequently, cDNA synthesis was performed using a RevertAid First Strand cDNA Synthesis kit (Thermo Fisher Scientific, Inc.) according to the manufacturer's protocols. RT-qPCR was performed using iQ^™ SYBR Green Supermix (Bio-Rad Laboratories, Inc.). The following primers were employed for the qPCR analysis: MRPL35 forward, 5′-TTGGCATCTTCAACCTACCGC-3′ and reverse, 5′-GGAGGAAACAACTGGTGTCTGA-3′; GAPDH forward, 5′-ACAACTTTGGTATCGTGGAAGG-3′ and reverse, 5′-GCCATCACGCCACAGTTTC-3′ and the thermocycling condition was 94°C 60sec, 37°C 60 sec, 72°C for 120 sec and 20–30 cycles. To calculate the relative expression level of the target gene, the 2^−ΔΔCq method was employed (29).

In the wound healing assay, HepG2 cells were seeded into a 6-well plate and allowed to reach 80% confluence. Under the deprivation of serum the cell monolayer was then scratched using a sterile micropipette tip. The floating cells were washed away with PBS to clear away any debris. Time points of 0, 12 and 24 h were chosen to observe the progression of wound healing within the gap with light microscope. Triplicate wells were examined for each condition to ensure reliability and reproducibility. The migration and closure of the wound were assessed by measuring the width of the wound at the two time points with those at 0 h.

The Transwell assay, also known as the Boyden chamber assay, is a widely used technique for the study of cell migration and invasion. In this assay, Transwell chambers were prepared, each comprising an insert with 8-µm pores. The upper side of the membrane in the Transwell chamber was coated with a diluted Matrigel matrix in a 37°C incubator for 30 min. Prior to the assay, the HepG2 cells were cultured in 2% FBS-DMEM for 12 h to induce starvation. The cells were then suspended in FBS-free DMEM and added to the upper chamber of the Transwell chamber at a concentration of 1×10⁵ cells/chamber. Simultaneously, DMEM supplemented with 10% FBS was added to the lower compartment. The Transwell chambers were then incubated for 48 h at 37°C. Cells that had migrated to the lower surface of the filter membrane were fixed with 4% paraformaldehyde (20 min at room temperature) and stained with 0.5% crystal violet (10 min at room temperature). Cells remaining on the upper surface of the filter membrane were gently scraped off with a cotton swab. The lower surfaces of the membrane were then observed and images captured using an inverted light microscope, and the counting process was repeated three times.

Drug Response CellMiner (https://discover.nci.nih.gov/cellminer/) is as an integration platform for both molecular and pharmacological data concerning the NCI-60 cancer cell lines (30). This resource was employed to explore the correlation between MRPL35 expression and drug response.

Statistical analyses were performed using the aforementioned online databases and R software (version 3.6.3). Data are presented as mean ± SD. For the statistical analysis of experimental data, GraphPad Prism 9.0 (GraphPad Software; Dotmatics) was employed. Unpaired Student's t and Wilcoxon rank-sum test and one-way ANOVA tests (with LSD as post hoc test) were used to calculate the significance of differences in data between and among groups. P<0.05 was considered to indicate a statistically significant result.

Following the extraction of gene expression data and clinical information associated with HCC from TGCA, data on 362 cancerous samples with total prognostic information and 50 normal samples were obtained. These data were subjected to pre-processing cleanup and quality control prior to being subjected to WGCNA to investigate the associations between gene clusters and clinical characteristics associated with HCC. Based on scale-free topology, a soft-thresholding power of β=4 was chosen, indicating the reliability of the chosen power value. The rationality test yielded positive results as shown in Fig. S1A-C and following the identification of 44 gene modules (Fig. S1D), the skyblue3 and purple modules were found to exhibit significantly positive connections with M stage (cor=0.86, P<0.05) and N stage (cor=0.89, P<0.05) (Fig. S1E and F). As a result, the two modules, consisting of 957 genes, were recognized as the modules most relevant to tumor-distant metastasis and lymph node metastasis in HCC.

There were seven genes in the area of overlap of the metastasis-associated WGCNA genes with the FA metabolism-related genes from MSigDB, and genes identified to be associated with the risk of HCC progression and OS by univariate Cox regression analysis (Fig. 2A). The LASSO-Cox regression model was used to further screen the seven core genes in HCC prognosis, and three genes were identified as preferable parameters (Fig. 2B and C). The relationship between the levels of the three genes, namely phospholipase C d3 (PLCD3), MRPL35 and cytochrome P450 family 2 subfamily S member 1 (CYP2S1), and outcomes is shown in Fig. 2D. A significant difference in OS was observed between the high and low WF score groups by Kaplan-Meier analysis (P=0.009; Fig. 2E).

The univariate Cox regression analysis of the seven genes is shown in Fig. 3A, each with P<0.05 and hazard ratio (HR)>1 for OS. MRPL35 had the highest HR among the three genes of particular interest (HR=1.461, P<0.05). Survival analyses were also conducted for the three preferred genes, CYP2S1 (HR=1.5, P=0.024, Fig. 3B), PLCD3 (HR=1.4, P=0.042, Fig. 3C and MRPL35 (HR=1.9, P<0.001, Fig. 3D). Subsequent expression analysis further indicated that PLCD3 and MRPL35 were significantly upregulated in cancerous samples compared with normal tissue (P<0.05; Fig. 3E-G).

These results indicate that the analysis of samples based on PLCD3, MRPL35 and CYP2S1 expression may be utilized to predict the OS of HCC. Therefore, the mechanism underlying the difference is worthy of exploration.

The effective infiltration of lymphocytes and ability of tumor cells to escape them affects the antitumor efficacy of immunity and the outcome of patients. With the CIBERSORT algorithm, the immune cells in the microenvironment of each sample were analyzed and predicted (Fig. 4A). In the high WF score group, it was predicted that the infiltration of lymphocytes was lower (P=0.0014) and the infiltration of macrophages was higher (P=0.047) than that in the low WF score group (Fig. 4B). The specific types of immune cell were analyzed, which indicated the presence of more M0 macrophages and fewer CD8⁺ T cells in the high WF score group compared with the low WF score group (P<0.05; Fig. 4C). The impaired immune efficacy might be a reason for the worse outcome of patients with a high WF score.

To investigate the downstream pathways and characteristics, GSEA was performed for the high and low WF score groups, which demonstrated that several tumor promotive pathways were evidently activated, namely TGF-β signaling, hypoxia, KRAS signaling, IL6/JAK/STAT3 signaling, epithelial-mesenchymal transition (EMT) and G2M checkpoint pathways, each with FDR <0.25 and normalized enrichment score >0 (Fig. 5). These pathways are known to promote tumor development and immune dysfunction.

Given the clinical value and potential mechanism based on the WF score, focus was given to the MRPL35 gene, which had the highest HR value, based on the assumption that it is the main factor facilitating tumor progression. In agreement with its mRNA level, the protein level of MRPL35 in HCC tumor tissue was higher than that in normal tissue (Fig. 6A), and a significant association with metastasis status was identified using TNMplot (P<0.05, Fig. 6B). Hypermethylation of a gene promoter region is known to regulate transcriptional silencing, and promoter hypomethylation can result in the upregulation of downstream genes. The analysis performed in the present study found that the promoter methylation levels of MRPL35 in HCC were significantly lower compared with those in tumor samples (P<0.05, Fig. 6C).

LinkedOmics analysis was performed to evaluate the potential biological function of elevated MRPL35, and the results included activation of ‘citrate cycle (TCA cycle)’ and ‘fatty acid metabolism’ (Fig. S2A and B).

The associations of MRPL35 with clinical characteristics are shown in Fig. S2C-G. The results indicate that the level of MRPL35 was significantly higher in cases with residual tumor status R1 (vs. R0; Fig. S2E) and female patients (Fig. S2G), both P<0.05.

When considering whether a gene possesses value as a therapeutic target or predictive parameter, the precise location of the gene at the cell level is important. The scRNA-seq method allows the expression level of a gene to be explored among different types of cells. In the microenvironment of HCC, MRPL35 was found to be mainly expressed by proliferating T cells and malignant cells (Fig. 6D). Similarly, a negative association with CD8⁺ T cells and positive association with myeloid-derived suppressor cells (MDSCs) were observed in HCC (Fig. 6E). Furthermore, the comparison of MRPL35 IHC between tumor and normal human tissues supported the oncogenic role of MRPL3 (Fig. 6F-I).

The aforementioned findings suggest that inhibiting the expression of MRPL35 in tumor cells may impair their proliferation, migration and immune-escape abilities.

To validate the role of MRPL35 as a risk factor and its potential as a novel therapeutic target, assays were conducted to evaluate the effect of its knockdown on the invasion and migration competence of HepG2 cells. Following the confirmation of the efficacy of knockdown by the MRPL35 siRNA compared with the siCtrl (P<0.001; Fig. 7A), negative control (NC), siCtrl and MRPL35 siRNA groups were evaluated in vitro. In the Transwell assay, the two siRNAs against MRPL35 showed a significant inhibitory effect on invasion after 48 h compared with siCtrl (P<0.001; Fig. 7B and C). In addition, the wound healing assay indicated that liver cancer cells with MRPL35 knockdown possessed significantly impaired ability to migrate compared with the controls (P<0.001; Fig. 7D and E). These experimental results suggest that MRPL35 may contribute to liver cancer metastasis and the blockade of MRPL35 could be a promising therapeutic target.

To investigate the potential value of MRPL35 in drug treatment, an analysis of the correlation between the expression of MRPL35 and the response to the drug in 60 cell lines was conducted. The results indicated that MRPL35 was negatively correlated with the IC₅₀ of 10 in 16 drugs: LY-3023414, Everolimus, spebrutinib, PQR-620, JNJ-42756493, M2698, INK-128, Rapamycin, futibutinib and WORTMANNIN, all with a correlation coefficient <0 and P<0.05 (Fig. 8). This provides preliminary guidance for strategies that may be applicable to patients with high MRPL35 levels.