Endoplasmic Reticulum Membrane Protein Complex Regulates Cancer Stem Cells and is Associated with Sorafenib Resistance in Hepatocellular Carcinoma

Introduction

Hepatocellular carcinoma (HCC), the most common primary hepatic malignancy, is characterized by high recurrence, metastasis, and mortality due to its complex pathophysiology.1,2 Over half of the global incidence occurs in China, placing a significant burden on the local healthcare system.3 Sorafenib, a multi-tyrosine kinase inhibitor, remains the standard first-line therapy for advanced-stage HCC.4 However, the overall response rate to sorafenib is suboptimal due to frequent primary and acquired resistance.5–7 Advances in immunology and molecular biology have highlighted dysregulation of molecular signaling pathways and immunosuppression as key factors in sorafenib resistance, though the precise mechanisms are not yet fully understood.8 T Thus, elucidating the mechanisms underlying HCC invasion, metastasis, and sorafenib resistance is crucial for improving patient prognosis.

The endoplasmic reticulum (ER) membrane protein complex (EMC) was first identified in yeast studies as crucial for cellular transmembrane transport.9,10 EMC is composed of nine proteins arranged around a highly conserved core, notably the EMC3-EMC6 subcomplex.11 Advances in structural resolution have enhanced our understanding of membrane protein biogenesis across all life forms, revealing that EMC abnormalities are linked to various human diseases.9,12 In cancer research, EMC6 has emerged as a potential prognostic biomarker for pancreatic, lung, and gastric cancers,13,14 underscoring the growing recognition of EMC’s role in tumors.

However, research on EMC in the context of tumors remains sparse. To address this, a comprehensive analysis of EMC in cancer is crucial. In the present study, we examined the expression heterogeneity of genes encoding 10 EMC subunits across various cancers. Our analysis revealed that genomic alterations, particularly copy number variations (CNVs), impact EMC levels. We specifically investigated the role of EMC in HCC. Through Friends analysis, we identified EMC1—a gene not previously studied in HCC—as significantly associated with the disease.

Further investigations involving single-cell analysis and in vitro experiments demonstrated that EMC1 plays a critical role in maintaining cancer cell stemness and enhancing the tricarboxylic acid (TCA) cycle. Additionally, using multiple immune infiltration algorithms and spatial transcriptome analysis, we established that EMC1 is essential for M2 macrophage infiltration. Multiplex immunofluorescence (mIF) further validated the association of EMC1 with β-catenin and M2 macrophage markers. Consequently, we have delineated the EMC1/cancer stem cells (CSCs)/M2 macrophage axis.

Materials and Methods Public Datasets and Corresponding Processing

For pan-cancer analyses, we utilized the GSCALite web tool to process data from The Cancer Genome Atlas (TCGA).15 For hepatocellular carcinoma (HCC) data specifically (TCGA-LIHC), we employed the R package “TCGAbiolinks” to download RNA sequencing data (RNA-seq, FPKM values) and converted these values to Transcripts Per Million (TPM).16 Pathology image data were obtained from the Genomic Data Commons (GDC) portal.17 For datasets from the Gene Expression Omnibus (GEO),18 we assigned each HCC sample an EMC score using single-sample gene set enrichment analysis (ssGSEA) and conducted subsequent analyses using the Biomarker Exploration of Solid Tumors (BEST) database.19,20 Specifically, for sorafenib resistance data (GSE109211, GSE94550, and GSE73571), we downloaded the raw data via the “GEOquery” package and processed them into expression matrices for further analysis.21–24 The single-cell RNA sequencing (scRNA-seq) data used in this study (GSE166635) were sourced from a prior HCC study.25

scRNA-Seq Data Analyses

The filtered output of Cellranger from hepatocellular carcinoma (HCC) single-cell RNA sequencing (scRNA-seq) data, which included cells from two patients26, was analyzed in this study. A total of 22,631 cells were processed using Seurat V4.27 We utilized the uniform manifold approximation and projection (UMAP) coordinates and cell annotation information provided by the Tumor Immune Single Cell Hub (TISCH) without additional processing. Malignant cells were then distinguished from epithelial cells. The “FindAllMarkers” function with default parameters was employed to identify markers for each malignant cluster28, and dot plots were generated using the “Nebulosa” package.29

To analyze communication between malignant and myeloid cells, we utilized the ligand-receptor information from the “CellChat” package.30 Additionally, cell stemness and metabolism scores were calculated using the “CytoTRACE” and “scMetabolism” packages, respectively.31,32

Spatial Transcriptome (ST) Analysis

We downloaded spatially resolved transcriptomics data sets VISDS000514 and VISDS000511 from the CROST database, a comprehensive repository for spatial transcriptomics.33 M2 macrophage scores for each spot were assessed using the xCell method, while EMC+ epithelial cell scores were determined using single-sample gene set enrichment analysis (ssGSEA).34 The spatial overlap of EMC+ epithelial cells, M2 macrophages, and macrophage migration inhibitory factor (MIF) was visualized using the “SpatialPlot” function in Seurat V4.

Estimation of Immunological Characteristics in the Tumor Microenvironment (TME) of HCC

To analyze the correlation between EMC levels and immune profiles, we utilized the BEST database and visualized the results with heatmaps.

Calculation of the Stemness Level of Malignant Cells

In TCGA-LIHC, we calculated the mRNA-based stem index (mRNAsi) value (ranging from 0 to 1) for each HCC sample, which was strongly correlated with the characterization of the level of dedifferentiation, with higher values representing stronger stemness. Further, the correlation between the mRNAsi index and EMC1 was analyzed.

Enrichment Analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for EMC were conducted using the BEST web tool.35,36 Additionally, Metascape was employed to compare signaling differences among various clusters of malignant cells.37

Reagents

All employed antibodies and reagents are listed in Table S1.

Cell Culture

Human HCC HuH-7 and SK-Hep1 cells were obtained from Wuhan Pu-nuo-sai Life Technology Co. Ltd. Hepa1-6-luc cell obtained from Cellcook Biotech Co., Ltd. Sorafenib-resistant hepatoma cell lines (Huh7 sorafenib-resistant cells and SK-Hep1 sorafenib-resistant cells) were acquired from Shanghai Zhong Qiao Xin Zhou Biotechnology Co., Ltd. Cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM) with 10% FBS at 37°C in a 5% CO2 humid incubator. Sorafenib was added to the cell lines (Huh-7/Sorafenib and SK-Hep1/Sorafenib) at a concentration of 2 μg/mL to maintain resistance as per the manufacturer’s instructions. THP-1 cells (human monocytic cells) were purchased from the cell bank of the Chinese Academy of Sciences (Shanghai, China). THP-1 cells were cultured in Roswell Park Memorial Institute-1640 (RPMI-1640) medium supplemented with 10% fetal bovine serum (FBS) in incubator at 37°C in 5% CO2.

Lentiviral Vector Generation and Cell Transfection

Lentiviral vectors were employed for EMC1 knockdown experiments. The GeneChem Corporation synthesized the viruses. The specific information is shown in Table S2. For transient transfection, cells were transfected according to the manufacturer’s indications. In short, HCC cells underwent a 24-hour recombinant lentiviral transduction via 2 μg/mL polybrene, and the successfully transfected cells were screened with 1.5 μg/mL puromycin. Lastly, Western blot assays were used to determine EMC1 knockdown efficiencies as well as outcomes.

Western Blotting (WB)

In a previous report, the WB protocol was described.38 Proteins were extracted from lysed cells using radioimmunoprecipitation (RIPA) buffer, followed by Bradford assays. Subsequently, sodium dodecyl-sulphate polyacrylamide gel electrophoresis (SDS-PAGE) was conducted with either a 10% or 8% gel concentration to obtain 20 µg of each sample. The samples were then transferred onto a polyvinylidene fluoride (PVDF) membrane. Subsequently, the membrane was blocked using a 5% bovine serum albumin (BSA) solution. Following this, primary antibodies were incubated with the membrane overnight at 4°C. The membrane was then washed three times with an incubation solution containing Tris-buffered saline and 0.05% Tween-20. Subsequently, the membrane was labeled using secondary antibodies. The β-actin were utilized as references.

CCK8 Assay

We evaluated the IC50 of HCC using the CCK8 assay. We put 5×103 HCC cells in 96-well plates and cultured them with 5% CO2 at 37 °C for 24 hours. The cells were then treated with sorafenib in culture medium for another 24 hours. Optical density (OD) was then measured at 450 nm using a microplate reader (BioTek Synergy HT).

Animal Experiments

For the subcutaneous xenograft tumor model, 36 male BALB/c nude mice (6 weeks of age) were randomly divided into six groups. Huh-7 and Huh7 sorafenib-resistant cells (Huh7/sorafenib) were infected with/without NC or sh-EMC1 for 24 hours, then collected for subcutaneous injection (1×106 cells/injection) into athymic BALB/c nude mice. Once the subcutaneous nodules grew to the size of a rice grain (approximately one week), the subcutaneous xenograft model was successfully established. The nude mice were then treated with sorafenib by gavage at a dosage of 30 mg/kg/day for 21 days. On day 28, the mice were euthanized, and subcutaneous tumor tissues were harvested for further analysis. Tumors were excised, fixed in 10% formalin, and embedded in paraffin.

For the orthotopic liver tumor model, Hepa1-6-luc cells were infected with/without NC or sh-EMC1 for 24 hours, then collected for injection (5 × 10^5 cells/injection). Eighteen male 6-week-old C57BL/6 mice were divided into three groups: control group, NC group, and sh-EMC1 group. Mice were anesthetized using isoflurane (oxygen delivered at 0.5 L/min with 3% isoflurane for induction and 1.5% isoflurane for maintenance). After sterilizing the site, a 2 mm transverse incision was made below the xiphoid, perpendicular to the median line, and was 1–1.5 cm long. The left liver lobes were carefully pulled out of the abdominal cavity with a sterile cotton swab. A 0.2 cm incision was made in Glisson’s capsule to serve as the transplantation site. The cells were directly inoculated into the left liver lobe of the mice through a 27-gauge needle. To prevent tumor cell leakage from the injection point, a gelatin sponge patch, 2 mm in diameter, was attached to the site of injection for a few minutes after needle withdrawal.Seven days after Hepa1-6-luc cell inoculation, the mice were administered sorafenib orally at a dose of 30 mg/kg/day for 21 days. On day 28, the mice were anesthetized and intraperitoneally injected with 1.5 mg D-Luciferin in 100 μL PBS for the live imaging assay. The fluorescence images were analyzed using AniView SE Living Imaging software (BLT, Guangzhou, China). All mice were sacrificed 28 days after Hepa1-6-luc cell inoculation, and the tumors were fixed for further experiments. Additionally, to account for potential variations in tumor size, we calculated a “tumor ratio” as the percentage of tumor area to the entire liver for each mouse. This study adhered to the COPE guidelines and complied with ARRIVE guidelines 2.0 (https://arriveguidelines.org/arrive-guidelines). All animal experimental procedures were conducted in accordance with the Basel Declaration and were approved by the Animal Ethics Committee of Jiangsu Province Hospital of Chinese Medicine (Ethics approval number: 2022DW-10-01).

Hematoxylin and Eosin (HE) Staining

The histopathology of liver tissues was evaluated using HE staining. After dehydration in an ethanol gradient, liver samples were embedded in paraffin and sectioned into 4-μm pieces following immersion in a 10% formaldehyde solution. The sections were then deparaffinized, stained with hematoxylin and eosin, mounted, and examined under an upright epifluorescent microscope (Nikon, Eclipse Ni-E, Japan).

Patient’s Specimens and Clinical Data Collection

Primary liver cancer tissues and adjacent non-tumor samples were collected from 80 patients with HCC at the Affiliated Hospital of Nanjing University of Chinese Medicine. Tumors were selected based on the criterion that no prior chemotherapy or radiation therapy had been administered before surgical resection. Upon extraction, tissue specimens were rinsed with cold phosphate-buffered saline, fixed immediately in 4% paraformaldehyde, and subsequently embedded in paraffin. Histopathological examination confirmed the diagnosis of HCC in all patients following hepatectomy. Informed consent was obtained from each patient participating in this study. Additionally, clinical data were retrieved from the electronic medical records of the same hospital. This study was conducted in accordance with ethical standards and received approval from the Ethics Committee of the Affiliated Hospital of Nanjing University of Chinese Medicine (Ethics approval number: 2022NL-162-02). Table S3 shows the clinical and pathological data of patients.

Immunohistochemical (IHC) and Immunofluorescence (IF) Staining

The protocols for IHC and IF staining were adopted from previous investigations.39,40 The final concentrations of antibodies were established based on previous studies or as recommended by the antibody providers. Tissue sections were first blocked with a protein blocking solution to minimize non-specific binding. Following this, the slides were incubated with the primary antibody. Immunohistochemical (IHC) scoring, which evaluated both staining intensity and extent, was performed independently by two researchers. Images were captured using a NIKON Eclipse Ni-E microscope (NIKON, Japan). Staining intensity was classified as follows: 0 for negative, 1 for weak, 2 for moderate, and 3 for strong. The extent of staining was assessed based on the percentage of positively stained cells: 0 for no positive staining, 1 for less than 10% positive staining, 2 for 10–50% positive staining, and 3 for more than 50% positive staining. The overall IHC score, or H-SCORE, was computed by combining the staining intensity and the proportion of stained cells, using the formula: H-SCORE = ∑(PI × I) = (percentage of cells with weak intensity × 1) + (percentage of cells with moderate intensity × 2) + (percentage of cells with strong intensity × 3). In this context, PI denoted the percentage of positively stained cells relative to the total number of cells in a specified field of view, while I represented the intensity of staining. The H-SCORE ranged from 0 to 300, with a higher value indicating a more pronounced positive staining.

Tissue sections or cell climbing slides were first subjected to blocking with either 5% bovine serum albumin (BSA) in phosphate-buffered saline (PBS) or 0.1% Tween-20 to reduce non-specific binding. Following blocking, the slides were incubated overnight at 4 °C with the primary antibody. After this, the slides were treated with secondary antibodies at 25 °C for 1 hour. Nuclear staining was performed using 4’,6-diamidino-2-phenylindole (DAPI) for 3 minutes in the dark. The slides were then rinsed four times with PBS for 5 minutes each. Finally, the sections were sealed with an immunofluorescence quencher solution. Immunofluorescence staining was observed, and images were captured using an inverted fluorescence microscope (Olympus CKX-4, Japan).

Sphere Production

To evaluate sphere formation capabilities, 5 × 10³ HCC cells were seeded per well in 6-well ultralow attachment plates. The cells were incubated for 7 days in serum-free medium supplemented with 20 ng/mL human recombinant fibroblast growth factor (FGF), 20 ng/mL human recombinant epidermal growth factor (EGF), and 2% B27. The number and size of the formed spheres were subsequently analyzed using light microscopy (Olympus BX53, Japan).

Colony Formation Assays

We seeded 500 HCC cells per well in 6-well plates and incubated them for approximately 14 days. Following incubation, colonies were stained with 0.5% crystal violet at room temperature (20–25°C) for 10 minutes. Colony formation was then assessed and counted using a stereomicroscope (Nikon, MZ1500, Japan).

Transwell Assay

To assess the invasive ability of the cells, a Transwell invasion assay was conducted. In 24-well Transwell chambers, the upper compartment was filled with 200 μL of DMEM medium, and 1 × 104 cells were seeded. The lower compartment was supplemented with 500 μL of medium containing 10% FBS as a chemoattractant. A layer of Matrigel was applied to the upper compartment to simulate the extracellular matrix. After a 24-hour incubation period, the upper compartment was washed with 1% PBS. Cells were then fixed with 4% paraformaldehyde for 15 minutes and stained with 0.1% crystal violet. For analysis, cells were photographed in 5 random fields under a light microscope (Olympus BX53, Japan) and counted.

Wound Healing Assay

The migratory ability of HCC cell lines was evaluated using a wound healing assay. In each well of six-well plates, 4×105 cells were seeded and incubated for 24 hours in serum-free medium. Following incubation, the medium was removed, and wounds were created in the cell monolayer using a 200-μL pipette tip. The rate of wound closure was assessed at three time points (0, 24, and 48 hours) using an Olympus CKX-41 inverted microscope with a 200× magnification.

Transmission Electron Microscopy (TEM)

The morphology of mitochondria was examined using transmission electron microscopy (TEM). Cells were prepared according to the manufacturer’s protocol. Specimens were stained with 0.3% lead citrate and observed under an electron microscope (Hitachi, Tokyo, Japan) at magnifications of 2500× or 30,000×.

Measurement of ATP Content

HCC cells were seeded at a density of 2×105 cells per well in six-well plates. After 24 hours of incubation, the growth medium was replaced with 0.1 mL of glucose- and FBS-free DMEM, and cells were further incubated for 12 hours. ATP levels were measured using an ATP assay kit. Cells were lysed and centrifuged at 12,000 × g for 5 minutes. In a 96-well plate, 100 μL of cell lysate was mixed with 100 μL of ATP detection working solution. ATP levels were quantified using the BioTek Synergy HT Microplate Reader. Protein concentration in the lysates was determined using a BCA assay kit. The ATP level was calculated by normalizing the ATP concentration to the total protein concentration (ATP level = ATP concentration / total protein concentration).

Establishment of a Co-Culture System

To induce macrophage differentiation, THP-1 cells (1 × 105 cells/mL) were cultured in 6-well plates and treated with phorbol 12-myristate 13-acetate (PMA) for 24 hours.41 The PMA-containing medium was replaced with serum-free medium and the cells were cultured for another day. Two days prior to the co-culture experiment, cells (1 × 105 cells/mL) from the control group, knock-down group (sh-EMC1), and negative control (NC) group were seeded onto 0.4-μM transwell inserts. For co-culture, the culture medium in the inserts with HCC cells was removed to the top of the 6-well plates with differentiated THP-1 cells (macrophage). The cells were co-cultured for an additional 48 hours, and macrophage were harvested for immunofluorescence staining.

Enzyme Linked Immunosorbent Assay (ELISA)

Considering that Macrophage migration inhibitory factor (MIF) is a secreted protein, its concentration within the co-culture system was assessed utilizing the MIF ELISA Kit. The supernatant was carefully transferred from the culture into sterile vessels, followed by a 20-minute centrifugation at 2000 rpm to isolate the supernatant. The ELISA assay was conducted in adherence to the manufacturer’s instructions, with the optical density of each well measured at 450 nm using a Bio-Tek Synergy HT microplate reader immediately thereafter.

Retrospective Study

This retrospective study included 20 patients with HCC who received sorafenib monotherapy following transarterial chemoembolization (TACE) at the Affiliated Hospital of Nanjing University of Chinese Medicine from March 2022 to March 2024. All patients were pathologically diagnosed through hepatic aspiration biopsy. Lesion samples from the biopsies were evaluated by experienced pathologists in the Department of Pathology. The patient clinicopathological characteristics are summarized in Table S4. Inclusion Criteria: Patients were eligible if they: (1) had a confirmed diagnosis of HCC; (2) had a Barcelona Clinic Liver Cancer (BCLC) stage B or C; (3) had good liver function, classified as Child-Pugh scores A to B7; and (4) had an Eastern Cooperative Oncology Group (ECOG) performance status of ≤ 2. Extrahepatic metastases were permitted. Exclusion Criteria: Exclusion criteria were: (1) liver transplantation at any time; (2) presence of only nodal or distant metastases without viable liver lesions; (3) secondary malignancies; (4) prior or concurrent use of other targeted agents, chemotherapy, or immunotherapy; and (5) serious adverse reactions or drug allergies.

Treatment and Monitoring: Patients were treated with sorafenib starting at 400 mg twice daily (b.i.d). Dose adjustments were made based on toxicity, with a minimum dose of 200 mg b.i.d. on alternate days. The study protocol was approved by the ethics committee of the Affiliated Hospital of Nanjing University of Chinese Medicine (Ethics approval number: 2022NL-162-02), and written informed consent was obtained from all participants. Assessment and Response Evaluation: Radiological assessments were conducted every two months during treatment using liver computed tomography (CT) or magnetic resonance imaging (MRI), based on the Response Evaluation Criteria in Solid Tumors version 1.1. Tumor response was classified as stable if there was no significant change. Non-response was indicated by tumor progression greater than 30% over three consecutive months or evidence of extensive metastases.

Statistical Analysis

The potential relationships between variables were determined using Spearman correlation coefficient. The t-test (normally distributed data) and Wilcoxon rank-sum test (non-normally distributed data) were used for two group comparison. Additionally, Kruskal–Wallis and one-way analysis of variance (ANOVA) were employed for within-group and multiple-group comparisons. A two-tailed P-value of 0.05 was used as the significance criterion for all data analyses.

Results The Transcriptional Variations and Genetic Characteristics of 10 EMC Molecules

Information on the 10 EMC subunit-coding genes was obtained from the HUGO Gene Nomenclature Committee (HGNC) web tool. These genes include EMC1, EMC2, EMC3, EMC4, MMGT1, EMC6, EMC7, EMC8, EMC9, and EMC10 (Table S5). Previous studies have highlighted functional differences of EMC molecules across various cancers. For instance, overexpression of EMC6 in gastric cancer cells has been shown to inhibit cell growth, migration, and invasion, while knockdown of EMC3 in melanoma has led to enhanced immune responses.14,42 Thus, EMC molecules may exhibit either cancer-promoting or cancer-inhibiting functions depending on the specific environment. To examine the transcriptional variation of EMC molecules in cancer, we analyzed HCC tissue samples from the TCGA-LIHC cohort alongside normal tissues from both TCGA-LIHC and GTEx databases. As illustrated in Figure S1A, most EMC molecules were upregulated in cancer samples compared to normal controls, suggesting that EMC function may be activated in cancer.

To elucidate the cause of abnormal EMC expression, we analyzed genetic variations, including single nucleotide variants (SNVs) and copy number variations (CNVs), in pan-cancer data. We observed that EMC1 was mutated at a high frequency, particularly in endometrioid cancer (UCEC, 37%), colon cancer (COAD, 13%), skin cutaneous melanoma (SKCM, 10%), and stomach adenocarcinoma (STAD, 10%) (Figure S1B). Among 342 samples with at least one mutation in the 10 EMC genes within TCGA pan-cancer data, EMC1 exhibited the highest mutation frequency (36%), followed by EMC2 (21%) and EMC7 (12%), with missense mutations being the most prevalent. UCEC was the most frequently mutated cancer type (Figure S1C). These findings suggest that aberrant EMC expression may be linked to SNVs, particularly in EMC1.

Next, we investigated the impact of CNVs on EMC expression in pan-cancer datasets. Analysis of CNV frequency revealed substantial variability in the CNV of the 10 EMC genes across different cancers. EMC2 and EMC10 displayed the highest CNV frequencies (Figure S1D), primarily involving copy number heterozygous amplification and deletion (Figure S1E and F). We hypothesize that CNVs, in addition to SNVs, may contribute to the transcriptional dysregulation of EMC molecules. Correlation analysis between CNV and mRNA levels confirmed this hypothesis, showing a positive correlation between CNV and mRNA levels in most cancer types (Figure S1G).

In this study, we focused on the role of EMC in HCC. To accurately compare the differential expression of EMC molecules in various tissue samples, we analyzed paired mRNA expression data from HCC tumor and normal samples. Consistent with earlier findings, we confirmed that all EMC genes were expressed at higher levels in HCC tumor tissues compared to normal tissues (Figure 1A). Additionally, a high degree of transcriptional correlation among the EMC molecules was observed (Figure 1B). To further explore the expression patterns of EMC molecules in tumor versus normal samples, we performed principal component analysis (PCA). As shown in Figure 1C, principal components PC1, PC2, and PC3 effectively distinguished between normal and tumor samples in the TCGA-LIHC and GTEx datasets. Notably, paracancerous tissue samples (TCGA-LIHC) were more similar to HCC tissue samples than to normal liver tissue samples (GTEx). This suggests that EMC activation may precede malignant transformation in normal tissues.Finally, Cox regression analysis indicated that most EMC genes were associated with adverse outcomes in HCC (Figure 1D, HR > 1). These results suggest that dysregulation of EMC molecules may contribute to the development and progression of HCC.

Figure 1 Abnormal patterns of the endoplasmic reticulum (ER) membrane protein complex (EMC) subunits expression in cancer. (A) Differential expression of EMC-related molecules between paired normal and cancer samples in TCGA-liver hepatocellular carcinoma (LIHC) (paired t-test). (B) The chord diagram demonstrates the correlation among EMC-related molecules (TCGA-LIHC). Red color represents positive correlation. (C) The three-dimensional map (principal component analysis) was provided to classify the normal liver samples (GTEx) and normal and tumor samples in TCGA-LIHC based on EMC-related molecules. (D) The Cox regression analysis of the average OS rate of 10 EMC-related molecules in the TCGA-LIHC cohort. ***P < 0.001.

EMC is Activated in HCC and Associated with Poor Prognostic Features

Next, we examined the clinical and biological significance of EMC in HCC. To achieve this, we constructed an EMC score using the ssGSEA algorithm across nine publicly available datasets, which allowed us to comprehensively evaluate EMC activity (Figure 2A). We classified samples into high and low EMC groups based on the median EMC level. Consistent with the pancancer analysis, the EMC score was significantly elevated in tumor tissues across multiple datasets (E_TABM_36, GSE144269, GSE14520, and TCGA-LIHC) (Figure 2B). Additionally, patients in the high-EMC group exhibited significantly poorer outcomes compared to those in the low-EMC group. This included overall survival (OS) in TCGA-LIHC and GSE144269, progression-free survival (PFS) in TCGA-LIHC, disease-specific survival (DSS) in TCGA-LIHC, disease-free survival (DFS) in TCGA-LIHC, and recurrence-free survival (RFS) in E_TABM_36 (Figure 2C). Moreover, data from the Human Protein Atlas (HPA) database corroborated these findings, showing that the protein levels of EMC subunits were upregulated in HCC tissues (Figure 2D). Pathological examination further revealed that samples with high EMC1 expression exhibited more pronounced cytologic atypia compared to those with low EMC1 expression (Figure S2), indicating a higher degree of malignancy in the former.

Figure 2 Single sample gene set enrichment analysis (ssGSEA) of EMC. (A) The box plots show the EMC levels measured by ssGSEA across the nine publicly available hepatocellular carcinoma (HCC) datasets. (B) Differences of EMC levels between cancer samples and normal controls in five public HCC datasets. (C) Kaplan-Meier survival curves of OS (TCGA-LIHC and GSE144269), PFS (TCGA-LIHC), DSS (TCGA-LIHC), DFS (disease-free survival, TCGA-LIHC), and RFS (release-free survival, E_TABM_36) analyses according to the median expression level of EMC were generated. (D) The representative immunohistochemical (IHC) staining images of EMC subunits in normal and HCC tissues from the online Human Protein Atlas (HPA) database. (E) The Oncoplot displays the differences of mutation and CNV profiles between EMC-high and -low groups (cutoff value: median expression level of EMC). (F and G) Gene ontology (GO, (F)) and Kyoto Encyclopedia of Genes and Genomes (KEGG), (G) enrichment analysis of EMC. *P < 0.05, **P < 0.01.

We then explored the potential role of EMC in carcinogenesis. Analyzing mutational landscapes and copy number variations in the high- and low-EMC groups using the TCGA-LIHC cohort revealed that the high-EMC group exhibited a higher proportion of mutations, gains, and losses compared to the low-EMC group (Figure 2E). Notably, CTNNB1 mutations were significantly more frequent in the high-EMC group (P < 0.05). Given that CTNNB1 mutations lead to aberrant β-catenin expression—a major downstream effector of the Wnt pathway associated with poorer overall survival (OS) in various cancers43—this finding underscores the relevance of the β-catenin pathway in EMC-associated tumorigenesis. To further elucidate the functional implications of EMC, we performed biological process (BP) enrichment analysis, which indicated that EMC is primarily involved in cellular respiration, aerobic respiration, oxidative phosphorylation, and ATP synthesis coupled with electron transport. Cellular component (CC) enrichment analysis revealed that EMC’s role in tumorigenesis is associated with the mitochondrial envelope, mitochondrial inner membrane, mitochondrial protein-containing complex, and mitochondrial membrane. Molecular function (MF) enrichment analysis highlighted that EMC is enriched in protein binding, oxidoreduction-driven active transmembrane transporter activity, electron transfer activity, and structural constituent of ribosomes (Figure 2F). Additionally, KEGG pathway analysis identified that EMC is involved in several pathways, including oxidative phosphorylation, metabolic pathways, and lysosomes (Figure 2G). These findings suggest that EMC may play a significant role in tumorigenesis through its involvement in key cellular processes and pathways.

EMC Drives β-Catenin Expression and is Associated with Sorafenib Resistance

CTNNB1 mutations are recognized as aberrations associated with acquired drug resistance to several chemotherapeutic agents, including enzalutamide, platinum, and sorafenib.44–46 To explore this further, we investigated the sensitivity to chemotherapeutic drugs across different EMC levels using the Genomics of Drug Sensitivity in Cancer (GDSC) database. Our analysis revealed a positive correlation between EMC levels and the 50% inhibitory concentration (IC50) values for 30 drugs, including semagacestat, sorafenib, and nilotinib (Figure 3A). Given the significant clinical role of sorafenib in treating HCC, we focused specifically on its relationship with EMC. Further validation in the GSE109211 dataset showed that EMC levels were significantly higher in the sorafenib non-responder group compared to the responder group (Figure 3B). Additionally, we observed that EMC was highly expressed in malignant cells (Figure 3C and D and Figure S3), where CTNNB1 expression was also upregulated (Figure 3E). These findings suggest that EMC may contribute to sorafenib resistance through CTNNB1-related signaling pathways.

Figure 3 Continued.

Figure 3 Resistance-related analysis of EMC. (A) The heatmap shows the correlations between estimated IC50 values of 60 representative drugs and EMC levels. Sorafenib is highlighted by the red line. (B) The boxplot shows the differences in EMC levels between sorafenib-responder (R) and non-responder (NR) groups. (C) The Uniform Manifold Approximation and Projection (UMAP) graph shows the cell types identified in the single-cell RNA-seq (scRNA-seq) dataset GSE166635. (D) UMAP revealing the joint level of EMC genes assessed by “Nebulosa” package. Higher level is indicated with a redder color. (E) UMAP revealing the expression level of CTNNB1. Higher level is indicated with a lighter color. (F) Friends analysis of EMC-related molecules. EMC1 was identified as the key member. (G) IHC staining of our own samples for analyzing EMC1 expression in HCC versus normal liver tissues. Wilcoxon rank sum test. (H) The expression level of EMC1 among different pathologic T stages, histological grades, and pathologic stages was analyzed in TCGA-LIHC. Kruskal–Wallis Test. (I) DSS, OS, and progression-free interval (PFI) analyses according to the medium expression level of the EMC1 gene were performed using HCC cases in the TCGA-LIHC cohort. (J) Box plots of EMC1 expression in parental cells and sorafenib-resistant cells (GSE94550), non-treated spheres and sorafenib-resistant spheres (GSE73571), and sorafenib-sensitive xenograft and sorafenib-resistant xenograft (GSE73571). Wilcoxon rank sum test. (K) Western blot for confirming the differential protein levels of EMC1 between sorafenib-resistant cell lines and controls (Huh-7 and SK-Hep). T test. (L) Western blot of β-catenin in EMC1-disrupted Huh-7 and SK-Hep cells. T test. (M) Dose-response curves of HCC treated with sorafenib. The IC50 values for the control, NC, and sh-EMC1 groups. Each data point represents the mean ± SD of three independent experiments. Knockdown of EMC1 significantly increased the sensitivity of both HCC cell lines (Huh-7 and SK-Hep1) and sorafenib-resistant cells (Huh-7/sorafenib and SK-Hep1/sorafenib) to sorafenib treatment, as evidenced by the lower IC50 values in the sh-EMC1 groups compared to the NC groups. (N and O) Representative images of tumors harvested from Huh-7 (N) and Huh-7/ sorafenib (O) xenografts in BALB/c nude mice treated with sorafenib for 21 days. Tumors from control, NC, and sh-EMC1 groups are shown. Tumor volume and weight analysis at day 28 post-inoculation. Tumor volumes (left) were measured over time, and the tumor weights (right) were recorded at the end of the study. Knockdown of EMC1 significantly reduced tumor volume and weight compared to NC groups (ANOVA). (P) Representative bioluminescence images of mice with Hepa1-6-luc orthotopic liver tumors. Mice were divided into three groups: Control, NC, and sh-EMC1. The bioluminescence signal, indicative of tumor burden, was measured 28 days post-inoculation. Quantification of the total flux of bioluminescence (p/sec/cm²/sr) in the different groups. The sh-EMC1 group showed a significant reduction in bioluminescence signal compared to the NC groups. (Q) Representative images of livers harvested from the control, NC, and sh-EMC1 groups showing visible metastatic foci (indicated by white arrowheads). HE staining of liver sections from the control, NC, and sh-EMC1 groups. Dashed lines outline the metastatic regions. Scale bar = 50 μm. Quantification of the number of metastatic foci per mouse and Ratio of tumor area to total liver area (%) in the control, NC, and sh-EMC1 groups. The sh-EMC1 group showed a significant reduction in tumor area ratio compared to the NC groups. (R) Multiple immunofluorescence (mIF) staining images of EMC1, β-catenin, and CK in the HCC biopsy. Scale bar, 50 μM; nuclei (DAPI) in blue. A representative view of co-staining was shown in the enlarged images at right (scale bar, 20 μM). Co-localization (EMC1+CK versus β-catenin+CK) was quantified using the Spearman correlation coefficient (n = 80). (S) Spearman correlation between EMC1 and CTNNB1 in TCGA-LIHC. The experiment was repeated three times. ns, not significant, *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

To further investigate the molecular mechanisms of EMC, we performed an analysis using Friends to identify the most critical EMC members. As depicted in Figure 3F, EMC1 emerged as the most strongly correlated with other EMC subunits, suggesting it plays a pivotal role in EMC-related biological signaling. The protein-protein interaction (PPI) network analysis also identified EMC1 as a potential hub gene (Figure S4, red circle). Consequently, we focused on EMC1 for further investigation.

We assessed EMC1 expression in 80 pairs of HCC and adjacent non-cancerous tissues using immunohistochemistry (IHC). The results revealed a significant upregulation of EMC1 in HCC tissues (Figure 3G). Clinical analysis indicated that EMC1 expression positively correlated with T stage, tumor grade, and pathological stage. Survival analysis demonstrated that overall survival (OS), disease-specific survival (DSS), and progression-free interval (PFI) were significantly reduced in HCC patients with high EMC1 expression compared to those with low EMC1 expression (Figure 3H and I, TCGA-LIHC). These findings underscore the importance of EMC1 in HCC.

We then investigated the association between EMC1 and sorafenib resistance. As shown in Figure 3J, EMC1 was significantly upregulated in drug-resistant cell lines (GSE94550), drug-resistant spheres (GSE73571), and drug-resistant xenograft tumors (GSE73571), suggesting a strong correlation between EMC1 and secondary sorafenib resistance. For in vitro validation, we compared EMC1 protein levels between resistant and non-resistant cell lines. We found that EMC1 protein expression was elevated in sorafenib-resistant Huh-7 and SK-Hep1 cell lines (Figure 3K).

We then transfected HCC cells with sh-EMC1 and confirmed the transfection efficiency (Figure 3L). Western blot analysis showed that sh-EMC1 reduced intracellular β-catenin levels (Figure 3L), indicating that EMC1 regulates β-catenin expression. Additionally, CCK-8 assays revealed a significant decrease in the IC50 value of sorafenib in both parent and drug-resistant Huh-7 and SK-Hep1 cells treated with sh-EMC1 (Figure 3M). These results suggest that targeting EMC1 may enhance sorafenib sensitivity in HCC.

To investigate the role of EMC1 in sorafenib resistance in vivo, we established both xenograft and orthotopic transplantation models. EMC1-disrupted or control cells were subcutaneously injected into BALB/c nude mice using Huh-7 and Huh-7-SR cell lines, or orthotopically transplanted into C57BL/6 mice using the Hepa1-6-luc cell line. The mice were subsequently treated with sorafenib. Our results demonstrated that knockdown of EMC1 enhanced sensitivity to sorafenib, as evidenced by reduced tumor sizes and areas at the endpoint analysis (Figure 3N–Q).

To further validate the association between EMC1 and β-catenin, we performed multiplex immunofluorescence (mIF) on HCC tissues. The analysis revealed strong expression of β-catenin in regions with high EMC1 expression and minimal β-catenin expression in areas with low EMC1 levels (Figure 3R; EMC1+CK versus β-catenin+CK, R = 0.7465). This correlation was also supported by bulk data from TCGA-LIHC (R = 0.580, Figure 3S). In conclusion, our findings indicate that EMC1 mediates resistance to sorafenib in HCC cells.

Down-Regulation of EMC1 Impairs HCC Cell Stemness and Inhibits the Tricarboxylic Acid (TCA) Cycle

Given that EMC was upregulated in malignant cells, we conducted a clustering analysis of all malignant cells to explore the precise role of EMC. As shown in Figure 4A, the “FindClusters” function identified three distinct cell subgroups (malignant C1-C3). Differentially expressed gene (DEG) analysis highlighted the top 5 genes with high and low expression levels in these subgroups (Figure 4B and Table S6). Enrichment analysis revealed that malignant C1 and C2 were associated with a broader array of cancer-related signaling pathways compared to malignant C3, suggesting that C3 may be less malignant (Figure 4C). Notably, malignant C1 exhibited strong correlations with KRAS signaling, adipogenesis, angiogenesis, mTORC1 signaling, and xenobiotic metabolism—pathways not prominently observed in C2 and C3.

Figure 4 Downregulation of EMC impairs cancer stem cell (CSC) phenotype and inhibits the tricarboxylic acid (TCA) cycle. (A) Malignant cell populations were re-clustered into three subclusters (color coding). (B) Differential gene analysis of three malignant cell clusters (malignant C1, C2, and C3). The top five upregulated (red) and downregulated (blue) genes was shown. (C) GO and KEGG enrichment analysis of the three malignant cell clusters. (D) UMAP revealing the stemness level of malignant cells assessed by “CytoTRACE” package. Higher level is indicated with a redder color. (E) UMAP revealing the joint level of epithelial marker (EPCAM, CDH1, and KRT18). Higher level is indicated with a redder color. (F) UMAP revealing the expression level of EMC1. Higher level is indicated with a lighter color. (G) Spearman correlation between EMC1 and the stemness index (mRNAsi) in TCGA-LIHC. (H) 3D sphere-formation assays in sh-EMC1 Huh-7 and SK-Hep-1 cells. Scale bars, 50μm. (I–K) Colony formation (I), wound-healing (J), and transwell migration (K) assays in EMC1-disrupted Huh-7 and SK-Hep-1 cells. Scale bar for J, 50 μm; Scale bar for K, 20 μm. (L) Ultrastructures of mitochondria (sh-EMC1 Huh-7 and SK-Hep-1 cells) were observed by transmission electron microscope (TEM). (M) Bar plots show the changes in ATP after sh-EMC1 treatment. (N) The bubble plot demonstrates the metabolic differences between different malignant clusters. (O) UMAP revealing the levels of TCA cycle in different malignant clusters. Higher level is indicated with a redder color. (P) Spearman correlation between EMC1 and the key enzymes of TCA cycle (TCGA-LIHC). (Q) Western blot of the two most relevant enzymes in EMC1-disrupted Huh-7 and SK-Hep cells. T test. The experiment was repeated three times. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Abbreviation: ns, not significant.

CytoTRACE analysis, which assesses dedifferentiation and stemness properties, showed that malignant C1 and C2 had significantly higher CytoTRACE scores than C3 (Figure 4D). This suggests that C1 and C2 may represent potential cancer stem cell (CSC) populations, further supported by their low levels of epithelial markers (Figure 4E, EPCAM+CDH1+KRT18). Additionally, EMC1 was upregulated in malignant C1 and C2, indicating its potential role in maintaining HCC cell stemness (Figure 4F). Bulk data also showed a weak positive correlation between EMC1 transcript levels and tumor stemness (Figure 4G, R = 0.12).

To further investigate the impact of EMC1 on stemness characteristics, we performed in vitro experiments. Sphere formation assays demonstrated that downregulation of EMC1 impaired cancer cell stemness and reduced sphere formation in Huh-7 and SK-Hep1 cells (Figure 4H). Moreover, clone formation (Figure 4I), wound healing (Figure 4J), and transwell migration (Figure 4K) assays revealed that EMC1 interference significantly inhibited cell proliferation, migration, and invasion.

There is evidence that β-catenin deficiency leads to mitochondrial dysfunction in hepatocytes, including impaired tricarboxylic acid (TCA) cycle activity and oxidative phosphorylation (OXPHOS), ultimately resulting in reduced adenosine triphosphate (ATP) production.47 Tumor cells often exhibit impaired TCA cycle function and rely primarily on glycolysis for ATP production. However, they still depend on mitochondria-associated reductive carboxylation to metabolize glutamine.48,49 Mitochondria-endoplasmic reticulum interactions are crucial for regulating mitochondrial respiration and bioenergetics in cancer cells.50 This prompted us to investigate whether downregulation of EMC1 could influence mitochondrial structure and function.

Transmission electron microscopy revealed that sh-EMC1 treatment impaired mitochondrial structure (Figure 4L), and suppression of EMC1 resulted in reduced ATP production in HCC cells (Figure 4M). We then assessed the metabolic profiles of each malignant subpopulation using the “scMetabolism” package and visualized the scores for the top 30 metabolic pathways. The levels of glycolysis, pyruvate metabolism, oxidative phosphorylation, lipid metabolism (eg, steroid biosynthesis), and the TCA cycle were notably higher in malignant C1 compared to other malignant clusters (Figure 4N), consistent with the characteristics of CSCs.

Focusing on the TCA cycle (Figure 4O), we found a significant positive correlation between EMC1 levels and the transcript levels of two TCA enzymes: Aconitase 2 (ACO2, R = 0.493) and Oxoglutarate Dehydrogenase (OGDH, R = 0.473) (Figure 4P, TCGA-LIHC). Western blotting further confirmed that downregulation of EMC1 resulted in decreased expression of ACO2 and OGDH in Huh-7 and SK-Hep1 cells (Figure 4Q). These findings suggest that EMC1 regulates mitochondrial metabolic pathways in HCC cells.

EMC1 Promotes M2 Macrophage Phenotype

CSCs interact with the immune system, thereby promoting immune evasion.51 In this study, we further explored the potential oncogenic mechanisms of EMC1 through alterations in the immune microenvironment. Initially, we employed eight immune estimation algorithms to assess the association between EMC and immune cell infiltration. As illustrated in Figures 5A and S5, there is a notable correlation between EMC levels and the abundance of M2 macrophages (highlighted in the red box). Subsequently, we conducted a co-expression analysis using HCC datasets to determine the correlations between EMC levels and immunomodulator-related genes. The heatmap (Figure 5B and Figure S5) clearly shows that β2-Microglobulin (B2M) and C-C Motif Chemokine Ligand 15 (CCL15) are co-expressed with EMC (highlighted in the red box). Both B2M and CCL15 have been widely demonstrated to induce M2-like polarization of macrophages. To further investigate, we performed cell communication analyses to study the interactions between EMC-represented CSC populations and myeloid cells. Our findings suggest that malignant C1 cells may interact more intensively with myeloid cells compared to other malignant groups (Figure 5C). Further ligand-receptor analyses indicated that malignant C1 cells may frequently communicate with myeloid cells through MIF signaling, leading to phenotypic changes in these cells (Figure 5D, highlighted in the red box). Numerous studies have emphasized the role of MIF signaling in the TME, particularly in suppressing innate immunity.52 Our analysis also identified that SPP1 signaling predominantly originates from malignant C1 cells (Figure 5E and F). To validate the cellchat analysis results, we obtained spatial transcriptional data from the CROST database to characterize the spatial overlap of EMC+ epithelial cells, M2 macrophages, and MIF signaling within HCC cancer tissues. Consistently, a portion of EMC+ epithelial cells, M2 macrophages, and MIF were found to be concentrated at the same locations (Figure 5G–I), indicating their physical proximity. To examine whether EMC1 contributes to the maintenance of the M2-like phenotype of macrophages in vitro, we established a co-culture model as depicted in Figure 5J. After 48 hours of co-culture with sh-EMC1-treated HCC cells, macrophages exhibited a reduced M2-like profile (Figure 5K, indicated by decreased CD163 expression). Additionally, ELISA analysis of the co-cultures showed that MIF secretion levels in the sh-EMC1 group were significantly lower than those in the control/NC group (Figure 5L). For in vivo validation, we performed IHC on the aforementioned transplant tumors. As shown in Figure 5M, sh-EMC1 treatment significantly reduced CD163 expression levels. Multiplex immunofluorescence (mIF) staining of HCC tissue samples further demonstrated that the infiltration level of M2 macrophages near tumor tissues with high EMC1 expression was also higher (Figure 5N; EMC1+CK versus CD163, R = 0.9441). In conclusion, EMC1 appears to help maintain the M2-like phenotype of macrophages through MIF signaling.

Figure 5 EMC1 promotes M2 macrophage phenotype through MIF. (A) Scatter plots show the correlation between EMC (ssGSEA) and M2 macrophages estimated by multiple algorithms. (B) Scatter plots show the correlation between EMC (ssGSEA) and C-C Motif Chemokine Ligand 15 (CCL15) as well as Beta-2-Microglobulin (B2M). (C) The circle plots show the number (left) or strength (right) of interactions in the cell-cell communication network. (D) The bubble plot shows that malignant C1 may communicate with myeloid cell frequently through MIF signals (red box). (E and F) The heatmap (E) and chord diagram (F) demonstrate that malignant C1 is the major sender of MIF signaling. (G) Hematoxylin and eosin (H&E) staining of tissue sections of ST spots in HCC samples VISDS000514 (left) and VISDS000511 (right). (H) Unbiased clustering of ST spots in HCC samples VISDS000514 (above) and VISDS000511 (below). (I) Spatial overlap of signature score of M2 macrophages, EMC+ Epithelial cells, and MIF in tissue sections. (J) Co-culture model diagram of cancer cells (above) and macrophages (below). (K) IF staining of M2 marker CD163 in macrophages with the intervention of EMC1-disrupted Huh-7 and SK-Hep cells. T test. (L) Enzyme-linked immunosorbent assay (ELISA) of MIF in the co-culture system. (M) Representative IHC images showing CD163 expression in Huh-7, Huh-7/sorafenib, and Hepa1-6-luc mouse model tumor tissues. The control, NC, and sh-EMC1 groups are displayed for each cell line. Brown staining indicates CD163 expression, while blue staining represents the nuclei. Scale bar = 10 μm. Quantitative analysis of CD163 expression levels in mouse model tumor tissues, demonstrating a significant reduction in CD163expression in the sh-EMC1 group compared to the NC groups. Data are presented as H-scores (mean ± SD) for each group. (N) mIF staining images of EMC1, CK, and CD163 in the HCC sections. Scale bar, 50 μM; nuclei (DAPI) in blue. A representative view of co-staining was shown in the enlarged images at right (scale bar, 20 μM). Co-localization (EMC1+CK versus CD163) was quantified using the Spearman correlation coefficient (n = 80). *P < 0.05, **P < 0.01, ***P < 0.001.

Abbreviation: ns, not significant.

The Expression of EMC1 and CTNNB1 is Positively Correlated with CD163 Activity in Human HCC Patients

We further examined the correlation between EMC1, CTNNB1, CD163 activity, and sorafenib resistance. In the TCGA-LIHC dataset, we observed positive correlations between EMC1, CTNNB1, and CD163 (Figure 6A). Interestingly, this positive correlation persisted in the non-responder group of the GSE109211 dataset, while it disappeared in the responder group (Figure 6B and C). To verify this phenomenon, we performed multiplex immunofluorescence (mIF) staining of EMC1, CTNNB1, and CD163 in 80 cases of human primary HCC. As shown in Figure 6D and E, the subjects were divided into two groups based on their response to sorafenib: non-responders and responders (see Methodology section for details). The results indicated a strong positive correlation between the expression of EMC1, CTNNB1, and CD163 in the non-responder group, whereas no significant correlation was observed in the responder group (Figure 6F). Additionally, the degree of co-localization of these three markers was higher in the non-responder group compared to the responder group (Figure 6G). We also examined the correlation between the combined expression levels of these markers and the response to sorafenib treatment. The integrated expression levels (EMC1 and CTNNB1 expression in CK-positive regions and overall CD163 expression) were lower in responder patients than in non-responders (Figure 6H). Patients with high integrated expression levels had lower response rates to sorafenib than those with low integrated expression levels (Figure 6I). Figure 6J demonstrates the imaging of representative lesions.In conclusion, the positive correlation between EMC1, CTNNB1, and CD163, particularly in non-responders, suggests that their combined expression may be associated with sorafenib resistance in HCC.

Figure 6 The expressions of EMC1, β-catenin, and CD163 are positively correlated in Sorafenib-resistant HCC samples. (A–C) Circle plots show the Spearman correlation among EMC1, CTNNB1, and CD163 in TCGA-LIHC (A), GSE109211 non-responder samples (B), and GSE109211 responder samples (C). Red represents positive correlation and blue represents negative correlation. (D and E) mIF s

留言 (0)

沒有登入
gif