Mitochondrial Homeostasis–Related lncRNAs are Potential Biomarkers for Predicting Prognosis and Immune Response in Lung Adenocarcinoma

The prognosis of the most common histological subtype of lung cancer, lung adenocarcinoma (LUAD), is relatively poor. Mitochondrial homeostasis depends to a great extent on the coordination between mitophagy and mitochondrial biogenesis, the deregulation of which causes various human diseases, including cancer. There is accumulating evidence that long noncoding RNAs (lncRNAs) are critical in predicting the prognosis and immune response in carcinoma. Therefore, it is critical to discern lncRNAs related to mitochondrial homeostasis in LUAD patients. In this study, we identified mitochondrial homeostasis–related lncRNAs (MHRlncRNAs) by coexpression analysis. In order to construct a prognostic signature composed of three MHRlncRNAs, univariate and multivariate Cox regression analyses were performed. Kaplan–Meier analysis, stratification analysis, principal component analysis (PCA), receiver operating characteristic (ROC) curve, gene set enrichment analysis (GSEA), and nomogram were applied to evaluate and optimize the risk model. Subsequently, we identified the mitochondrial homeostasis–related lncRNA signature (MHLncSig) as an independent predictive factor of prognosis. Based on the LUAD subtypes regrouped by this risk model, we further investigated the underlying tumor microenvironment, tumor mutation burden, and immune landscape behind different risk groups. Likewise, individualized immunotherapeutic strategies and candidate compounds were screened to aim at different risk subtypes of LUAD patients. Finally, we validated the expression trends of lncRNAs included in the risk model using quantitative real-time polymerase chain reaction (qRT-PCR) assays. The established MHLncSig may be a promising tool for predicting the prognosis and guiding individualized treatment in LUAD.

Introduction

As the most commonly diagnosed pathological subtype of lung cancer, the incidence of lung adenocarcinoma (LUAD) is increasing globally every year. LUAD patients have poor prognosis which commonly results from late detection and individual treatment differences. Compared with the conventional treatments including surgery and chemotherapy, the development of multiple agents targeting driver gene mutations showed appreciable promise in LUAD treatment due to the advances of cancer genomics (Denisenko et al., 2018; Joseph et al., 2018). Unfortunately, secondary mutations in tumors usually contribute to targeted therapeutic resistance (Liu et al., 2018). In recent years, advance of cancer immunology makes immunotherapy a hotspot in the clinical treatment of LUAD. Classical immune checkpoint inhibitors (ICIs) including antiprogrammed cell death 1 (PD-1) and antiprogrammed cell death-ligand 1 (PD-L1) agents exert a persisting and powerful antitumor effect in LUAD patients (Forde et al., 2018). However, only a proportion of patients could benefit from immunotherapy because of the relatively low overall response rate of ICI (Li et al., 2018a). For the aforementioned reasons, the 5-year survival rate of LUAD remains inadequate (Wang et al., 2021; Zhu et al., 2022). Therefore, novel molecular biomarkers are urgently required to predict the prognosis and therapeutic response for LUAD patients.

Mitochondria are essential organelles that regulate ATP production and energy transformation via oxidative phosphorylation (OXPHOS) and the tricarboxylic acid cycle (TCA), which also modulate iron metabolism, Ca²⁺ signaling, innate immunity, and apoptotic cell death in mammalian cells (Zong et al., 2016; Pathak and Trebak, 2018). Mitochondria, generally defined as highly motile and plastic organelles, constantly undergo processes of fusion and fission and update through mitophagy and mitochondrial biogenesis to maintain homeostasis (Ma et al., 2020). In most cases, mitochondrial DNA (mtDNA) mutations, deleted, or damaged DNA replication induce the dysbiosis of mitochondrial homeostasis and consequent mitochondrial dysfunction (Butow and Avadhani, 2004; Wallace, 2012). Mitochondrial dysfunction caused by significant abnormalities in the mtDNA copy number is intimately associated with many diseases such as age-related pathologies, mtDNA depletion syndrome, and several carcinomas (Park et al., 2009; Imanishi et al., 2011; Greaves et al., 2012). The results from previous studies indicated that low copy numbers of mtDNA are widely observed in multiple cancers including colon, breast, hepatocellular carcinomas, prostate cancer, and astrocytoma (Horton et al., 1996; Lee et al., 2005; Petros et al., 2005; Tseng et al., 2006; Correia et al., 2011). Moreover, mtDNA depletion induced by experimental methods promotes aggressive phenotype in prostate and colorectal cancer cells (Moro et al., 2009; Guo et al., 2011).

Long noncoding RNAs (lncRNAs), a subclass of noncoding RNAs with a length of >200 nucleotides, have been confirmed to be involved in the tumorigenesis and progression of various tumors, including LUAD (Zheng et al., 2021). Recent studies have revealed that dysregulation of specific lncRNAs was inextricably associated with tumor proliferation, metastasis, and drug resistance in lung cancer (Mao et al., 2018; Wang et al., 2019a; Yang et al., 2020). Nonetheless, the specific function of lncRNAs in mitochondrial homeostasis remains to be clarified. Therefore, investigating the potential mechanism of mitochondrial homeostasis–related lncRNAs (MHRlncRNAs) in LUAD may be valuable for prognostic biomarkers.

In our study, we extracted the expression profiles of 1,499 mitochondrial homeostasis-related genes (MHRGs) from the publicly available dataset: The Cancer Genome Atlas (TCGA) dataset. Using Pearson’s correlation analysis, we identified 2,850 lncRNAs coexpressed with MHRGs. Next, mitochondrial homeostasis–related lncRNA signature (MHLncSig) was constructed to forecast the survival of LUAD patients utilizing differential expression analysis, univariate and multivariate Cox regression analyses. This personalized and robust prognostic signature is not only an independent indicator of overall survival, but also significantly related to tumor microenvironment, tumor mutation burden (TMB), immune infiltration, immunotherapeutic efficacy, and drug sensitivity.

Materials and Methods

Data Acquisition and Study Design

We downloaded the RNA transcriptome sequence data, corresponding clinical features, and mutation data for the patients with LUAD from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov/) database. The patients with missing survival information were removed from this study. The study flowchart is presented in Figure 1.

FIGURE 1

FIGURE 1. A concise flow diagram of this study.

Selection of Mitochondrial Homeostasis-Related Genes and lncRNAs

1,499 MHRGs with the relevance score >4 were obtained from the GeneCards (https://www.genecards.org/) database, including mitochondrial DNA-encoded genes and genomic DNA-encoded genes that regulate mitochondrial homeostasis via multiple pathways. We extracted the expression profile of MHRGs from the downloaded transcriptome data and performed differential expression analysis between the normal samples and LUAD samples. 294 differentially expressed MHRGs (|Log fold change (FC)| > 1 and false discovery rate (FDR) adjusted p < 0.05) were considered to be closely related to the pathogenesis and progression of LUAD. Next, we used Pearson’s correlation analysis and identified 2,850 MHRlncRNAs (|Pearson R| > 0.4 and p < 0.001). Further differential expression analysis based on 2,850 MHRlncRNAs was conducted to discern candidate lncRNAs for model building.

Establishment of the Risk Signature

We identified 131 differentially expressed MHRlncRNAs as candidate lncRNAs from the 2,850 MHRlncRNAs (|log FC| > 1 and FDR <0.05). The TCGA set was randomly separated into training (252 patients) and testing sets (252 patients), and no significant differences in the clinical properties was observed. Table 1 summarized the clinical features of these three sets. We used the training set to construct the prognostic signature The testing set and the entire TCGA set were utilized to validate the established signature. Univariate Cox regression analysis was performed to screen nine lncRNAs with significant prognostic value from 131 differentially expressed MHRlncRNAs. Multivariate Cox regression analysis was utilized to analyze and further confirm the prognostic significance from the previous steps, and three lncRNAs were retained during multiple computing. Thereafter, according to a linear combination of expression levels weighted with the regression coefficients calculated by multivariate Cox regression analysis, we established a risk formula that can calculate the risk score for each patient with LUAD

$MHLncSig\left(risk score\right)={\displaystyle \sum _{i=1}^{n}coefficient\left(\text{lncRNAi}\right)\ast expression level\left(\text{lncRNAi}\right) },$

where n indicates the amount of MHRlncRNAs included in the prognostic model, and the multivariate Cox regression analysis provided the coefficient for lncRNAi. According to the median risk score, we further distinguished the high- and low-risk subgroups of LUAD.

TABLE 1

TABLE 1. Clinical features of the three LUAD patients sets.

Function Enrichment Analysis

To discern the potential biological roles and metabolic pathways in differentially expressed MHRGs, the Gene Ontology (GO) categories and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were identified in the R program. We also used the Gene Set Enrichment Analysis (GSEA) (version 4.1.0) to investigate the biological processes and metabolic pathways involved in the different risk subgroups of LUAD stratified by MHLncSig. In terms of the reference file, c2.cp.kegg.v7.4.symbols.gmt was used and FDR <0.05 was considered significant.

Tumor Microenvironment and Tumor Mutation Burden (TMB)

R “estimate” package was utilized to calculate the stromal score, immune score, and ESTIMATE score for each LUAD patient, which represents the infiltration of both the stromal and immune cells in tumor tissues. We downloaded the somatic mutation data from the TCGA database and calculated the TMB score for each LUAD sample using a Perl script.

Immune Landscape and Immunotherapeutic Response Based on the MHLncSig

In order to calculate the relative abundance of 22 kinds of tumor-infiltrating immune cells (TICs) for each sample in the TCGA cohort, we used the CIBERSORT method. Subsequently, the Wilcoxon test was employed to verify the differentiation of 22 types of TICs between the low- and high-risk groups. Next, we converted the gene expression profile of tumor samples into immune function–related score and performed differential analysis of immune function between different risk subgroups using R “limma,” “GSVA,” “GSEABase,” “ggpubr,” and “reshape2” packages. Finally, we performed the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to access the response to immunotherapy.

Screening the Potential Compounds for Clinical Treatment of LUAD Risk Subgroups

We computed the IC50 of antitumor drugs that are commonly recommended for LUAD treatment based on gene expression profile of LUAD patients using the “pRRophetic” package in R software. The Wilcoxon test was adopted to evaluate the difference in the IC50 levels between distinct subgroups.

RNA Isolation and qRT-PCR

Patient samples comprising a cohort of 15 paired LUAD and adjacent normal tissues were collected between October 2021 and December 2021 from the Second Affiliated Hospital of Harbin Medical University. With the approval for experiments from the Ethics and Scientific Committees of the Second Affiliated Hospital of Harbin Medical University (Approval Number: KY2021-375), written informed consent was provided by all enrolled patients. The total RNA was extracted using the TRIzol reagent (Invitrogen, CA, United States). 1 µg total RNA and ReverTra Ace qPCR RT Master Mix (TOYOBO) were used for reverse transcriptase reaction. Next, 1 µl synthesized cDNA was used in PCR amplification. The levels of three lncRNAs were measured quantitatively by the SYBR Green Master Mix Kit (TOYOBO). GAPDH was selected as an internal reference. The relative expression was calculated based on the comparative Ct (2^−ΔΔCt) method, and Student’s t-test (two-tailed) was utilized to assess the significance of lncRNA expression differences in GraphPad Prism (version 8.0). The qRT-PCR primers in this study were as follows:

FENDRR

Forward: 5′-GCC​TCA​GAG​TGG​GCT​AGA​TT-3′

Reverse: 5′-TAA​CGA​TCC​CAC​CAA​CAC​CA-3′

AL590666.2

Forward: 5′-ACA​GAA​TGA​TCC​AGG​CAC​CA-3′

Reverse: 5′-AGG​ACA​AGA​TGG​ACG​CAG​AT-3′

AC090559.1

Forward: 5′-TGC​TAG​GCA​ATT​CTG​GAA​GC-3′

Reverse: 5′-TTG​CTG​TTG​CCA​CAA​AGT​GA-3′

GAPDH

Forward: 5′-CAT​GTT​CGT​CAT​GGG​TGT​GAA-3′

Reverse: 5′-GGC​ATG​GAC​TGT​GGT​CAT​GAG-3′

Statistical Analyses

We performed all bioinformatics analyses in R-version 4.1.1. Survival curves were plotted using Kaplan–Meier analysis to estimate the difference of survival for both the subgroups of LUAD in the R “survival” and “survminer” packages. Receiver operating characteristic (ROC) curves and the area under the curve (AUC) values were applied to evaluate the degree of robustness and accuracy of our MHLncSig. Principal component analysis (PCA) was utilized to reduce the dimension of high-dimensional data and visualize the differentiation between the risk subgroups in terms of entire genes expression in TCGA, 2,850 MHRlncRNAs, and the risk model of the 3 MHRlncRNAs. The R “limma” and “scatterplot3d” packages were employed to enable this process. In addition, we constructed a nomogram for predicting the overall survival of 1, 3, and 5 years, and a correction curve was applied to assess the consistency between the model prediction outcome and practical outcome by using the R “survival,” “survminer,” “rms,” and “regplot” packages, respectively. We set the statistically significant threshold as p-value <0.05 (*p < 0.05, **p < 0.01, and ***p < 0.001).

Results

Identification of Mitochondrial Homeostasis–Related lncRNAs in LUAD Patients

First, we downloaded the coding gene list associated with mitochondrial homeostasis from the GeneCards database. Based on the expression profile of LUAD and normal subjects in TCGA, differential expression analysis of 1,499 MHRGs was conducted to discern the genes that may serve a role in tumorigenesis. A heatmap and a volcano map were plotted to illustrate the 294 differentially expressed genes (205 upregulated genes and 89 downregulated genes in LUAD) (Figures 2A,B), and the detailed data are summarized in Supplementary Table S1. Furthermore, GO and KEGG analyses were utilized to explore the molecular functions and pathways behind these 294 genes. Figure 2C visualized the most highly significant cellular components, biological processes, and molecular functions. Moreover, biological pathways were mainly enriched in processes associated with mitochondrial synthesis and metabolism (Figures 2D,E).

FIGURE 2

FIGURE 2. Differential expression analysis and corresponding functional enrichment of the 1,499 MHRGs. (A) A heatmap of the 294 differentially expressed genes. (B) A volcano plot based on the 205 upregulated and 89 downregulated genes. (C) The most highly significant 10 results in the GO analysis. (D,E) The corresponding genes and the distribution of enriched biological pathways in the KEGG analysis.

Next, 2,850 lncRNAs coexpressed with the 294 differentially expressed genes were identified as MHRlncRNAs by Pearson’s correlation analysis. Furthermore, differential expression analysis based on the 2,850 lncRNAs was utilized to screen the candidate MHRlncRNAs for model construction. The result of differential analysis is summarized in Supplementary Table S2 and visualized in Figures 3A,B.

FIGURE 3

FIGURE 3. Construction and validation of the MHLncSig. (A) A heatmap of differential expression analysis based on 2,850 MHRlncRNAs. (B) A volcano plot based on the 109 upregulated MHRlncRNAs and 22 downregulated MHRlncRNAs. (C) A forest plot of nine MHRlncRNAs significantly related to survival in LUAD. (D,G,J) Different Kaplan–Meier curves in the two risk groups stratified by MHLncSig in the training set, testing set, and the entire TCGA set. (E,H,K) ROC curve of the MHLncSig at 3 years in the training set, testing set, and the entire TCGA set. (F,I,L) The expression patterns of MHLncSig, the risk grade distribution, and the survival status between the low- and high-risk groups in the training set, testing set, and the entire TCGA set.

Construction and Validation of the MHLncSig in LUAD Patients

To investigate the prognostic value of 131 candidate MHRlncRNAs, we performed univariate Cox analysis in the training set (252 patients) and identified nine MHRlncRNAs which are significantly associated with the survival of LUAD (Figure 3C). According to the multivariate Cox proportional hazards regression analysis, AL590666.2, FENDRR, and AC090559.1 were further established as components of the risk model. The details of multivariate Cox regression analysis are shown in Table 2. Finally, we designed a risk-score formula for LUAD patients’ survival prediction. The risk score formula is as follows: risk score = 0.063 × expression quantity of AL590666.2 + (−0.527) × expression quantity of FENDRR + (−0.158) × expression quantity of AC090559.1. AL590666.2 with a positive coefficient tended to be a detrimental factor. However, the other two lncRNAs including FENDRR and AC090559.1 tended to be protective factors. Next, grouping was conducted based on the median risk score of the training set; the Kaplan–Meier plot showed entirely different survival curves between different groups (p = 0.001, Figure 3D). The AUC value of the receiver operating characteristic (ROC) curve was 0.696, indicating that MHLncSig is equipped with quite accurate prediction performance for the prognosis (Figure 3E). In addition, the expression standards of MHLncSig for each patient, the risk grade distribution, and the survival status between the two groups are depicted in Figure 3F.

TABLE 2

TABLE 2. Multivariate Cox regression analysis of three prognostic lncRNAs.

Next, we validated our MHLncSig in the testing set (252 patients) and the entire TCGA set (504 patients). The survival analysis indicated that the patients in the low-risk group had better overall survival (OS) than that in the high-risk group (Figure 3G). The AUC value of MHLncSig in the testing set reached 0.695 (Figure 3H). Figure 3I shows the expression patterns of MHLncSig, the risk grade distribution, and the survival of LUAD patients. Notably, the entire TCGA set also manifested similar results as the aforementioned findings (Figures 3J–L).

Relationship Between MHLncSig and Clinical Traits

We next performed correlation analysis between the risk score calculated by MHLncSig and clinical traits including age, gender, T stage, N stage, M stage, and pathological stage. The results indicated that there were high correlations between MHLncSig and all these clinical features in the entire TCGA set (p < 0.05, Supplementary Figure S1). Furthermore, we evaluated the predictive performance of our MHLncSig in different subgroups classified by clinical traits. In the vast majority of subgroups of LUAD patients, the low-risk group continued to maintain supremacy in OS (Supplementary Figure S2). Moreover, univariate and multivariate Cox regression analyses were performed to investigate whether the MHLncSig had independent prognostic effect for LUAD. After correction for other clinical variables, we observed that the MHLncSig retained independent significance for the prediction of OS in the training set, the testing set, and the entire TCGA set (Table 3).

TABLE 3

TABLE 3. Univariate and multivariate Cox regression analysis of the MHLncSig and prognosis.

Evaluation of the Prognostic MHLncSig in Terms of PCA and GSEA Analyses

We conducted the PCA analysis to estimate the distributions of the two different risk groups based on the total gene expression profiles, 2,850 MHRlncRNAs, and the MHLncSig categorized by the expression profiles of the three risk lncRNAs (Figure 4A). Compared with the relatively scattered distributions of the two different risk groups based on the total gene expression profiles and 2,850 MHRlncRNAs, the results based on our prognostic signature indicated that the MHLncSig had excellent grouping ability to a certain extent. Furthermore, we performed GSEA enrichment analysis to reveal the underlying mechanisms and pathways behind the high-risk group characterized by dismal prognosis. A total of 37 pathways were significantly enriched in the high-risk group (FDR < 0.05, Supplementary Table S3). The pathways involving cell cycle, cysteine and methionine metabolism, glutathione metabolism, oxidative phosphorylation, proteasome, purine metabolism, and pyrimidine metabolism are visualized in Figure 4B.

FIGURE 4

FIGURE 4. PCA and GSEA analyses of the MHLncSig and construction of a prognostic nomogram. (A) PCA based on the total gene expression profiles, 2,850 MHRlncRNAs, and the MHLncSig categorized by three risk MHRlncRNAs. (B) GSEA analysis of the high-risk group calculated by MHLncSig. (C) A nomogram predicting 1-, 2-, and 3-year OS of LUAD cancer. (D) The calibration plot of the nomogram. (E) ROC curves of the prognostic nomogram at 1, 3, and 5 years.

Establishment and Performance Evaluation of the Prognostic Nomogram

We constructed a prognostic nomogram comprising the risk signature and clinical traits to predict the survival rate of 1, 2, and 3 years (Figure 4C). A total score was assigned to each LUAD patient by combining six individual scores in the nomogram, where a lower total point was related to a better outcome. Furthermore, we utilized ROC curve and calibration plot analyses to evaluate the validity of the nomogram (Figures 4D,E). The prediction curve was very close to the actual survival curve in the calibration plot and the AUCs of 1-, 2-, and 3-year curve were 0.728, 0.71, and 0.696, respectively, suggesting that the nomogram predicted the prognosis for LUAD patients well.

Evaluation of Tumor Microenvironment and TMB Using MHLncSig

Tumor microenvironment comprising stromal and immune components together with the secreted factors provide an immunosuppressive and protumoral environment for tumor development, which also correlates closely with immunotherapy positive response (Fridman et al., 2012; Gajewski et al., 2013). We made a comparison of the levels of stromal score, immune score, and ESTIMATE score between the two risk groups. The violin plots (Figure 5A) showed that all these three indicators were significantly decreased in the group with high risk. Further survival analysis suggested that lower levels of stromal score, immune score, and ESTIMATE score were significantly correlated with the worse prognosis in LUAD cohorts (Figure 5B). Moreover, we utilized R “maftools” package to analyze and visualize the mutational landscape of LUAD patients. The top 20 genes with the highest mutation frequency between the two risk subgroups are shown in Figures 5C,D. It is worth noting that the TMB score in the high-risk group exceeded that in the low-risk group, suggesting that MHLncSig had a high degree of TMB relevance (Figure 5E).

FIGURE 5

FIGURE 5. Evaluation of tumor microenvironment, TMB, and immune cell infiltration between different risk groups. (A) The differences of stromal score, immune score, and ESTIMATE score between the low- and high-risk groups. (B) The survival differences between the LUAD subgroups classified by stromal score, immune score, and ESTIMATE score. (C) The mutational landscape of the top 20 genes with the highest mutation frequency in the low-risk group. (D) The mutational landscape of the top 20 genes with the highest mutation frequency in the high-risk group. (E) The differences of TMB score between the low- and high-risk groups. (F) Estimation of immune cell infiltration in LUAD patients using the CIBERSORT algorithm.

Immune Landscape and Immunotherapeutic Response Based on the Prognostic Signature

To further explore the relationship between our prognostic signature and immune infiltration, we used the CIBERSORT algorithm to analyze the relative abundance of 22 TICs in each LUAD patient. Figure 5F showed that the contents of nearly half of these TICs were significantly different between the two different risk groups. Of these, the plasma cells, T-cell CD8, B-cell naive, macrophages M0, and T-cell CD4 memory activated were dramatically increased in the high-risk group, while monocytes, macrophages M2, T-cell CD4 memory resting, mast cells resting, and dendritic cells resting were significantly decreased in the high-risk group. Moreover, we observed a comprehensive suppression of immune functions in the high-risk group (Figure 6A). To elucidate the relationship between the MHLncSig and immune checkpoint genes, the expression of immune checkpoint genes were compared between the low-risk (n = 252) and high-risk (n = 252) groups. The results indicated that 40 types of immune checkpoint genes were significantly upregulated in the low-risk group (Figure 6B). Previous literatures have reported that the TIDE algorithm is a powerful tool to assess the immunotherapeutic response (Jiang et al., 2018). We discovered that LUAD patients in the high-risk group are more likely to benefit from immunotherapy, indicating that the MHLncSig may serve as a potential indicator for predicting immunotherapeutic response (Figure 6C).

FIGURE 6

FIGURE 6. Immune landscape and immunotherapeutic response based on the MHLncSig. (A) The difference of immune function between the low- and high-risk groups. (B) Expression levels of immune checkpoint genes in the low- and high-risk patients. (C) TIDE difference between the low- and high-risk groups.

Screening of Candidate Compounds Targeting MHLncSig and Model Validation in Human Surgical Resection Specimens

According to the pRRophetic algorithm, we calculated the IC50 of nine chemotherapeutic and targeted agents in the low- and high-risk patients, which are recommended by the National Comprehensive Cancer Network (NCCN) guidelines for LUAD treatment. Wilcoxon test analysis indicated that five of these agents (paclitaxel, docetaxel, erlotinib, pemetrexed, and carboplatin) had lower IC50 in high-risk patients (p < 0.05, Figures 7A–I), which suggested that LUAD patients with high risk were more sensitive to these five agents. Finally, based on human paired LUAD tissues obtained by surgery, we validated the differential expression of three risk lncRNAs included in the MHLncSig by qRT-PCR assays (Figures 7J–L). The differential analysis revealed that AL590666.2 was significantly upregulated in LUAD tissues, while FENDRR and AC090559.1 were remarkably downregulated in LUAD tissues.

FIGURE 7

FIGURE 7. Screening of candidate compounds and model validation in human LUAD specimens. (A–I) Sensitivity analysis of nine chemotherapeutic and molecular-targeted agents in the low- and high-risk patients. (J–L) Validation of the expression differences of three risk MHRlncRNAs in human LUAD tissues and adjacent normal tissues.

Discussion

In actual clinical work, pathologic staging is still the mainstream procedure to predict long-term survival and guide treatment modality for LUAD patients. However, a population with the same pathologic staging would also present with distinct clinical outcomes, indicating that the traditional staging system could not accommodate the individualized prediction and treatment. More recently, next-generation sequencing brings a revolutionary change in predicting the prognosis of cancer (Nair et al., 2018). A growing number of prognostic signatures on the basis of encoding genes and noncoding RNAs were established to forecast the survival outcome of cancer patients (Peng et al., 2021; Zhou et al., 2021). To the best of our knowledge, we are the first to identify the prognostic MHLncSig and comprehensively analyze the tumor immunity of this risk model in LUAD, which provides valuable novel options for the guidance of clinical personalized treatment.

Due to their pivotal functions in the resistance of cancer cells to regulated cell death (RCD) induced by treatment, mitochondria have drawn considerable attention for their potential value in developing novel anticancer agents. Actually, mitochondria affect immunosurveillance through both the malignant cell-extrinsic and malignant cell-intrinsic mechanisms. From one side, the mitochondria of malignant cells release many danger signals when RCD occurs, and these signals are paramount for activating the dendritic cells to promote antitumor immune responses. From another aspect, mitochondrial metabolism involves a number of functions related to antitumor immunity, including the establishment of protective immunological memory, inflammasome activation, and the differentiation of macrophage subsets (Porporato et al., 2018). It is worth noting that dysregulation of mitochondrial metabolism is found to be associated with metastatic cascade and worse prognosis in the nine types of tumors (Gaude and Frezza, 2016). Moreover, the dysbiosis of mitochondrial homeostasis is related to the excessive generation of reactive oxygen species (ROS) and consequent metastatic cascade (Fu et al., 2017). Recently, accumulated evidence has suggested that lncRNAs play important roles in tumorigenesis and tumor development. The dysfunctional lncRNAs may be presumed as key factors in the active phases of tumors (Guan et al., 2019). Zheng et al. (2015) discovered novel mechanisms of lncRNA HOTAIR in maintaining mitochondrial function in the HeLa cells. In addition, nuclear-encoded lncRNA MALAT1 was reported to regulate the metabolic reprogramming through the mitophagy pathway in hepatocellular carcinoma cells (Zhao et al., 2021). Nonetheless, prognostic biomarkers of LUAD and studies on the protumor mechanisms concerning MHRlncRNAs are still lacking. Therefore, we attempted to establish an independent signature for predicting the prognosis of LUAD based on MHRlncRNAs.

In our study, we identified 2,850 MHRlncRNAs based on the TCGA data. Subsequent differential analysis and univariate Cox analysis confirmed nine MHRlncRNAs with prognostic value, and three of these were utilized to construct a MHLncSig to predict the OS of LUAD patients.

Among them, FENDRR is commonly considered to be a tumor suppressor, which has been extensively studied in multiple cancer types, including cholangiocarcinoma, malignant melanoma, non–small cell lung cancer, hepatocellular carcinoma, and breast cancer (Li et al., 2018b; Zhang et al., 2018; Wang et al., 2019b; Qin et al., 2019; Chen et al., 2020). Nevertheless, specific mechanisms of two other lncRNAs (AL590666.2 and AC090559.1) during tumorigenesis and progression have never been previously reported and our research unmasked their underlying roles in LUAD.

TMB is defined typically as the total number of somatic coding mutations, which is associated with the production of neoantigens that induce antitumor immunity (Allgäuer et al., 2018). Recent studies have indicated that the higher TMB score was related to durable clinical benefit and improved objective response in tumor immunotherapy (Rizvi et al., 2015). Notably, the LUAD patients with high risk tend to possess higher levels of TMB, indicating that they are more likely to benefit from immunotherapy. This finding was subsequently confirmed by the TIDE algorithm. Moreover, another significant finding of our study suggested that the MHLncSig was significantly associated with the immune infiltration in LUAD, which further supported the fact that mitochondrial homeostasis plays a critical role in the tumor immune microenvironment.

There are still several limitations of the present study that need to be considered. First, it would be better to validate our prognostic model with several external datasets. In fact, we were unable to find an ideal Gene Expression Omnibus (GEO) database owing to the lack of complete expression profiles of the three risk MHRlncRNAs and detailed clinical information. Second, although we have initially screened a number of candidate compounds targeting MHLncSig in our study, the sensitivity of LUAD patients to specific immunotherapy agents requires further investigation due to a lack of immunotherapy drugs for LUAD in the R “pRRophetic” package.

Finally, functional experiments should be performed in vivo and in vitro to further corroborate our findings.

In summary, using multiple bioinformatics approaches and qRT-PCR experiments, we established and validated a 3-MHRlncRNA signature to independently predict the OS of LUAD patients. In addition, the MHLncSig shows advantages in terms of guiding individualized treatment of LUAD.

Conclusion

In conclusion, we successfully developed an accurate prognostic model of mitochondrial homeostasis, which was closely associated with tumor microenvironment, TMB, immune infiltration, and the response to immunotherapy in LUAD. Meanwhile, the findings in our study may provide clues for further elucidating the molecular mechanism of MHRlncRNAs in the tumorigenesis and LUAD progression.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics Statement

The studies involving human participants were reviewed and approved by the Ethics and Scientific Committees of the Second Affiliated Hospital of Harbin Medical University. The patients/participants provided their written informed consent to participate in this study.

Author Contributions

BP conceived the study, participated in the bioinformatics analysis, and wrote the manuscript. HL, CC, LW, and HX participated in the study design and performed experiments. HL, TL, RN, RX, TX, LY, and KW conducted the bioinformatics work. The manuscript was revised and approved by LZ and XL. All authors contributed to the article and approved the submitted version.

Conflict of Interest

The 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.

Acknowledgments

We sincerely thank professor Yong Zhang (Department of Pharmacology, College of Pharmacy, Harbin Medical University, China), for his precious suggestions contributed to the revision of the manuscript. We also thank all participants for their contribution to the study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.870302/full#supplementary-material

References

Allgäuer, M., Budczies, J., Christopoulos, P., Endris, V., Lier, A., Rempel, E., et al. (2018). Implementing Tumor Mutational Burden (TMB) Analysis in Routine Diagnostics-A Primer for Molecular Pathologists and Clinicians. Transl. Lung Cancer Res. 7 (6), 703–715. doi:10.21037/tlcr.2018.08.14

PubMed Abstract | CrossRef Full Text | Google Scholar

Butow, R. A., and Avadhani, N. G. (2004). Mitochondria and cell signalling. Mol. Cell. 14 (1), 1–15. doi:10.1016/s1097-2765(04)00179-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X.-e., Chen, P., Chen, S., Lu, J., Ma, T., Shi, G., et al. (2020). Long Non-coding RNA FENDRR Inhibits Migration and Invasion of Cutaneous Malignant Melanoma Cells. Biosci. Rep. 40 (3), BSR20191194. doi:10.1042/BSR20191194

PubMed Abstract | CrossRef Full Text | Google Scholar

Correia, R. L., Oba-Shinjo, S. M., Uno, M., Huang, N., and Marie, S. K. N. (2011). Mitochondrial DNA Depletion and its Correlation with TFAM, TFB1M, TFB2M and POLG in Human Diffusely Infiltrating Astrocytomas. Mitochondrion 11 (1), 48–53. doi:10.1016/j.mito.2010.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Denisenko, T. V., Budkevich, I. N., and Zhivotovsky, B. (2018). Cell Death-Based Treatment of Lung Adenocarcinoma. Cell. Death Dis. 9, 117. doi:10.1038/s41419-017-0063-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Forde, P. M., Chaft, J. E., Smith, K. N., Anagnostou, V., Cottrell, T. R., Hellmann, M. D., et al. (2018). Neoadjuvant PD-1 Blockade in Resectable Lung Cancer. N. Engl. J. Med. 378, 1976–1986. doi:10.1056/NEJMoa1716078

PubMed Abstract | CrossRef Full Text | Google Scholar

Fridman, W. H., Pagès, F., Sautès-Fridman, C., and Galon, J. (2012). The Immune Contexture in Human Tumours: Impact on Clinical Outcome. Nat. Rev. Cancer 12 (4), 298–306. doi:10.1038/nrc3245

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, L., Dong, Q., He, J., Wang, X., Xing, J., Wang, E., et al. (2017). SIRT4 Inhibits Malignancy Progression of NSCLCs, through Mitochondrial Dynamics Mediated by the ERK-Drp1 Pathway. Oncogene 36 (19), 2724–2736. doi:10.1038/onc.2016.425

PubMed Abstract | CrossRef Full Text | Google Scholar

Gajewski, T. F., Schreiber, H., and Fu, Y. X. (2013). Innate and Adaptive Immune Cells in the Tumor Microenvironment. Nat. Immunol. 14 (10), 1014–1022. doi:10.1038/ni.2703

PubMed Abstract | CrossRef Full Text | Google Scholar

Gaude, E., and Frezza, C. (2016). Tissue-specific and Convergent Metabolic Transformation of Cancer Correlates with Metastatic Potential and Patient Survival. Nat. Commun. 7, 13041. doi:10.1038/ncomms13041

PubMed Abstract | CrossRef Full Text | Google Scholar

Greaves, L. C., Reeve, A. K., Taylor, R. W., and Turnbull, D. M. (2012). Mitochondrial DNA and Disease. J. Pathol. 226 (2), 274–286. doi:10.1002/path.3028

PubMed Abstract | CrossRef Full Text | Google Scholar

Guan, H., Zhu, T., Wu, S., Liu, S., Liu, B., Wu, J., et al. (2019). Long Noncoding RNA LINC00673-V4 Promotes Aggressiveness of Lung Adenocarcinoma via Activating WNT/β-catenin Signaling. Proc. Natl. Acad. Sci. U.S.A. 116 (28), 14019–14028. doi:10.1073/pnas.1900997116

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, J., Zheng, L., Liu, W., Wang, X., Wang, Z., Wang, Z., et al. (2011). Frequent Truncating Mutation of TFAM Induces Mitochondrial DNA Depletion and Apoptotic Resistance in Microsatellite-Unstable Colorectal Cancer. Cancer Res. 71 (8), 2978–2987. doi:10.1158/0008-5472.CAN-10-3482

PubMed Abstract | CrossRef Full Text | Google Scholar

Horton, T. M., Petros, J. A., Heddi, A., Shoffner, J., Kaufman, A. E., Graham, S. D., et al. (1996). Novel Mitochondrial DNA Deletion Found in a Renal Cell Carcinoma. Genes. Chromosom. Cancer 15 (2), 95–101. doi:10.1002/(sici)1098-2264(199602)15:2<95:aid-gcc3>3.0.co;2-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Imanishi, H., Hattori, K., Wada, R., Ishikawa, K., Fukuda, S., Takenaga, K., et al. (2011). Mitochondrial DNA Mutations Regulate Metastasis of Human Breast Cancer Cells. PLoS One 6 (8), e23401. doi:10.1371/journal.pone.0023401

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Joseph, N. A., Chiou, S.-H., Lung, Z., Yang, C.-L., Lin, T.-Y., Chang, H.-W., et al. (2018). The Role of HGF-MET Pathway and CCDC66 cirRNA Expression in EGFR Resistance and Epithelial-To-Mesenchymal Transition of Lung Adenocarcinoma Cells. J. Hematol. Oncol. 11, 74. doi:10.1186/s13045-018-0557-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H.-C., Yin, P.-H., Lin, J.-C., Wu, C.-C., Chen, C.-Y., Wu, C.-W., et al. (2005). Mitochondrial Genome Instability and mtDNA Depletion in Human Cancers. Ann. N. Y. Acad. Sci. 1042, 109–122. doi:10.1196/annals.1338.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Shao, C., Shi, Y., and Han, W. (2018). Lessons Learned from the Blockade of Immune Checkpoints in Cancer Immunotherapy. J. Hematol. Oncol. 11, 31. doi:10.1186/s13045-018-0578-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., Zhang, W., Liu, P., Xu, Y., Tang, L., Chen, W., et al. (2018). Long Non-coding RNA FENDRR Inhibits Cell Proliferation and Is Associated with Good Prognosis in Breast Cancer. Ott Vol. 11, 1403–1412. doi:10.2147/OTT.S149511

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Q., Yu, S., Zhao, W., Qin, S., Chu, Q., and Wu, K. (2018). EGFR-TKIs Resistance via EGFR-independent Signaling Pathways. Mol. Cancer 17, 53. doi:10.1186/s12943-018-0793-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, K., Chen, G., Li, W., Kepp, O., Zhu, Y., and Chen, Q. (2020). Mitophagy, Mitochondrial Homeostasis, and Cell Fate. Front. Cell. Dev. Biol. 8, 467. doi:10.3389/fcell.2020.00467

PubMed Abstract | CrossRef Full Text | Google Scholar

Mao, C., Wang, X., Liu, Y., Wang, M., Yan, B., Jiang, Y., et al. (2018). A G3BP1-Interacting lncRNA Promotes Ferroptosis and Apoptosis in Cancer via Nuclear Sequestration of P53. Cancer Res. 78, 3454–3496. doi:10.1158/0008-5472.CAN-17-3454

PubMed Abstract | CrossRef Full Text | Google Scholar

Moro, L., Arbini, A. A., Yao, J. L., di Sant'Agnese, P. A., Marra, E., and Greco, M. (2009). Mitochondrial DNA Depletion in Prostate Epithelial Cells Promotes Anoikis Resistance and Invasion through Activation of PI3K/Akt2. Cell. Death Differ. 16 (4), 571–583. doi:10.1038/cdd.2008.178

PubMed Abstract | CrossRef Full Text | Google Scholar

Nair, M., Sandhu, S. S., and Sharma, A. K. (2018). Cancer Molecular Markers: A Guide to Cancer Detection and Management. Seminars Cancer Biol. 52 (Pt 1), 39–55. doi:10.1016/j.semcancer.2018.02.002

CrossRef Full Text | Google Scholar

Park, J. S., Sharma, L. K., Li, H., Xiang, R., Holstein, D., Wu, J., et al. (2009). A Heteroplasmic, Not Homoplasmic, Mitochondrial DNA Mutation Promotes Tumorigenesis via Alteration in Reactive Oxygen Species Generation and Apoptosis. Hum. Mol. Genet. 18 (9), 1578–1589. doi:10.1093/hmg/ddp069

PubMed Abstract | CrossRef Full Text | Google Scholar

Pathak, T., and Trebak, M. (2018). Mitochondrial Ca2+ Signaling. Pharmacol. Ther. 192, 112–123. doi:10.1016/j.pharmthera.2018.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Peng, B., Li, H., Na, R., Lu, T., Li, Y., Zhao, J., et al. (2021). Identification of a Novel Prognostic Signature of Genome Instability-Related LncRNAs in Early Stage Lung Adenocarcinoma. Front. Cell. Dev. Biol. 9, 706454. doi:10.3389/fcell.2021.706454

PubMed Abstract | CrossRef Full Text | Google Scholar

Petros, J. A., Baumann, A. K., Ruiz-Pesini, E., Amin, M. B., Sun, C. Q., Hall, J., et al. (2005). mtDNA Mutations Increase Tumorigenicity in Prostate Cancer. Proc. Natl. Acad. Sci. U.S.A. 102 (3), 719–724. doi:10.1073/pnas.0408894102

PubMed Abstract | CrossRef Full Text | Google Scholar

Porporato, P. E., Filigheddu, N., Pedro, J. M. B.-S., Kroemer, G., and Galluzzi, L. (2018). Mitochondrial Metabolism and Cancer. Cell. Res. 28 (3), 265–280. doi:10.1038/cr.2017.155

PubMed Abstract | CrossRef Full Text | Google Scholar

Qin, X., Lu, M., Zhou, Y., Li, G., and Liu, Z. (2019). LncRNA FENDRR Represses Proliferation, Migration and Invasion through Suppression of Survivin in Cholangiocarcinoma Cells. Cell. Cycle 18 (8), 889–897. doi:10.1080/15384101.2019.1598726

PubMed Abstract | CrossRef Full Text | Google Scholar

Rizvi, N. A., Hellmann, M. D., Snyder, A., Kvistborg, P., Makarov, V., Havel, J. J., et al. (2015). Mutational landscape determines sensitivity to PD-1 blockade in non–small cell lung cancer. Science 348 (6230), 124–128. doi:10.1126/science.aaa1348

PubMed Abstract | CrossRef Full Text | Google Scholar

Tseng, L.-M., Yin, P.-H., Chi, C.-W., Hsu, C.-Y., Wu, C.-W., Lee, L.-M., et al. (2006). Mitochondrial DNA Mutations and Mitochondrial DNA Depletion in Breast Cancer. Genes. Chromosom. Cancer 45 (7), 629–638. doi:10.1002/gcc.20326

PubMed Abstract | CrossRef Full Text | Google Scholar

Wallace, D. C. (2012). Mitochondria and Cancer. Nat. Rev. Cancer 12 (10), 685–698. doi:10.1038/nrc3365

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, B., Xian, J., Zang, J., Xiao, L., Li, Y., Sha, M., et al. (2019). Long Non-coding RNA FENDRR Inhibits Proliferation and Invasion of Hepatocellular Carcinoma by Down-Regulating Glypican-3 Expression. Biochem. Biophysical Res. Commun. 509 (1), 143–147. doi:10.1016/j.bbrc.2018.12.091

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., Mao, C., Ouyang, L., Liu, Y., Lai, W., Liu, N., et al. (2019). Long Noncoding RNA LINC00336 Inhibits Ferroptosis in Lung Cancer by Functioning as a Competing Endogenous RNA. Cell. Death Differ. 26, 2329–2343. doi:10.1038/s41418-019-0304-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Tan, H., Yu, T., Chen, X., Jing, F., and Shi, H. (2021). Potential Immune Biomarker Candidates and Immune Subtypes of Lung Adenocarcinoma for Developing mRNA Vaccines. Front. Immunol. 12, 755401. doi:10.3389/fimmu.2021.755401

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, Y., Tai, W., Lu, N., Li, T., Liu, Y., Wu, W., et al. (2020). lncRNA ZFAS1 Promotes Lung Fibroblast-To-Myofibroblast Transition and Ferroptosis via Functioning as a ceRNA through miR-150-5p/SLC38A1 axis. Aging 12, 9085–9102. doi:10.18632/aging.103176

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M. Y., Zhang, Z. L., Cui, H. X., Wang, R. K., and Fu, L. (2018). Long Non-coding RNA FENDRR Inhibits NSCLC Cell Growth and Aggressiveness by Sponging miR-761. Eur. Rev. Med. Pharmacol. Sci. 22 (23), 8324–8332. doi:10.26355/eurrev_201812_16530

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, Y., Zhou, L., Li, H., Sun, T., Wen, X., Li, X., et al. (2021). Nuclear-Encoded lncRNA MALAT1 Epigenetically Controls Metabolic Reprogramming in HCC Cells through the Mitophagy Pathway. Mol. Ther. - Nucleic Acids 23, 264–276. doi:10.1016/j.omtn.2020.09.040

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, J., Zhou, Z., Qiu, Y., Wang, M., Yu, H., Wu, Z., et al. (2021). A Prognostic Ferroptosis-Related lncRNAs Signature Associated with Immune Landscape and Radiotherapy Response in Glioma. Front. Cell. Dev. Biol. 9, 675555. doi:10.3389/fcell.2021.675555

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, P., Xiong, Q., Wu, Y., Chen, Y., Chen, Z., Fleming, J., et al. (2015). Quantitative Proteomics Analysis Reveals Novel Insights into Mechanisms of Action of Long Noncoding RNA Hox Transcript Antisense Intergenic RNA (HOTAIR) in HeLa Cells. Mol. Cell. Proteomics 14 (6), 1447–1463. doi:10.1074/mcp.M114.043984

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, P., Lu, Y., Zhang, Y., and Wang, L. (2021). Construction of an Immune-Related Six-lncRNA Signature to Predict the Outcomes, Immune Cell Infiltration, and Immunotherapy Response in Patients with Hepatocellular Carcinoma. Front. Oncol. 11, 661758. doi:10.3389/fonc.2021.661758

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, M., Zhang, L., Cui, H., Zhao, Q., Wang, H., Zhai, B., et al. (2022). Co-Mutation of FAT3 and LRP1B in Lung Adenocarcinoma Defines a Unique Subset Correlated with the Efficacy of Immunotherapy. Front. Immunol. 12, 800951. doi:10.3389/fimmu.2021.800951

PubMed Abstract | CrossRef Full Text | Google Scholar

Zong, W.-X., Rabinowitz, J. D., and White, E. (2016). Mitochondria and Cancer. Mol. Cell. 61 (5), 667–676. doi:10.1016/j.molcel.2016.02.011

PubMed Abstract | CrossRef Full Text | Google Scholar