Identification of PRODH as a mitochondria- and angiogenesis-related biomarker for lung adenocarcinoma

Introduction

Lung cancer is the leading cause of cancer death, being attributable for an estimated 1.8 million deaths (18% of all deaths) in 2020 worldwide (1). Metastasis is the leading cause of cancer-related death, with more than 70% of patients exhibiting local or distant metastases at the time of initial diagnosis. Non-small cell lung cancer (NSCLC), which accounts for approximately 80–85% of lung cancers, is the most common pathological type and often and early leads to systemic metastasis, the mechanism of which is not yet fully understood (2). Since lung adenocarcinoma (LUAD) is prone to metastasis at an early stage, and since two-thirds of patients with LUAD are already at an advanced stage (stage IIIB/IV) at the time of diagnosis, their prognosis is poor, with an average 5-year survival rate of less than 20% (3).

Mitochondria are responsible for the bulk of cellular adenosine triphosphate (ATP) production through the process of oxidative phosphorylation. Mitochondrial membrane potential occurs in mitochondria establishes, which is known as the “powerhouse” of the cell (4). Therefore, being an essential organelle in cellular metabolism and cell death, mitochondria are a promising target for developing anticancer therapy (5). Indeed, defects in normal mitochondrial function are associated with a variety of human malignancies. Mitochondria have been suggested as a critical point of a key metabolic switch in normal cells in acquiring a malignant phenotype (6).

Angiogenesis is an integral part of tumor development and plays a key role in tumor growth and metastasis (7). In the 1970s, Folkman proposed that tumor growth and metastasis depend on angiogenesis and that inhibition of angiogenesis could be a therapeutic strategy for tumor treatment (8), and it has since been confirmed that the development of lung cancer relies on angiogenesis. In recent years, targeting angiogenic genes has become a research hotspot as a potential radiation-related treatment of lung cancer (9). We present this article in accordance with the REMARK reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2109/rc).

Methods

Data collection

We obtained the messenger RNA (mRNA) expression profiles of normal and LUAD tissues from The Cancer Genome Atlas (TCGA) database to screen for differentially expressed genes (DEGs). This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Differential expression analysis

We used a Venn diagram to screen for the DEGs of LUAD [log fold change (FC) >1.5] and mitochondrial-related genes that are also associated with angiogenesis. We also validated the expression levels of these selected genes using the GSE27262 dataset from the Gene Expression Omnibus (GEO). All the data analysis was performed with R version 4.2.2 software (The R Foundation of Statistical Computing).

Analysis of DEGs

Based on TCGA database, we screened for the hub genes downstream of proline dehydrogenase (PRODH). The data were analyzed by using pair plots. According to the expression levels of PRODH in lung cancer samples in TCGA were split into groups of high and low PRODH expression based on the median PRODH expression score. The R packages “limma” and “ggplot” were used to conduct DEG analysis between these two groups, with an adjusted P value <0.05 and |logFC| >1 set as the thresholds of DEGs.

Cancer stage and PRODH correlation analysis

Clinical data were downloaded from TCGA, and survival analysis was performed on PRODH. Univariate Cox and multivariate Cox analyses were performed to analyze the risk of downstream DEGs after grouping was performed based on PRODH expression levels. Correlation analysis was also conducted between PRODH and factors such as age, gender, and clinical stage.

PRODH-associated protein-protein interaction (PPI) network

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) website (10) (https://string-db.org/) was used to acquire the proteins related to PRODH under the following parameters: minimum required interaction score, “high confidence (0.700)”; meaning of network edges, “evidence”; max number of interactors to show, “no more than 5 interactors” in the first shell; and active interaction sources, “Experiments, Text mining, Databases, Co-expression, Neighborhood, Gene Fusion, Co-occurrence”. After the results were obtained from the STRING online database, they were then imported into the Cytoscape version 3.9.1 to identify the critical nodes for the visualization of the molecular interaction networks. Based on our constructed PPI network, the essential genes were identified via the CytoHubba plugin.

The enrichment analysis of PRODH

Subsequently, we obtained the top 289 PRODH expression-related genes using Cytoscape. The “ClusterProfiler” function package in R software was used for Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the DEGs. A threshold of adjusted P value <0.05 was applied to filter the significantly enriched GO terms and KEGG pathways. Bioconductor package “org.Hs.eg.db” version 3.13 was used for GO enrichment while “Release100.0” from the KEGG database was used for KEGG enrichment.

Acquisition of somatic mutation data

We obtained the somatic mutation data from the publicly available TCGA database through the Genomic Data Commons (GDC) Data Portal (https://portal.GDC.cancer.gov/). We selected “Masked Somatic Mutation” data from the data files of four subtypes and processed the data using VarScan software. We prepared mutation annotation format (MAF) for somatic variants and implemented the “maftools” R package, which provides multiple analysis modules, to perform visualization processes. Additionally, we downloaded transcriptome profiles of all available LUAD samples for comparison with normal tissue using the high-throughput sequencing (HTSeq) fragments per kilobase per million mapped fragments (FPKM) workflow. We also obtained corresponding clinical information from the GDC portal, including clinical variables such as age, gender, tumor grade, and pathological stage.

Calculation of tumor mutational burden (TMB) scores and prognostic analysis

TMB was defined as the total number of somatic gene coding errors, base substitutions, and insertions or deletions detected per million bases. In our study, we calculated the variant number or exon length for each sample using a Perl script based on the Java 8 software. We divided LUAD samples into a high-TMB group and a low-TMB group based on the median value (4.421053). We then merged TMB data with corresponding survival information using the sample identification number. Kaplan-Meier analysis was used to compare the difference in survival rate between the high-TMB group and the low-TMB group. Additionally, we further evaluated the association between TMB level and clinical features.

DEGs and functional pathways analysis

According to the TMB level, we divided the LUAD sample transcriptome data into a high-TMB group and low-TMB group using R software and used the R “limma” package to detect the DEGs between the two groups with FC =2 and false discovery rate (FDR) <0.05. Additionally, we obtained a list of immune-related genes from the Immunology Database and Analysis Portal (ImmPort) and selected differentially expressed immune genes between the two groups using the “VennDiagram” package in R.

Tumor Immune Estimation Resource (TIMER) database and CIBERSORT algorithm

Based on the “SCNA” module from the TIMER database (11), we further evaluated the mutation types of central immune genes with immune infiltrates in LUAD. Known mutation types of 20 central genes are displayed in the lower right corner, and the distribution of each immune cell subset under each mutation status in LUAD is represented by a box plot. The differences between each category and normal infiltration levels were compared using the two-sided Wilcoxon rank sum test, with the P values being calculated. We additionally obtained the transcriptome profiles for the two groups of patients with LUAD and normalized them using the “limma” package in R. We then input the prepared data into subsequent analysis, evaluating the immunological irregularity of each sample using the CIBERSORT algorithm, which provides estimates of member cell type abundance in a mixed cell population using gene expression data. CIBERSORT is still based on a known reference set, providing a set of gene expression features for 22 subtypes of white blood cells. The “heatmap” package in R was used to display the distribution of two immune cell subsets. Differences in immune infiltrate abundance between the high-TMB group and low-TMB group were compared using the Wilcoxon rank sum test, and the “vioplot” package in R was used to output the P values.

Drug sensitivity analysis and molecular docking

To determine the sensitivity of drugs potentially targeting PRODH, we analyzed the correlations between the PRODH expression levels and drug sensitivity, including obtaining the half maximal inhibitory concentration (IC₅₀) to cisplatin from the Genomics of Drug Sensitivity in Cancer (GDSC) database (http://www.cancerrxgene.org/). Regarding molecular docking, we searched the PubChem database to determine the name, molecular weight and 3D structure of cisplatin and then downloaded the 3D structure corresponding to the PRODH gene from the RCSB Protein Data Bank (PDB) (http://www.rcsb.org/). We subsequently used the AutoDock Vina software (http://vina.scripps.edu/) to prepare ligands and proteins for molecular docking. PyMOL (Schrödinger Inc., New York, NY, USA) was used to visualize the results (Figure 1).

Figure 1 Drug sensitivity of cisplatin in the GDSC database and molecular docking targets of the interaction between PRODH and cisplatin. GDSC, Genomics of Drug Sensitivity in Cancer; PRODH, proline dehydrogenase.

Cell culture

Five lung AD cell lines (NCI-H1975, PC-9, NCI‐H1299, A549, and H2126) and a normal human bronchial epithelial (HBE) cell line were purchased from the IMMOCELL (Xiamen, China). A549 and NCI‐H1299 cells were cultured in RPMI 1640 medium (BDBIO, China). PC-9, NCI‐H1975, H2126, and HBE cells were cultured in Dulbecco’s Modified Eagle Medium (DMEM; BDBIO, China) supplemented with 10% fetal bovine serum (FBS; BDBIO, China), 100 U/mL of penicillin, and 100 mg/mL of streptomycin (Invitrogen, Thermo Fisher Scientific, Waltham, MA, USA) at 37 ℃ in a 5% CO₂ atmosphere.

RNA interference with small interfering RNA (siRNA)

PC-9 and H1975 cells were plated and cultured in growth media until the cell density reached 60%, which was followed by transfection with siRNA (HanBio Therapeutics, Shanghai, China). The sequences of the siRNAs were GCACCUACUUCUACGCCAATT UUGGCGUAGAAGUAGGUGCTT, CCAAAUGGCUGUGGAGCAATT UUGCUCCACAGCCAUUUGGTT, and GGAAGUUCAAUGUGGAGAATT UUCUCCACAUUGAACUUCCTT. At 72 hours posttransfection, cells were harvested for Western blot analysis.

Western blotting

Total cellular proteins were extracted using a total protein extraction kit (Beyotime Institute of Biotechnology, Haimen, China). Cell lysates were separated via 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis and transferred to a nitrocellulose membrane. The membranes were blocked with 5% nonfat milk, incubated with primary antibodies, and then incubated with species-specific secondary antibodies. The following antibodies were used at the indicated concentrations: PRODH (#ER1915-44; HUABIO, Hangzhou, China), beta-actin (#HA1006; HUABIO), E-cadherin (#EM0502; HUABIO), N-cadherin (#ET1607-37; HUABIO), Vimentin (#ET1610-39; HUABIO), and Snail (#ER1706-22; HUABIO).

Cell Counting Kit-8 (CCK-8) assay

We conducted inoculation at 1,000 cells per well in a 96-well plate. The plates were cultured at 37 ℃ in 5% CO₂ for 24, 48, and 72 h; additionally, 100 µL of CCK-8 solution was added to each well, and the cells were cultured for 2 h. Finally, the cell absorption value of the plates was detected with a microplate reader at 450 nm.

Cell migration assay

Transfected cells (5×10⁴ cells/well) suspended in serum-free DMEM were seeded into the upper chamber of Transwell inserts. The completed DMEM was placed into the lower chamber. After 24 h, any unmigrated cells were removed, and the indicated cells were fixed and stained with 0.5% crystal violet and then photographed and counted under microscopy.

Statistical analysis

R software version 4.2.2 was used for statistical analysis, with P values of <0.05 indicating a statistical significance.

Results

Identification of DEGs

A total of 13,404 DEGs were identified from TCGA dataset. The top 100 DEGs are displayed in the heatmap in Figure 2A, and the DEGs between the LUAD group and the control group are displayed in the volcano plot in Figure 2B. Mitochondrial- and angiogenesis-related datasets were selected from the gene set enrichment analysis (GSEA) database for visualization analysis to identify corresponding genes (Figure 2C,2D). After TCGA and GSEA database analysis, two DEGs, LCAT1 and PRODH, were selected via a Venn diagram (Figure 2E).

Figure 2 Mitochondrial- and angiogenesis-related genes in lung adenocarcinoma. (A,B) Differential genes for lung adenocarcinoma in TCGA. (C,D) Mitochondrial pathway and angiogenic pathway in the GSEA. (E) Intersection of mitochondrial-related genes, angiogenic-related genes, and differential genes for lung adenocarcinoma. FC, fold change; TCGA, The Cancer Genome Atlas; GSEA, gene set enrichment analysis.

The analysis and validation of DEGs

GSE27262 from the GEO database was used to detect the differential expression of two DEGs, and the results showed that PRODH was highly expressed and had a significant difference in lung cancer tissue. We used pair and scatter plots to assess PRODH expression levels (Figure 3A-3C) and then used the Kaplan-Meier plotter to analyze the predictive value of PRODH expression levels for the prognosis of patients with LUAD. The results showed that the 10-year survival rate of patients with high PRODH expression was higher than that of those with low expression (P<0.05) (Figure 3D).

Figure 3 The expression and the Kaplan-Meier survival curve of PRODH. (A) Differential expression of two differential genes in the GEO (***, P ≤0.001). (B,C) PRODH in lung adenocarcinoma. (D) The Kaplan-Meier survival curve of PRODH. PRODH, proline dehydrogenase; LCAT, lecithin-cholesterol acyltransferase; LUAD, lung adenocarcinoma; HR, hazard ratio; CI, confidence interval; GEO, Gene Expression Omnibus.

The relationship between PRODH expression levels and the clinicopathological characteristics of patients with LUAD

GSE27262 from the GEO database was used to investigate the relationship between PRODH expression levels and the clinicopathological characteristics of patients with LUAD. The univariate Cox regression analysis showed a significant correlation between gender and PRODH expression levels (Figure 4A). Furthermore, multivariate Cox regression analysis identified tumor stage [hazard ratio (HR) =1.54; P<0.001] as an independent prognostic factor for patients with LUAD (Figure 4B). These results suggest that the expression level of PRODH is closely related to the clinicopathological characteristics of patients with LUAD.

Figure 4 The association of PRODH expression with the clinicopathological features of patients. (A) UniCox analyses of clinical characteristics. (B) MultiCox analyses of clinical characteristics. ***, P≤0.001. PRODH, proline dehydrogenase; CI, confidence interval; AIC, Akaike information criterion.

GO and KEGG enrichment analysis

The identified DEGs were subjected to KEGG pathway enrichment and GO annotation analysis using a clustering analyzer to characterize their biological function. The GO and KEGG pathway results are shown in Figure 5. The GO analysis showed that these DEGs were enriched in biological processes, including icosanoid metabolic process, apical part of cell, and apical plasma membrane, while their molecular functions mainly included endopeptidase activity. Furthermore, the data from KEGG analysis revealed that out of the 22 DEGs, enrichment was primarily in Ras signaling pathway, arachidonic acid metabolism, and pancreatic secretion (Figure 5).

Figure 5 Enrichment analysis based on the content of PRODH. (A,B) Differential genes screened by PRODH content. (C-F) GO and KEGG enrichment analysis of the relevant differential genes. **, P<0.01; ***, P<0.001. PRODH, proline dehydrogenase; FC, fold change; BP, biological process; CC, cellular component; MF, molecular function; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

The PPI network of the DEGs and the clinical correlation analysis

We used the online STRING tool to create a PPI network and identified hub genes to investigate potential interactions among all identified DEGs. As shown in Figure 6, the network of DEGs was complex, with the top 14 hub genes being SFTPB, SFTPC, SFTPA1, SCGB1A1, SFTPD, NKX2-1, MUC1, MUC5B, MUC3A, MUC21, B3GNT8, SFTPA2, NAPSA, and SCGB3A2 (Figure 6A,6B). We performed univariate and multivariate analyses of these 14 hub genes, and the HRs were visualized via forest plots (Figure 7A,7B).

Figure 6 Analysis based on the content of PRODH. (A) The PPI. (B) The 14 hub genes. PRODH, proline dehydrogenase; PPI, protein-protein interaction.

Figure 7 Univariate and multivariate Cox regression analyses identified the top 12 genes with high scores analyzed via the CytoHubba plugin. (A) Univariate Cox regression analyses identified the genes. (B) Multivariate Cox regression analyses identified the genes. CI, confidence interval.

TMB analysis

We used the maftools algorithm to examine the mutations in the high-risk and low-risk groups and found that for most genes, the frequency of mutations was higher in the high-risk group than in the low-risk group (TP53: low-risk 38%, high-risk 54%; TTN: low-risk 38%, high-risk 50%; MUC16: low-risk 38%; high-risk 42%) (Figure 8A-8C). In addition, the difference in TMB between the high- and low-risk groups was significant (P<0.05) (Figure 8D). We further investigated the possible differences in survival between patients with high- and low-TMB and found that the overall survival (OS) was significantly longer in the high-TMB group than in the low-TMB group (P<0.05) (Figure 8E,8F).

Figure 8 The TMB analysis between the high- and low-risk groups in lung adenocarcinoma. (A) Immune-related functions of hub genes (*, P<0.1; **, P<0.01; ***, P<0.001). (B,C) Waterfall plot of the top 15 mutated genes in lung adenocarcinoma in the high- and low-risk groups. (D) Differences in TMB between the high- and low-risk groups in lung adenocarcinoma. (E,F) Survival curves for the high- and low-TMB groups in lung adenocarcinoma and combined TMB risk survival curves. IFN, interferon; HLA, human leukocyte antigen; APC, antigen presenting cell; CCR, chemokine receptor; MHC, major histocompatibility complex; TMB, tumor mutational burden.

Landscape of mutation profiles in LUAD

We downloaded somatic mutation profiles of 616 patients with LUAD from TCGA, including four types of data based on different processing software. We used the “maftools” package to visualize the results of variant data in VCF format. The mutation information of each gene in each sample is visualized via the waterfall plot in Figure 9A, with the various colors annotated at the bottom representing different mutation types. These mutations were further classified according to different classification categories, among which missense mutations accounted for the majority (Figure 9B); single-nucleotide polymorphisms appeared more frequently than did insertions or deletions (Figure 9B); and C>A was the most common single-nucleotide variant (SNV) in LUAD (Figure 9B). Additionally, we calculated the number of altered bases in each sample, and the mutation types of LUAD are displayed in different colors in the box plot of Figure 9B. Finally, we identified the top 10 mutated genes in LUAD, which included TP53 (50%), TTN (43%), MUC16 (41%), CSMD3 (39%), RYR2 (34%), LRP1B (32%), ZFHX4 (31%), USH2A (29%), and KRAS (26%) (Figure 9A). The co-occurrence and exclusivity relationships between mutated genes are shown in Figure 9C, where the green color represents co-occurrence and the red color represents a mutually exclusive relationship. Meanwhile, the gene cloud map in Figure 9C shows the mutation frequencies of the other genes.

Figure 9 Progress in the study of mutation profiles in lung adenocarcinoma samples. Mutation information for each gene in each sample is shown in a waterfall plot, where various colors are annotated at the bottom to represent different mutation types. The markers above the legend show the mutation burden. (A,B) Classification of mutation types according to different categories: missense mutations accounted for the majority of mutations, SNPs appeared more frequently than do insertions or deletions, and C>A was the most common SNVs. The tumor mutation load of a given sample is shown; the x-axes are the mutation status of each sample, mutation type, and the number of mutated genes. (C) The top 10 mutated genes in LUAD, and the concordance and exclusivity associations between mutated genes. TMB, tumor mutational burden; TNP, tri-nucleotide polymorphism; SNPs, single-nucleotide polymorphisms; ONP, oligo-nucleotide polymorphism; INS, insertion; DEL, deletion; SNVs, single-nucleotide variants; LUAD, lung adenocarcinoma.

TMB correlated with survival outcomes, pathological stage, and tumor grade

We calculated the mutation events per million bases as the TMB for 336 patients with LUAD and further divided them into two groups of high- and low-TMB levels using the median TMB as the cutoff value. In addition, a higher TMB level was found to be associated with age (P=0.003) (Figure 10A), and a higher tumor grade was associated with gender (P=0.033) (Figure 10B). However, there was no significant correlation of TMB with T, N, or M stage (Figure 10C-10E).

Figure 10 The association of TMB with the clinicopathological features of patients. (A) The relationship between TMB level and age. (B)The relationship between TMB level and gender. (C) The relationship between TMB level and T stage. (D) The relationship between TMB level and N stage. (E) The relationship between TMB level and stage. TMB, tumor mutational burden.

Differential abundance of immune cells in the high- and low-TMB groups

As it was demonstrated that the DEGs were involved in immune crosstalk and that the DEGs mutations were negatively correlated with immune infiltration, we sought to further compare the differential distribution of immune components between the high- and low-TMB groups. After filtering analysis using the “CIBERSORT” package for samples with P>0.05, a total of 535 samples were selected for immune cell analysis. The box plot in Figure 11A shows the specific proportions of 22 immune cells in each LUAD sample. In addition, the Wilcoxon rank sum tests showed that the infiltration levels of CD8⁺ T cells, activated memory CD4⁺ T cells, M0 macrophages, and M1 macrophages were higher in the high-TMB group than in the low-TMB group (Figure 11B).

Figure 11 The relationship between the gene mutations and immune infiltration in LUAD. (A) The correlations between risk score and immune infiltration of 22 immune cell types in patients with LUAD. The x-axis represents the 535 patients from The Cancer Genome Atlas. (B) Different immune cells in the high-risk and low-risk groups. Green represents the low mutation burden group and red represents the high mutation burden group. (C) Twenty DEGs with |FC| >1 in the Venn diagram. NK, natural killer; FC, fold change; LUAD, lung adenocarcinoma; DEGs, differentially expressed genes.

Comparison of the gene expression profiles between the high- and low-TMB groups

We used differential analysis to generate a list of 20 DEGs with |FC| >1 as displayed in the Venn diagram in Figure 11C. As TMB is associated with immune features or pathways in LUAD, we further identified the top 20 immune-related genes from the ImmPort database for further analysis. Additionally, we evaluated the potential relationship between these gene mutations and immune infiltration in the LUAD microenvironment. Different forms of mutations carried by these genes can often inhibit immune infiltration compared to samples with characteristic wild-type immune infiltrates, including CD8⁺ T cells, neutrophils, dendritic cells, macrophages, CD4⁺ T cells, and B cells.

Relationship between the CNV of immune genes and immune cell infiltration

The TIMER database was used to investigate the correlation between the CNVs of immune-related DEGs and immune cell infiltration in LUAD. When the 20 genes varied in arm-level gain, the infiltration of B cells, CD8⁺ cells, CD4⁺ T cells, macrophages, and neutrophils decreased significantly in LUAD (Figures 12-15).

Figure 12 TIMER analysis of ALB, APOD, APOH, CALCA, and CD1A. ·, P<0.05; *, P<0.1; **, P<0.01; ***, P<0.001. ALB, albumin; APOD, apolipoprotein D; APOH, apolipoprotein H; CALCA, calcitonin related polypeptide alpha; CD1A, CD1a molecule; TIMER, Tumor Immune Estimation Resource.

Figure 13 TIMER analysis of CTSE, CXCL17, FGFR3, GAL, and GREM1. ·, P<0.05; *, P<0.1; **, P<0.01; ***, P<0.001. CTSE, cathepsin E; CXCL17, C-X-C motif chemokine ligand 17; FGFR3, fibroblast growth factor receptor 3; GAL, galanin and GMAP prepropeptide; GREM1, gremlin 1; TIMER, Tumor Immune Estimation Resource.

Figure 14 TIMER analysis of NR0B1, NR0B2, PCSK1, SCGB3A1, and SCTR. ·, P<0.05; *, P<0.1; **, P<0.01; ***, P<0.001. NR0B1, nuclear receptor subfamily 0 group B member 1; NR0B2, nuclear receptor subfamily 0 group B member 2; PCSK1, proprotein convertase subtilisin/kexin type 1; SCGB3A1, secretoglobin family 3A member 1; SCTR, secretin receptor; TIMER, Tumor Immune Estimation Resource.

Figure 15 TIMER analysis of SFTPA1, SFTPA2, SFTPD, SPINK5, and UCN3. ·, P<0.05; *, P<0.1; **, P<0.01; ***, P<0.001. SFTPA1, surfactant protein A1; SFTPA2, surfactant protein A2; SFTPD, surfactant protein D; SPINK5, serine peptidase inhibitor kazal type 5; UCN3, urocortin 3; TIMER, Tumor Immune Estimation Resource.

Drug sensitivity and molecular docking of PRODH

In order to potentially inform clinical treatment with the relevant functions of the PRODH gene, the drug sensitivity of PRODH was calculated using the information from the GDSC drug sensitivity database. As cisplatin is already a well-established clinical class of drugs and as a large number of studies on its relevant effects have been published, we conducted computer molecular docking simulations of cisplatin with PRODH (Figure 1).

PRODH was highly expressed and critical to the malignant behaviors in LUAD cells

The results of Western blotting showed that PRODH was highly expressed in LUAD cells, especially in the PC-9 and NCI-H1975 cell lines (Figure 16A). Compared with that in the control siRNA-transfected cells, the expression of PRODH was dramatically decreased in the si-PRODH#2 transfected cells (Figure 16B).

Figure 16 PRODH was highly expressed and integral to the malignant behavior of LUAD. (A) The expression levels of PRODH in the normal HBE cell line and five LUAD cell lines (H2126, H1975, PC-9, A549, and H1299) were examined using Western blotting. (B) The expression of PRODH was knocked down using siRNA in PC-9 cells. (C) CCK-8 assay of cell proliferation after knocking down PRODH in PC-9 and H1975 cells (***, P<0.001 for the difference between si-RNA with si-NC according to analysis of variance with Dunnett correction for multiple comparisons). (D) Transwell assays were conducted to examine the effects of PRODH knockdown on lung adenocarcinoma cell migration. The Transwell chamber was placed under a microscope at ×10 magnification to count the number of cells, which were stained with crystal violet staining (***, P<0.001 with the Student t-test). (E) Western blot shows the expression levels of E-cadherin, N-cadherin, Vimentin, and Snail after transfection with PRODH siRNA. PRODH, proline dehydrogenase; HBE, human bronchial epithelial; NC, negative control; OD, optical density; LUAD, lung adenocarcinoma; siRNA, small interfering RNA; CCK-8, Cell Counting Kit 8.

Cell proliferation was measured using CCK-8 assay after PRODH knockdown, which showed that silencing PRODH significantly reduced the proliferation ability of cells (Figure 16C). Next, we used Transwell assay to investigate the effect of PRODH on the invasion of LUAD cells, which showed that transfected LUAD cells exhibited significantly reduced invasive ability (Figure 16D). Western blotting was conducted, which indicated that PRODH knockdown influenced the levels of the proteins related to epithelial-mesenchymal transition, including E-cadherin, N-cadherin, Snail, and Vimentin (Figure 16E). These results suggest that PRODH is integral to the proliferation, migration, and invasion of LUAD.

Discussion

Lung cancer is one of the most commonly diagnosed cancers worldwide and a leading cause of cancer-related death. LUAD is the most common pathological type of lung cancer, and it often metastasizes through lymphatic and hematogenous routes (12).

Angiogenesis is a key process in the development and progression of tumors (13), including LUAD. Mitochondria plays an important role in regulating the energy metabolism in cancer cells (14), which is closely linked to angiogenesis. The increase in mitochondrial activity leads to the activation of various signaling pathways involved in angiogenesis. Conversely, angiogenic factors can also affect mitochondrial metabolism. Therefore, targeting both mitochondrial metabolism and angiogenesis yield significant therapeutic benefits for patients with LUAD.

In this comprehensive bioinformatics analysis, we identified genes related to vascular and mitochondrial functions in LUAD. PRODH is a pro-oxidant gene located in the inner mitochondrial membrane which directly transfers electrons to coenzyme Q1 (CoQ1) (15). Proline metabolism is related to ATP synthesis, protein and nucleotide synthesis, and redox homeostasis in tumor cells. The degradation of proline involves an oxidative step catalyzed by PRODH/proline oxidase (PRODH/POX) (16).

It has been reported that PRODH is involved in p53-induced reactive oxygen species (ROS)-dependent apoptotic response (17). One study demonstrated that PRODH is involved in regulating cyclooxygenase-2 (COX-2) (18). COX-2 is an enzyme involved in the biosynthesis of prostaglandins, and its expression is associated with poor prognosis in several malignant tumors (19). Reports indicate that high proline concentrations in cancer cells are associated with poor histological differentiation and an advanced clinical stage of malignancy (20,21). PRODH also plays a role in lung cancer (22) and has been identified as a potential target for developing anticancer drugs (17).

This study has several limitations that should be addressed. First, our microarray data were obtained from open access databases, which inevitably introduces systematic bias due to the sample heterogeneity across different studies. Second, there was a lack of direct evidence for the correlation in PRODH expression, human cancer prognosis, and immune cell infiltration. Third, a large amount of clinical samples are still needed to validate the prognostic effect of TMB and its potential relationship with immune infiltration. Further studies on genetic variations and large-scale clinical trials should be conducted in the future. Fourth, additional research on the related mechanism is still required to verify and explain the effect of PRODH on targeted therapy response. Finally, this study only examined the role of PRODH in tumors, and further investigations are needed to explore its roles in other diseases and biological processes.

Conclusions

This analysis study examined the role of the vascular- and mitochondrial-related gene, PRODH, in LUAD and statistically correlated PRODH expression with clinical prognosis, molecular characteristics, immune cell infiltration, immune-related genes, and TMB. In addition, possible mechanisms downstream of PRODH were clarified. Although further research is needed to validate these results, they suggest that PRODH may play a role in cancer prognosis and immune therapy response and thus may be a promising therapeutic target.

Acknowledgments

Funding: This work was supported by the Natural Science Research Program (key) of the Educational Commission in Anhui Province (No. KJ2021A0841) and Introducing Talents Natural Science Foundation of Yijishan Hospital of Wannan Medical College (No. YR202118).

Footnote

Reporting Checklist: The authors have completed the REMARK reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2109/rc

Data Sharing Statement: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2109/dss

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2109/prf

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-23-2109/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2019. CA Cancer J Clin 2019;69:7-34. [Crossref] [PubMed]
3. Lin JJ, Cardarella S, Lydon CA, et al. Five-Year Survival in EGFR-Mutant Metastatic Lung Adenocarcinoma Treated with EGFR-TKIs. J Thorac Oncol 2016;11:556-65. [Crossref] [PubMed]
4. Sharma P, Sampath H. Mitochondrial DNA Integrity: Role in Health and Disease. Cells 2019;8:100. [Crossref] [PubMed]
5. Porporato PE, Filigheddu N, Pedro JMB, et al. Mitochondrial metabolism and cancer. Cell Res 2018;28:265-80. [Crossref] [PubMed]
6. Tan YQ, Zhang X, Zhang S, et al. Mitochondria: The metabolic switch of cellular oncogenic transformation. Biochim Biophys Acta Rev Cancer 2021;1876:188534. [Crossref] [PubMed]
7. Viallard C, Larrivée B. Tumor angiogenesis and vascular normalization: alternative therapeutic targets. Angiogenesis 2017;20:409-26. [Crossref] [PubMed]
8. Proceedings Folkman J. Tumor angiogenesis factor. Cancer Res 1974;34:2109-13. [PubMed]
9. Li Y, Lin M, Wang S, et al. Novel Angiogenic Regulators and Anti-Angiogenesis Drugs Targeting Angiogenesis Signaling Pathways: Perspectives for Targeting Angiogenesis in Lung Cancer. Front Oncol 2022;12:842960. [Crossref] [PubMed]
10. Szklarczyk D, Gable AL, Nastou KC, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res 2021;49:D605-12. [Crossref] [PubMed]
11. Li T, Fan J, Wang B, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res 2017;77:e108-10. [Crossref] [PubMed]
12. Thai AA, Solomon BJ, Sequist LV, et al. Lung cancer. Lancet 2021;398:535-54. [Crossref] [PubMed]
13. Akbarian M, Bertassoni LE, Tayebi L. Biological aspects in controlling angiogenesis: current progress. Cell Mol Life Sci 2022;79:349. [Crossref] [PubMed]
14. Vasan K, Werner M, Chandel NS. Mitochondrial Metabolism as a Target for Cancer Therapy. Cell Metab 2020;32:341-52. [Crossref] [PubMed]
15. Hancock CN, Liu W, Alvord WG, et al. Co-regulation of mitochondrial respiration by proline dehydrogenase/oxidase and succinate. Amino Acids 2016;48:859-72. [Crossref] [PubMed]
16. Geng P, Qin W, Xu G. Proline metabolism in cancer. Amino Acids 2021;53:1769-77. [Crossref] [PubMed]
17. Burke L, Guterman I, Palacios Gallego R, et al. The Janus-like role of proline metabolism in cancer. Cell Death Discov 2020;6:104. [Crossref] [PubMed]
18. Liu Y, Borchert GL, Surazynski A, et al. Proline oxidase, a p53-induced gene, targets COX-2/PGE2 signaling to induce apoptosis and inhibit tumor growth in colorectal cancers. Oncogene 2008;27:6729-37. [Crossref] [PubMed]
19. Hashemi Goradel N, Najafi M, Salehi E, et al. Cyclooxygenase-2 in cancer: A review. J Cell Physiol 2019;234:5683-99. [Crossref] [PubMed]
20. Loayza-Puch F, Rooijers K, Buil LC, et al. Tumour-specific proline vulnerability uncovered by differential ribosome codon reading. Nature 2016;530:490-4. [Crossref] [PubMed]
21. Tołoczko-Iwaniuk N, Dziemiańczyk-Pakieła D, Celińska-Janowicz K, et al. Proline-Dependent Induction of Apoptosis in Oral Squamous Cell Carcinoma (OSCC)-The Effect of Celecoxib. Cancers (Basel) 2020;12:136. [Crossref] [PubMed]
22. Zhang L, Zhao X, Wang E, et al. PYCR1 promotes the malignant progression of lung cancer through the JAK-STAT3 signaling pathway via PRODH-dependent glutamine synthesize. Transl Oncol 2023;32:101667. [Crossref] [PubMed]