OPA1 supports mitochondrial dynamics and immune evasion to CD8⁺ T cell in lung adenocarcinoma

Acquisition and processing of TCGA database

Transcriptome RNA sequencing (RNA-seq) data of 1,111 NSCLC samples (LUAD: 513 tumor tissues and 59 matched adjacent tissues; lung squamous cell carcinoma, LUSC: 489 tumor tissues and 50 matched adjacent tissues) were obtained from the TCGA database (https://portal.gdc.cancer.gov/) with level 3. The association between OPA1 and MFN1 expression and survival in all tumor tissue samples was analyzed by Cox regression, which was obtained from http://www.oncolnc.org/cancer/.

Immunohistochemical analysis (IHC)

All of the tissues were handled with the following steps: (1) Deparaffinizing and rehydrating the paraffin section; (2) Antigen retrieval; (3) For cooling to room temperature before proceeding, the sections were placed in PBS (pH = 7.4) and shaken on a decolorization shaker 3 times for 5 min each; (4) Blocking endogenous peroxidase activity; (5) Primary antibody incubation and secondary antibody incubation; (6) DAB chromogenic reaction nuclear counterstaining, dehydration and mounting. Finally, staining of the tissues were visualized with a microscope, and images were acquired and analyzed. H-score: The depth and quantity of positivity was scored by Quant Center -an analysis software matched with a 3D scanner, which can only quantify the brown–yellow color of DAB. The larger the value, the stronger is the positivity. The assay refers to the methods in previously reported research (Yin et al., 2021).

ATP Assays

ATP was analized by enhanced ATP assay kit (Beyotime, S0027). The assay was performed ref to the protocols reported by previous studies (Kuang et al., 2021).

Protein isolation, Western blotting and antibodies

Proteins were extracted from cells with RIPA and protease inhibitor: phenylmethanesulfonyl fluoride (Beyotime, Haimen, China) and cocktail (MedChemExpress, Shanghai, China). The BCA method was used to measure protein concentrations. Western blotting was performed as previously described (Xie et al., 2018). The antibodies used were as follows: anti-OPA1 (1:1000, ab157457; Abcam); anti-cytochrome C (1:5000, ab133504; Abcam); and anti-GAPDH (1:1000, AG019; Beyotime).

Cell culture and transfection

Human LUAD cell lines (NCI-H1299 and NCI-A549) were bought from the Shanghai Institute of Biochemistry and Cell Biology. Cells were cultured in DMEM (Gibco, Carlsbad, MA) supplemented with 10% FBS (Gibco). These cells were incubated at 37 °C with 5% CO₂ in a humidified incubator.

The shRNAs specific for OPA1 were synthesized (RiboBio, Guangzhou, China; target sequence 1#: cgGGAGTTTGATCTTACCAAA and target sequence 2#: ccGGACCTTAGTGAATATAAA), and the overexpression plasmid targeting OPA1 was custom-designed (GENOME, Nanjing, China). The plasmid DNA or shRNA was transiently transfected into cells with Lipofectamine 2000 reagent (Invitrogen, Shanghai, China).

Cell proliferation assay

Cell proliferation was analyzed by CCK8 (Dojindo, Japan) and colony formation assays per the manufacturer’s instructions. A total of 4 ×10³ cells were seeded in 96-well plates (Corning Incorporated, Corning, New York, USA), and the culture medium was replaced with 10 µl CCK8 solution mixed with 100 µl RPMI 1640 or DMEM every 1 day and then incubated at 37 °C. The optical density value was measured using a Multimode microplate reader and obtained from five technical replicates (450 nm; Tecan, Mechelen, Belgium). 1 ×10³ cells were seeded in 6-well per plate (Corning Incorporated, Corning, New York, USA) and maintained for 14 days. Cell colony formation abilities were fixed with methanol and then stained with crystal violet (Beyotime, Haimen, China) for 30 min. All wells were photographed and counted. Each assay was performed at least three times. To evaluate the effects of mitochondrial target treatment on cell proliferation, cells were treated with N-(1,5-dimethyl-3-oxo-2-phenyl-2,3-dihydro-1H-pyrazol-4-YL)-3 methyl-1-PH+ (MYLS22, 50 µM) or Oligomycin A (10 µM) when the complete DMEM medium was added. Then, CCK8 and colony formation assays were performed to estimate the cell proliferation ratio of different treated cells.

Transwell assays

The migration and invasion capacity of lung cancer cells were investigated using Transwell plates (Corning Incorporated, Corning, NY, USA). Culture medium containing 10% FBS was added to the bottom compartment of the chamber as a chemoattractant. A total of 2 × 10⁴ cells/100 µl of lung cancer cells were planted on the upper chamber in serum free solution. Cells that adhered to the transwells were fixed with methanol for 15 min and then stained with crystal violet (Beyotime, Haimen, China), and quantified from five averaged fields via a Q-fired cooled CCD camera attached to an Olympus microscope. The assay was repeated three times in duplicate. To evaluate the effects of mitochondrial target treatment on cell migration, cells were treated with MYLS22 (50 µM) or Oligomycin A (10 µM) when migration was induced. Then, transwell based assays were performed to estimate the cell migration ability of different treated cells.

Mitochondrial membrane potential (Δψm)

Cells of different groups were incubated in tetramethylrhodamine methyl ester (TMRM; Thermo her Scientific, T-668) (50 nM). We obtained the images with LSM800 confocal microscope, equipped with a 60  × 1.40 NA oil immersion objective. Then the cell fluorescence intensity was analyzed with software ImageJ. Mitochondria with normal cells fluoresce bright red staining; however, red staining becomes diffuse or very light when mitochondrial membrane potential (Δψm) dissipates.

Mitochondrial Oxygen consumption rate (OCR) determination

Mitochondria oxygen consumption rate (OCR) was analyzed by using the Seahorse Extracellular Flux (XF96) Analyzer in group of both sh-OPA1 and control cells (Seahorse Bioscience Inc, USA). NCI-H1299 cells were plated at 4 × 10⁵ cells per well in the Seahorse V7 microplate for 24 h in the incubator. Then, the above cells were incubated in XF medium without NaHCO3 and FBS for 1 h/37 °C without CO2 input. OCR was calculated as percentage of the OCR value before the treatment of tested agents (pMoles/minute). The respiratory capacity was analyzed after the treatment with oligomycin and FCCP. OCR reduction after antimycin A, an inhibitor of ATP synthase, treatment represents ATP turnover under specific condition.

Metabolomics analysis

1 × 10⁷ cells in different groups (group of Control and shRNA-OPA1, n = 3/group) were seeded in 10 cm culture dishes. Then the cells were harvested and kept in liquid nitrogen. The metabolites were extracted and subjected to analysis of central carbon metabolism (Zhang et al., 2021).

Analysis of cell cycle progression and apoptosis using flow cytometry

Cells were harvested using 0.25% trypsin and washed with PBS. After resuspending in PBS and counting, 2 ×10⁵ cells were isolated and centrifuged to obtain a cell precipitate. Annexin V-FITC and propidium iodide (PI) staining solution were added with low-speed shaking and keep in the dark for 10 min.

Cells were re-suspended in 70% cold ethanol and 4 ° C overnight then centrifuged and washed in PBS. 150 µl RNase A and PI staining solution were then added to each tube and were then incubated at 4 °C. 20 min later, the results were analyzed by using flow cytometry and the ratio of cell among the G1, S and G2/M phases were analyzed (Zuo et al., 2014).

Quantitative real-time RT-PCR (qRT-PCR)

The cDNA obtained by reverse transcription was amplified on a QuantStudio™ 7 Flex Real-Time PCR System. The expressions of target genes were determined using the 2^−ΔCt method relative to GAPDH. The sequence of primers are presented in Table S1.

RNA sequencing and analysis

OPA1 knockdown and control LUAD cells were conducted using RNA-seq as described previously (Wang et al., 2018). The mRNA was extracted from total RNA using oligo (dT) beads to obtain the mRNA-seq library according to a standard protocol (Illumina, San Diego, USA) and then sequenced using an Illumina HiSeq X Ten system. The statistical power of this experimental design, was calculated in Source Package RNASeqPower in R (version 4.2) with the calculated power of 0.89. Three biological replicates were involved in both of the two groups for the RNA-seq analysis. 6G, 294 × depth raw data per sample were filtered and processed by Q30 and aligned to the mouse reference genome GRCh37/Gencode v19 with STAR software (v2.5.3a). Then, assembly and quantification of the transcripts were accomplished with RESM software (Version 1.3.3) guided by the Ensembl gtf gene annotation file (http://www.ensembl.org/). Read counts were used for the measurements of the relative abundance of the transcripts. The differential expressed genes (DEG) between groups were identified by the Wilcoxon rank sum test (Liu et al., 2021) (P < 0.05 and —log2[fold change]—>1). DAVID was further used to conduct an enrichment analysis of the DEGs (version 6.8; https://david.ncifcrf.gov/). The top 5 GO terms of the biological processes were examined. GSEA was performed using the R package cluster Profiler as previously reported (Chen et al., 2020).

Single-cell RNA sequencing (scRNA-seq) statistical analysis

The 10× Genomics scRNA-seq matrix of early stage LUAD tumor from primary lung tissues (tLUNG) were obtained from the previous study (Kim et al., 2020). The filtered and batched digital gene expression matrix (UMI counts per gene per cell) of 11 LUAD tLUNG samples was imported in R version 4.1 using Seurat v4.1.0. In our study, we used the t-Distributed Stochastic Neighbor Embedding (tSNE) algorithm for visual processing of this dimension-reduced dataset. Besides, main cell types were identified by scoring canonical cell type markers across clusters (Bischoff et al., 2021).

After identifying cell types in LUAD tLUNG samples, cell–cell communication was analyzed by implementing the Cell Chat R package v1.4.0 and used the Cell Chat website https://github.com/sqjin/CellChat. A new CellChat object was created from the merged Seurat object. The dataset of human CellChatDB was set as referencing database. Next, the communication probability was computed using the function computeCommunProb with the default parameters. After that, the cell–cell communication was inferred and the cell–cell communication network was aggregated with default parameters. The number of interactions was visualized to show the aggregated cell–cell communication network and signaling sent from each cell cluster.

Statistical analysis

We provided our raw numeric data for review and publication so the statistical analysis performed in our study can be checked (as shown in Tables S2–S7). All statistical analyses were performed by using R software. For experimental data, the Student’s t test or one-way ANOVA was used to assess differences in treatment groups using GraphPad Prism software (version 6.01). The results are expressed as the mean ± s.e.m. of three independent experiments. The significance threshold was set as a P value <0.05.

Identification of the clinical and prognostic value of mitochondrial dynamics-related genes in NSCLC from TCGA by bioinformatic analysis

To systematically evaluate the mechanism of mitochondria dynamics in NSCLC development, we performed the differential expression analysis of the five classical mitochondrial dynamics-related genes between the adjacent and tumor samples in the TCGA database (Fig. 1A). We identified that OPA1, MFN1, and DNM1L were significantly increased in both NSCLC tumor tissues (P < 0.05). However, Fis1 was decreased in tumors which implying that mitochondrial fusion is more active in tumors than mitochondrial fission. Fis1 is an important recruiting protein on the OMM. In prokaryotic cells, Fis1 assembles with DNM1L on the OMM to form a DNM1L ring and initiates mitochondrial fission by ring contraction and scission. Interestingly, in some cases, Fis1 is not necessary for mitochondrial fission, and the cells with Fis1 knockout can replicate and proliferate without obvious changes in mitochondrial morphology. In our study, we found that DNM1L was highly expressed in tumor tissues, while Fis1 was low expressed. Thus, we speculated that DNM1L may control fission through the regulation of other structural proteins than Fis1. To further explore the mechanism of the transcriptional activation of fusion related genes, we used a combined analysis with genomic data and found OPA1 and MFN1 have higher genomic alteration rate of 17% and 20%, with copy number variation (CNV) as the most common mutation type (Fig. 1B). The alteration type of the two genes is shown in Fig. 1C, mainly consisting of gain and amplification. Further analysis found the significant correlation of genomic amplification and gene expression of both OPA1 and MFN1 (Fig. 1D). It is worth noting that NSCLC patients with OPA1 and MFN1 gain or amplification tended to have co-occurrence implying their synergistic roles in mitochondrial dynamics. Immunohistochemical results also revealed that the protein expression of OPA1 was significantly increased in LUAD tumor tissues (Fig. 1E). In addition, Kaplan–Meier curves revealed that higher expression of OPA1 was significantly correlated with poor survival in LUAD patients, while other proteins does not (Fig. 1F). These results suggest that OPA1 plays a stronger role in tumorigenesis than other mitochondrial fusion and fission related genes, possibly because it has other functions independent of mitochondrial dynamics.

Loss of OPA1 induced mitochondrial dysfunction and metabolic reprogramming

We further examined transcriptional datasets from shOPA1 and control cells, and the PCA is shown in Fig. S1A. With a filter P value <0.05 and —log2Fold-change—>1, significantly different genes consisting of 1,156 down-regulated genes and 1,470 up-regulated genes in LUAD are shown in a volcano plot (Fig. 2A). Clustering analysis is shown in the heatmap in Fig. 2B. Specifically, we performed a GO analysis of these genes and revealed that they were involved in many signaling pathways, including mitochondrial function, ATP binding, ATP-dependent DNA helicase activity and membrane pathways (Fig. 2C), which provided clues about mitochondria. KEGG & COG_ONTOLOGY analysis also revealed dysfunction of mitochondrial metabolism (Figs. S1B & S1C). Furthermore, GSEA was performed to verify that the function of mitochondria and ATP production were inhibited in ShOPA1 cells (Figs. 2D and 2E). Furthermore, we conducted electron microscopy, and the representative images were shown in Fig. 2F, we identified the significantly reduced length/width ratio of mitochondrial in shOPA1 cells compared to the cell in control groups. To confirm the effect of OPA1 deficiency on ATP levels, we measured the production of ATP in shOPA1 and control cells, and the results showed that the ATP concentration of the OPA1 deficiency group was lower than that of the control group (Fig. 2G). These data inspired us to explore more about the aberrant metabolism caused by OPA1 deficiency.

OPA1 modulates the TCA metabolic process in LUAD cells

The mitochondrial metabolism of central carbon metabolism and oxygen consumption were further measured in NCI-H1299 cell lines. The obtained score plot of PCA analysis enabled a clear separation between the control group and the shOPA1 group (Fig. 3A). In addition, we identified 49 differentially abundant metabolites, consisting of 30 down-regulated and 19 up-regulated metabolites, based on a volcano plot (Fig. S1D). Variable importance in projection showed that significantly different metabolites were mainly involved in the processes of TCA and glycolysis (Fig. 3B). The relationship between the samples and metabolites is shown in Fig. S1E. When these significantly altered metabolites were subjected to unsupervised hierarchical clustering, a more defined pattern of metabolic alterations induced by exercise was observed. The results revealed significant alterations in nucleotide and amino acid metabolism, glycolysis, and the TCA cycle (Figs. 3C and 3D). A heatmap of all the metabolites is shown in Fig. S1F. In addition, the OCR was significantly lower in shOPA1 than in control cells, and proton leakage and ATP production were decreased (Figs. 3E and 3F). The above results demonstrated that the loss of OPA1 blocked the activity of mitochondria to balance metabolism.

OPA1⁺ tumor epithelial cells decreased its immune response to CD8⁺ T cell

To further explore the roles of OPA1 in tumorigenesis, we performed RNA-seq analysis and identified increased activity of immune response pathways (pathway of “immune response” and “MHC class II protein complex binding”) in OPA1 knocked down LUAD cells. GSEA analysis showed the higher activity of immune response and the mRNA expression of related genes, B2M, GBP2, CD24, HLA-F, CTSS, all increased in OPA1 knocked down LUAD cells. Moreover, we used sigle-cell sequencing data to further analysis the immune response of OPA1⁺ tumor epithelial cell to immune cells. All cells were divided into 19 clusters according to the t-SNE clustering algorithm and exhibited a higher expression of OPA1 in tumor epithelial cell (Figs. 4A–4B). Also, the proportion of OPA1 and MFN1 in various cells tend to be more distributed in tumor epithelial cells. Interestingly, it was found that OPA1 and MFN1 were strongly correlated in tumor epithelial cells (Fig. 4C). We further performed a GO analysis of the significantly different genes in OPA1⁺ and OPA1⁻ patients in different cell types (Fig. 4D). MKi67 expression is found within proliferating cells alone under general conditions and we found that MKi67 was higher expressed in OPA1⁺ cells (Fig. 4E). Go analysis of up-regulated and down-regulated genes showed that up-regulated genes mainly played a role in cell adhesion, while down-regulated genes mainly concentrated in immune response-related pathways (Figs. 4F–4G). By comparing the interaction network between OPA1⁺ and OPA1⁻ tumor tissues, we found that the communication of tumor epithelial with CD8⁺ T cells was significantly increased in OPA1⁻ tumor tissues (Fig. 4H). This suggests that OPA1⁻ tumor epithelial cells may regulate the immune responses strongly, while OPA1⁺ tumor epithelial cells are associated with immune escape (Fig. 4I).

OPA1 deficiency induced cell death via mitochondrial stress

Having confirmed the decreased activity of mitochondria, we further measured mitochondrial membrane potential. Mitochondrial membrane potential is a key mediator responsible for the activity of mitochondria (Zorov, Juhaszova & Sollott, 2006). As shown in Figs. 5A–5B, using TMRM as a fluorescent indicator for Δψm, we observed a 30% loss of mitochondrial membrane potential Δψm in shOPA1 compared with control cells. With the damaged mitochondrial membrane potential, we detected an increased amount of the cytochrome c protein in the whole cell (Fig. 5C). Cytochrome c is critical in the activation of programmed cell death pathway (Ow et al., 2008a). We further showed that shOPA1 treatment resulted in a significant increase in the percentage of early and late apoptotic cells (about 10.2%) compared to the control group through Annexin V-FITC (Fig. 5F). In addition, we detected the influence on the cell cycle of the damaged mitochondria. Flow cytometric analysis and RT-PCR verified that the cell cycle was blocked in G1 phase (Figs. 5D–5E), accelerating the process of apoptosis.

OPA1 overexpression shows enhanced proliferation and migration of LUAD cell lines

To evaluate the roles of OPA1 in LUAD, we knocked down OPA1 through virus infection (shRNAs) in the NCI-A549 cell line (Fig. 6A). We found that OPA1 knockdown exhibited significantly reduced cell growth rate (Fig. 6B) and resulted in fewer colonies (Fig. 6C). Moreover, transwell assays showed that knockdown of OPA1 inhibited cell migration compared with the negative control (Fig. 6D). Similarly, inhibition of cell growth, cloning and migration were observed in shOPA1 NCI-H1299 cells (Figs. 6E–6H). In addition, we overexpressed OPA1 in NCI-A549 cells (Fig. 6I), which significantly increased the cell growth rate (Figs. 6J) and colony numbers (6K). Additionally, transwell assays showed that OPA1 accelerated the speed of cell migration compared with the negative control (Fig. 6L).

We further tested whether MYLS22, the first in class safe and specific Opa1 inhibitor discovered by Stéphanie Herkenne et al. (Zamberlan et al., 2022), and Oligomycin A, a mitochondrial F0F1-ATPase inhibitor, could recapitulate the effects of genetic OPA1 silencing on breast cancer cell phenotype. A non-toxic dose of MYLS22 inhibited lung cancer cells migration and proliferation (Fig. S3). These results indicate that genetic or pharmacological OPA1 inhibition might reduce lung cancer cell proliferation, migration and invasiveness.