BNIP3 Upregulation Characterizes Cancer Cell Subpopulation With Increased Fitness and Proliferation.

BNIP3 is a BH3-only protein with both pro-apoptotic and pro-survival roles depending on the cellular context. It remains unclear how BNIP3 RNA level dictates cell fate decisions of cancer cells. Here, we undertook a quantitative analysis of BNIP3 expression and functions in single-cell datasets of various epithelial malignancies. Our results demonstrated that BNIP3 upregulation characterizes cancer cell subpopulations with increased fitness and proliferation. We further validated the upregulation of BNIP3 in liver cancer 3D organoid cultures compared with 2D culture. Taken together, the combination of in silico perturbations using public single-cell datasets and experimental cancer modeling using organoids ushered in a new approach to address cancer heterogeneity.

Introduction

The heterogeneity of cancer is a well-known phenomenon that poses a daunting challenge for effective treatment. Cell-to-cell variability in signaling pathways and transcription factor activities render the whole cancer cell population only partially responsive to most drugs (1, 2). The design of a better combination targeting strategy relies on the accurate identification of key genes and pathways that define cancer cell subpopulations with increased cancer hallmarks.

The ability of cancer cells to elicit uncontrolled proliferation and evade apoptosis requires a healthy mitochondrial network maintained through coordinated fission and mitophagy (3). BNIP3 is involved in cellular responses to a multitude of different stresses through either apoptotic or non-apoptotic cell death (4). BNIP3 also serves as an autophagy receptor that plays crucial roles in the removal of damaged mitochondria via interaction with ATG8. We have previously shown that phosphorylation of S17 and S24 in the LC3 interacting domains dictates whether BNIP3 signals apoptosis or pro-survival mitophagy (5). However, it is still unclear how the RNA expression level of BNIP3 dictates cell fate decisions of cancer cells at the single-cell level.

Single-cell RNA sequencing (scRNA-seq) has been harnessed to gain important insights into cancer heterogeneity and resulted in overwhelmingly rich datasets (6). Almost all solid tumors and hematological malignancies have been investigated with scRNA-seq. Those datasets enabled the possibility to perform in silico perturbation experiments with single-cell resolution to investigate the functional significance of genes of interest (7).

Here, we undertook a comprehensive analysis of BNIP3 expression and functions in single-cell datasets and The Cancer Genome Atlas (TCGA) datasets. We identified a cancer cell subpopulation characterized by upregulated BNIP3 in most epithelial malignancies. We also interrogated the pathway alterations in cancer cells with upregulated BNIP3 expression with a quantitative pathway enrichment approach using gene set variation analysis (GSVA) (8). Our study underscored the power to combine computational and experimental approach to address gene-centered cancer heterogeneity.

Results

BNIP3 expression was first investigated in the tumor and normal samples from the TCGA and the Genotype-Tissue Expression (GTEx) projects. Using transcripts per million reads normalization, BNIP3 expression was investigated in cancer samples and paired normal samples across different cancer types ( ). The highest BNIP3 expression was found in Kidney Renal Clear Cell Carcinoma (KIRC), while significant patient-to-patient variability in BNIP3 was noted. Those population averaged measurements were incapable of capturing the intratumoral heterogeneity reflected by cell-to-cell variability of cancer cells and the complex tumor ecosystem. Single-cell transcriptomic datasets were used to determine the heterogeneous BNIP3 expression in cancer cells. Due to the inherent technical constrains of scRNA-seq, dropouts (zero UMI detected) were common. Considering the technical dropouts, cancer cells were stratified based on whether at least one UMI is detected whenever UMI count datasets were available. BNIP3 positivity actually might reflect BNIP3 upregulation. In almost all patients, scRNA-seq data revealed a cancer cell subpopulation with BNIP3 positivity.

The survival analysis was performed with all cancer types in the TCGA project ( ), suggesting BNIP3 mRNA expression as a worse prognostic factor also for cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC), cholangiocarcinoma (CHOL), and sarcoma (SARC). However, BNIP3 upregulation appeared to be a better prognosis indicator in kidney renal clear cell carcinoma (KIRC) and low-grade glioma (LGG).

Figure 1

(A) Prognostic significance of BNIP3 in the TCGA cohort. Highlighted squares indicate p-value smaller than 0.1. (B) Volcano plot showing the differentially expressed genes between BNIP3-positive and -negative cancer cells. (C) Top pathways enriched for BNIP3 upregulated genes shown as barplot. (D) The top transcription factors enriched for BNIP3 upregulated genes. (E) Top protein–protein interaction modules enriched for BNIP3 upregulated genes.

The functional significance of BNIP3 in cancer cells was first investigated using a single-cell dataset derived from head and neck cancer (9). Cancer cells were stratified by BNIP3 RNA expression. The differentially expressed genes between BNIP3-positive and BNIP3-negative cancer cells were shown ( ). The top pathways enriched for BNIP3 upregulated genes included formation of the cornified envelope, Nuclear factor erythroid 2-related factor 2 (NRF2) pathway, and response to wounding ( ). NRF2 is a transcription factor associated with antioxidant responses in cells. Interestingly, the top transcription factor regulating BNIP3 upregulated was HIF1A ( ), in agreement with the involvement of BNIP3 in cellular response to hypoxia. Using BNIP3 upregulated genes, protein–protein interaction network was constructed and analyzed for core modules. NRF2 pathway and metabolic reprogramming were among the enriched core modules, suggesting that cancer cells with higher expression of BNIP3 might have achieved increased fitness by multiple pathways ( ).

To gain a quantitative insight into BNIP3-associated pathways, we employed GSVA to investigate the differential pathway activity of BNIP3-positive and -negative cancer cells ( ). BNIP3-positive cancer cells have an upregulated activity in reactive oxygen species pathway, oxidative phosphorylation, as well as MYC targets. The high ROS burden within BNIP3-positive cancer cells might explain the feedback activation of antioxidant transcription factor NRF2. Cell cycle phase at the single-cell level was inferred using single-cell RNA, suggesting that the percentage of cells in S phase is higher in BNIP3-positive cells ( ). BNIP3-positive cancer cell subpopulation was also detected in lung cancer ( ). Next, a lung cancer dataset with 42 patients was integrated with CCA or harmony algorithm and employed to obtain BNIP3 altered gene lists and pathways. The differentially expressed genes with p-value < 0.05 and detection rate higher than 25% were identified with Wilcoxon test. Interestingly, BNIP3 upregulated genes were enriched for response to hypoxia, response to oxygen levels, response to decreased oxygen levels, and response to oxidative stress after CCA integration ( ). Response to oxidative stress was also among the top pathways enriched after harmony integration ( ).

Figure 2

(A) Heatmap of the hallmark pathways at the single-cell level. Each row represents one pathway and each column represents one cell. (B) Distribution of cell cycle phases for BNIP3-positive and BNIP3-negative cancer cells. (C) Pathway enrichment for BNIP3 upregulated and downregulated genes in lung cancer, using CCA integration. (D) Pathway enrichment for BNIP3 upregulated and downregulated genes in lung cancer, using harmony integration.

Next, we investigated BNIP3 expression in a cervical cancer single-cell atlas. BNIP3-positive cervical cancer cells displayed a shifted transcriptional signature ( ). Cervical cancer patients with high BNIP3 expression in the TCGA cohort had a significantly decreased overall survival as compared with those with low BNIP3 expression ( ). The top 3 pathways enriched for BNIP3 upregulated genes were HIF1 TF pathway, response to wounding, and Vitamin D receptor pathway ( ). The top 3 transcription factors regulating the upregulated genes were HIF1A, SP1, and RELA ( ). The proportion of BNIP3-positive cells is higher in breast cancer cells as compared to normal breast epithelial cells ( ). HER2-positive and triple-negative breast cancers seemed to have an increased proportion of BNIP3-positive cancer cells as compared with ER-positive breast cancers ( ). Of note, BNIP3 is mostly expressed by epithelial cells, but not immune cells in the tumor microenvironment ( ). The prognostic significance of BNIP3 in breast cancer patients was also investigated in the TCGA breast cancer cohort. Patients with high BNIP3 expression had a significantly worse prognosis compared with patients with low BNIP3 expression ( ).

Figure 3

(A) The percentage of BNIP3 high and low cancer cells in cervical cancer and the volcano plot visualizing the differentially expressed genes between BNIP3 high and low cervical cancer cells. (B) Survival curve for cervical cancer patients in the TCGA cohort, stratified by BNIP3 mRNA expression. (C) Top pathways enriched for BNIP3 upregulated genes shown as barplot. (D) The top transcription factors enriched for BNIP3 upregulated genes. (E) Proportion of BNIP3-positive and -negative cells in cancer epithelia and normal epithelia of the breast shown visualized as stacked barplot. Each bar indicates one individual. (F) Proportion of BNIP3-positive cells and negative cells in major cell types in the breast cancer cell atlas. (G) Survival curve for breast cancer patients in the TCGA cohort, stratified by BNIP3 mRNA expression.

The expression of BNIP3 was also investigated in a normal liver cell atlas. BNIP3 was mostly expressed by hepatocytes in the liver, but not by immune cells or stromal cells ( ). BNIP3-positive hepatocytes appeared to have a more active cycling feature, as evidenced by an increased proportion of hepatocytes in S and G2M phase ( ). In the TCGA liver cancer cohort, we did observe that liver cancer patients with high expression of HIF1A or NRF2 (NFE2L2) tend to have a worse prognosis (p-value < 0.1) ( ). The expression of HIF1A and NRF2 was highly correlated in liver cancer samples from the TCGA cohort ( ).

Figure 4

(A) Liver cell atlas visualized in UMAP plot, the intensity of color indicating expression of BNIP3. (B) Dotplot visualization of BNIP3 in major cell types within the liver. (C) Distribution of cell cycle phases for BNIP3-positive and BNIP3-negative hepatocytes. (D) Survival analysis of hepatocellular carcinoma patients in the TCGA cohort, stratified by mRNA expression of HIF1A. (E) Survival analysis of hepatocellular carcinoma patients in the TCGA cohort, stratified by mRNA expression of NFE2L2. (F) Correlation between the expression of HIF1A and NFE2L2 in liver cancers in the TCGA cohort.

Cancer cells cultured as organoids could better represent cancer cells grown in vivo and were shown to harbor increased stemness compared with cancer cells in 2D culture. We hypothesized that cancer cells might upregulate BNIP3 as a means to increase fitness when monolayer cell lines were converted into organoid lines. To validate this hypothesis, the HepG2 cell line was used as a parental cell line to establish a liver cancer organoid line ( ). HepG2 2D culture and organoid culture were subjected to bulk RNA-seq. The similarity matrix derived from RNA-seq data indicated a global change of transcriptome from 2D culture to 3D organoid culture ( ). Both oxidative phosphorylation and reactive oxygen species pathway increased in HepG2 organoids compared with HepG2 2D culture ( ). The expression of BNIP3 was significantly upregulated in the 3D culture of HepG2 as compared with 2D culture (p < 0.05). This upregulation was not observed for BCL2 Interacting Protein 3 Like (BNIP3L). Interestingly, liver cancer cells cultured as organoids have a significantly upregulated CD24 expression (p < 0.05), which played important roles in evasion from phagocytosis of cancer cells from macrophages ( ).

Figure 5

(A) Images of HepG2 cancer cells cultured in 2D culture or organoid culture. (B) Heatmap of the correlation matrix between individual cancer transcriptomes derived from 2D culture or organoid culture. (C) GSVA of hallmark pathways for individual cancer samples. Each row represents one hallmark pathway and each column represents one sample. Both rows and columns were arranged by hierarchical clustering. (D) Boxplots showing the expression of indicated genes for HepG2 cultured in 2D or organoids.

Discussion

As carcinogenesis progresses, cancer cells are making important decisions of life and death constantly. Cancer cells have unlocked the secret of phenotypic plasticity represented by distinct subpopulations with genetic or epigenetic variability. Identification of key genes and pathways that serve as master regulators of cancer cell fate decisions is key for the design of optimal treatment strategy. Our study unraveled BNIP3 upregulation as a hallmark characterizing cancer cell subpopulation with increased fitness and proliferation.

Single-cell RNA sequencing has been applied by the research community to gain insights into cancer heterogeneity and cellular ecosystem. The enormous datasets generated so far would serve as a gold mine to identify key regulators of cancer cell fate decisions if carefully reanalyzed and integrated.

Interestingly, the cancer type with the highest BNIP3 expression is clear cell renal cell carcinoma (ccRCC). This is in agreement with the fact that HIF is no longer degradable due to the loss of tumor suppressor VHL in ccRCC (10). It has been demonstrated in vitro that siRNA-mediated downregulation of BNIP3 very effectively reduced the colony-forming capacity of RCC cells (11). BNIP3 overexpression has also been shown to enhance tumor growth for lung cancer (12) and liver cancer (13). In liver cancer cells, BNIP3 was proposed to be a therapeutic target for cancer metastasis as BNIP3 upregulation enhanced anoikis resistance of HCC cells.

Controversial results have been reported regarding the role of BNIP3 in breast cancer. BNIP3 deletion in triple-negative breast cancer promoted the metastasis of disease by deregulating mitophagy (14). On the contrary, it has been demonstrated that BNIP3 promoted the malignant phenotypes of breast cancer cells under hypoxia (15). Other studies made a distinction between ductal carcinoma in situ (DCIS) and invasive carcinoma, suggesting that BNIP3 upregulation was correlated with higher risk of relapse and shorter disease-free survival only in DCIS (16).

Cancer cells have harnessed the built-in cellular programs to adapt to hypoxia, which is a common feature of tumor microenvironment. The hypoxic niches typically render chemotherapy (17) or radiation therapy (18) ineffective. Targeting HIF-2a with belzutifan (MK-6482) has been quite successful in a recent phase II trial, achieving a 49% objective response rate in patients with renal cell carcinoma (19). Another key transcription factor, NRF2, underlying BNIP3 upregulated cancer cell subpopulation has also recently been indirectly targeted with a chemical proteomic approach (20).

Our study suggested that BNIP3 might be involved in the enhanced tumorigenicity of liver cancer cells. This is consistent with a previous report that BNIP3 protects HepG2 cells against etoposide-induced cell death under hypoxia (21). Furthermore, BNIP3 upregulated cancer cells might be armed with immune evasion arsenals. Our results have demonstrated that CD24, a “don’t eat me” signal, has been upregulated in liver cancer organoids together with BNIP3. It has also been shown that a hypoxia-inducible factor elevated the expression of PD-L1 in ccRCC cells (22).

Taken together, the systems biology approach marrying in silico perturbations using public single-cell datasets and experimental cancer modeling using organoids in our study unraveled a cancer cell subpopulation characterized by BNIP3 upregulation and revealed the potential druggable master regulators of enhanced fitness and proliferation.

Methods

Processing of Single-Cell Datasets

For single-cell datasets, annotations (meta data) from the original publications were used whenever possible. For GSE131907, “Malignant cells” as defined by original researchers were considered as cancer cells and used in our analysis. For GSE168652, cells with the number of detected genes (nFeature_RNA) between 500 and 7,500 were retained. The upper limit of total UMI count was set as 50,000 to remove potential doublets and multiplets. Cells with more than 20% of mitochondrial RNA detected were also removed from our analysis. For datasets without meta data, quality control and unsupervised clustering were performed with Seurat. The count data were normalized using the “LogNormalize” method with a scaling factor of 10,000. The top 2,000 most variable genes were selected using the “vst” method. For cancer cell grouping based on BNIP3 expression, cancer cells with at least one UMI detected for BNIP3 were considered as BNIP3-positive.

Identification of Differentially Expressed Genes

Upregulated genes in each cell cluster were identified using the “FindMarkers” function with the statistical test method “wilcox”. Only genes expressed in more than 25% of cells and altered with log2FC higher than 0.25 were retained for further analysis.

Inference of Cell Cycle Phase From Single-Cell Data

Cell cycle scoring with single-cell transcriptomic data was performed with the “CellCycleScoring” function in Seurat. Each cell is assigned a score based on expression of G2/M markers and S phase markers. Cell cycle phase was predicted based on the respective cell cycle scores (G1, S, and G2M). The genes used for cell cycle scoring can be found in cc.genes.updated.2019 originally derived from a melanoma study (23).

Gene List Analysis

Differentially expressed genes with |log2FC| higher than 1 and p-value smaller than 0.05 were subjected to gene list analysis using metascape (24), including pathway enrichment, analysis of protein–protein interaction, and inference of transcription factors. Default parameters were used for implementation.

Gene Set Variation Analysis

GSVA was implemented with the GSVA package in R. The hallmark pathways and KEGG pathways were retrieved from MSigDB. For transcript per million reads (TPM) expression data, “Gaussian” was used as the kernel for the non-parametric estimation of the cumulative distribution function of expression levels. For single-cell datasets, the normalized data slot from the RNA assay was used as input for GSVA implemented also using “Gaussian” as the kernel for the non-parametric estimation of the cumulative distribution function of expression levels.

SCENIC Analysis

SCENIC (25) was implemented with pySCENIC software. Transcription factors and corresponding target genes (regulon) were inferred based on co-expression of genes across cells. In brief, SCENIC infers TFs and their target genes from correlations between the expression of genes across cells. A TF and its target genes are defined as a regulon. The regulons are then refined by pruning targets based on enriched motifs. Finally, the activity of a regulon is measured by an AUCell value in each single cell. A high AUCell value indicates high activity and enrichment of a regulon in a cell.

Transcription Factor Scoring

The bulk RNA-seq data from HepG2 2D culture and organoid culture were analyzed by a method previously developed for global transcription factor activity scoring (26). For each transcription factor, the target genes with known regulation modes were extracted from the TTRUST database (27), resulting in a list of genes activated by the transcription factor and a list of genes repressed by the transcription factor. The ratio between the median expression level of an activated target gene and the median expression level of a repressed target gene was calculated for each transcription factor and log2 transformed to obtain a final transcription factor score.

TCGA/GTEx Data Mining

Investigation of BNIP3 expression in cancer samples and normal samples from TCGA or GTEx consortium was performed with GEPIA2 (http://gepia2.cancer-pku.cn/) (28). For survival map analysis, the significance level of 0.05 was used and the median expression was used to stratify patients into a high-expression group and a low-expression group. In total, 33 different cancer types from the TCGA project were investigated.

Cell Culture

HepG2 cells were seeded in a 10-cm culture dish and maintained in DMEM medium (L110KJ, BasalMedia) supplemented with 10% FBS. Medium was renewed every 2 days. For derivation of organoid line, HepG2 cells were centrifuged at 500 g for 5 min at 4°C. The cell pellet was resuspended in Matrigel (R&D, 3533-005-02). For one well of a 24-well plate, 50 μl of cell suspension with 10,000 cells was seeded for the Matrigel to solidify. After Matrigel solidification, 1 ml of medium was added to each well. The organoid medium A contained 1% PS, 1% Glutamax, 10 mM HEPES, B27 (1:50), N2 (1:100), 1.25 mM n-Acetyl-L-cysteine, 10 mM nicotinamide, 10 nM recombinant human Gastrin I, 50 ng/ml recombinant human EGF, 100 ng/ml recombinant human FGF10, 25 ng/ml recombinant human HGF, 10 μM Forskolin, 5 μM A8301, 10 μM Y27632, and 3 nM dexamethasone. The organoid medium B contained 1% PS, 1% Glutamax, 10 mM HEPES, B27 (1:50), N2 (1:100), 1.25 mM n-Acetyl-L-cysteine, 10% Rspo-1 supernatant, 10 mM nicotinamide, 10 nM recombinant human Gastrin I, 50 ng/ml recombinant human EGF, 100 ng/ml recombinant human FGF10, 25 ng/ml recombinant human HGF, 10 μM Forskolin, and 5 μM A8301.

RNA Sequencing

HepG2 cultures were subjected to RNA extraction. After quality control with gel electrophoresis and Agilent 2100, mRNA were captured with beads coupled with oligo(dT) and fragmented before priming with random hexamers. First-strand and second-strand cDNA were synthesized and purified. The purified double-stranded cDNA were subjected to end repairing, A-tailing, and adapter ligation. The products were purified and size-selected before final PCR amplification. The PCR products were purified to obtain the final libraries, which were sequenced with Nova-seq 6000 to obtain 6G data for each sample. The raw reads were pre-processed and filtered before alignment to hg38 reference genome. Stringtie was employed to derive TPM expression matrix (29).

Statistical Analysis

All statistical analyses were performed with R. No statistical analysis was employed to estimate the sample size for desired statistical power. The identification of markers distinguishing different clusters of cells was performed with “wilcox” test, with 0.05 set as the cutoff for statistical significance. Multiple tests were corrected with the “BH” method. The statistical difference between survival curves for different patient groups stratified by BNIP3 expression level was tested with log-rank test, with 0.05 set as the cutoff for statistical significance. The difference in gene expression between 2D and 3D HepG2 cultures was tested with t-test, using 0.05 as the significance level.

Data Availability Statement

The data presented in the study are deposited in the Zenodo repository, accession number 6481921. The single cell datasets analyzed in this study can be accessed from GEO repository with the following accession numbers: head and neck cancer (GSE103322) (9), lung cancer (GSE131907) (30), lung cancer (GSE148071) (31), breast cancer (GSE176078) (32), cervical cancer (GSE168652) (33), normal liver (GSE124395) (34).

Author Contributions

YZ, YW, and XH designed and supervised the study. YZ, BC, JY, WZ, NS, YW, and XH performed data analysis and result interpretation. PD performed the experiments. YZ, YW, and XH wrote the manuscript. All authors participated in discussion and manuscript editing. All authors contributed to the article and approved the submitted version.

Funding

This study is supported by the National Science Foundation (Grant ID: 81903106).

Conflict of Interest

XH, BC, PD, and NS were employed by Intelliphecy.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.