A nicotinamide metabolism-related gene signature for predicting immunotherapy response and prognosis in lung adenocarcinoma patients

Data collection and preprocessing

The somatic mutation information and clinical phenotype data of LUAD were collected from the Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The RNA-seq expression profile was converted to FPKM format and log2-converted. The somatic gene mutation information was processed by mutect2 (Saba et al., 2020). After removing samples that lacked survival time or survival status, a total of 500 tumor samples with a survival time longer than 0 day were obtained.

The chip data of GSE31210 (Shahrajabian & Sun, 2023) was acquired from Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE31210) (Song et al., 2023). The probes were transformed to genes symbol based on the annotation information. A total of 226 tumor samples from GSE31210 were recruited for subsequent study after removing normal samples and those without clinical follow-up or overall survival (OS) data.

TCGA-LUAD and GSE31210 datasets served as the training dataset and independent validation dataset, respectively.

Consensus clustering analysis

Firstly, 42 NMRGs obtained from a previous study (Cui et al., 2023) were subjected to univariate Cox regression analysis (p < 0.05) to filter prognostically significant NMRGs in LUAD. Then, consensus clustering was performed to delineate molecular subtypes of TCGA-LUAD using ConsensusClusterPlus R package (the parameters were set as clusterAlg = “hc” and distance = “Spearman”) (Li et al., 2020). Sampling was repeated 500 times with 80% samples involved each time. The optimal cluster number (between 2 to 10) was determined according to the CDF. Finally, the prognosis between different subtypes in TCGA-LUAD cohort was evaluated and validated in GSE31210 dataset.

Somatic gene mutation landscape analysis

Somatic genes with mutation frequencies higher than 2 were selected. Then, genes with significantly high mutation frequencies within each molecular subtype were further identified using Fisher’s exact test, setting a threshold of p < 0.05. Next, the mutation landscape of the top 20 somatic genes in each molecular subtype was visually displayed using a waterfall diagram, generated with the maftools R package (Jia et al., 2024).

Identification of DEGs and GSEA

DEGs among different subtypes of TCGA-LUAD were selected using limma R package, with a screening threshold of |log2FoldChange (FC)|>log2(1.5) and p < 0.05. Following this, pathway enrichment analysis was conducted via GSEA R package to pinpoint significantly enriched pathways, using a threshold of p < 0.05 (Innis et al., 2021). The gene set “h.all.v2023.1.Hs.entrez.gmt” from the MsigDB (https://www.gsea-msigdb.org/gsea/msigdb/) served as the background gene sets (Li et al., 2023), encompassing the HALLMARK series of pathways.

Development and verification of a RiskScore model

Using univariate Cox regression analysis (p < 0.05), the prognostic genes were selected from the DEGs between the molecular subtypes in TCGA-LUAD. Subsequently, the model was refined by performing LASSO Cox regression analysis in the glmnet R package and 10-fold cross validation (Li et al., 2024). Then, the independent prognostic NMRGs were screened via stepwise multivariate regression analysis and applied to develop the RiskScore model as follow (Fan et al., 2024):

$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\sum \mathrm{\beta }\mathrm{i}\ast \mathrm{E}\mathrm{x}\mathrm{P}\mathrm{i}$$

βi indicates the coefficient of gene in Cox regression model, and ExPi represents gene expression value.

The RiskScore for all the samples in TCGA-LUAD and GSE31210 datasets was calculated, and LUAD patients were separated by the median RiskScore value into high- and low-risk groups. The ROC curve was constructed using timeROC R package (Zhang, 2022). The model robustness was validated using the GSE31210 dataset.

Evaluation of the TME

To investigate the association between the TME in LUAD and different risk groups, the scores of a total of 28 types of tumor infiltration immune cells (TIICs) were computed for each sample (Wang et al., 2020). Also, the scores of 10 types of immune cells were calculated by MCP-counter (Zhang et al., 2022). The TIMER online tool (http://cistrome.org/TIMER) was employed to calculate the scores of six immune cells (Cao et al., 2022). Moreover, ESTIMATE algorithm was utilized to output the ESTIMATEScore, StromalScore, ImmuneScore using estimated R package (Ke et al., 2021).

Immunotherapy response and drug sensitivity assessment

Firstly, the TIDE score was calculated based on the standardized transcriptome data in TIDE website (http://tide.dfci.harvard.edu/) (Zheng et al., 2024) to evaluate the immunotherapy response of LUAD patients, with higher TIDE scores indicating greater immune escape possibility and less immunotherapy benefit. Furthermore, the sensitivity of LUAD patients in TCGA-LUAD to 10 common chemotherapy drugs was compared based on the IC₅₀ value using pRRopetic R package.

Development and verification of a NAM metabolism-related nomogram

Firstly, independent prognostic factors were identified by subjecting clinical traits and RiskScore model to univariate and multivariate Cox regression analysis. Then, a nomogram integrating RiskScore with other clinicopathologic features was created to evaluate the 1-, 3-, and 5- year survival rate of patients with LUAD (Liu et al., 2024). The nomogram’s performance was assessed using a calibration curve and DCA. The predictive accuracy was further validated by calculating the C-index (Wang et al., 2022). Additionally, the diagnostic efficacy of the nomogram was evaluated through plotting receiver operating characteristic (ROC) curves using the timeROC R package (Lin et al., 2023).

Cell culture and transfection

Normal human lung epithelial cell line (BEAS-2B) and the LUAD cell line (A549) were procured from Wuhan Sunncell Biotechnology Co., Ltd (Wuhan, China). DMEM/F12 medium comprising 10% fetal bovine serum (FBS) and 1% penicillin-streptomycin was used in all cell culture under the incubation conditions of 37 °C, 5% CO₂, and saturated humidity.

Subsequently, the silencing of GJB3 was conducted by cell transfection. The small interfering (si) RNA specifically targeting GJB3 (si-GJB3) along with a negative control (si-NC) was obtained (Shanghai Sangon Biotechnology Co., Ltd., Shanghai, China). The sequence of si-GJB3 was AAGTTATGCAACTTTCGTTTTGG. Transfection of A549 cells was conducted with the use of Lipofectamine® 3000 (Invitrogen, Thermo Fisher Scientific) for 48 h, according to the instruction.

Quantitative real‑time PCR

Firstly, the total RNA of BEAS-2B and A549 cells obtained using TriZol reagent (Invitrogen, Waltham, MA, USA) was reverse-transcribed into cDNA with SuperScript IV reverse transcriptase (Invitrogen, Waltham, MA, USA). Next, to measure the mRNA expressions of independent prognostic NMRGs in BEAS-2B and A549 cells, quantitative real–time PCR (qRT–PCR) was conducted for amplification using SYBR® Premix Ex Taq™ II (Takara, Shanghai, China). See Table S1 for the primer sequences. The amplification was performed as follows: pre-denaturation at 95 °C for 1 minute (min), 40 cycles of 95 °C for 30 seconds (s), and annealing at 56 °C for 30 s, elongation at 72 °C for 40 s. The relative expressions of NMRGs were quantified by 2^−ΔΔCT method, applying GAPDH as the housekeeping gene (Song et al., 2023).

Wound-healing assay

The A549 cell migration was measured by wound-healing assay (Zhang et al., 2022). A549 cells at a concentration of 2 × 10⁵ cells/well were cultured in a 24-well plate and scratched using an aseptic pipette tip. After washing with phosphate buffer, the A549 cells were cultivated in serum-free medium. The BX53M upright metallographic microscope (Olympus, Tokyo, Japan) was applied to observe and capture the representative images at 0 and 48 h. Meanwhile, the wound closure (%) of A549 cells was counted.

Transwell assay

Transwell assay was conducted to evaluate the invasion of A549 cells (Zhou et al., 2021). Firstly, each Transwell chamber was added with 70 µL diluted Matrigel and placed at 37 °C for 2 h to coagulate the Matrigel. Then, the A549 cells (2 × 10⁵ cells/well) were suspended in 250 µL non-serum medium contained in the upper chamber and incubated, while 650 µL medium comprising 10% FBS was supplemented to the lower chamber. Following 24-h cultivation, 5% paraformaldehyde was used to fix the A549 cells invading the lower chamber for 30 min. Then, 0.1% crystal violet was employed to dye the fixed cells for 15 min. Finally, the invaded A549 cells were quantified under the BX53M upright metallic microscope (Olympus, Tokyo, Japan).

Statistical analysis

Bioinformatic analysis were performed using R software 4.2.1. Two-group differences were calculated by Wilcoxon rank test. The distribution of clinicopathologic features was analyzed by Chi-Squared test. Survival differences among different risk groups or subtypes were compared by Kaplan-Meier curves. Data were shown as mean ± standard deviation, and GraphPad Prism 7.0 was employed to conduct statistical analysis. Student’s t-test was employed to test two-group comparisons. And p < 0.05 was considered statistically significant.

Using NMRGs to classify two molecular subtypes of LUAD

Firstly, six NMRGs (ENPP1, NNT, NT5C1A, NT5C3A, NT5E, and PNP) were significantly correlated with the prognostic outcomes of LUAD (Fig. 1A). Using these six NMRGs, consensus clustering analysis was performed on the 500 samples in TCGA-LUAD dataset. According to the CDF and CDF Delta area curves, the clustering results were stable when there were two clusters (Fig. 1B). Thus, LUAD patients were divided into two subtypes when consensus matrix k = 2 (Fig. 1C). In TCGA-LUAD cohort, the OS probability in C2 subtype was more favorable than that in C1 subtype (Fig. 1D), indicating that C2 subtype patients exhibited a better prognosis than C1 subtype. In addition, 226 samples in GSE31210 dataset were similarly clustered into two stable subtypes, which also showed clear prognostic differences (Fig. 1E).

Clinical features and somatic gene mutation landscape between the two subtypes

Comparison on the distribution difference of clinicopathologic features between C1 and C2 in TCGA-LUAD cohort showed that compared to C2 subtype, more patients in C1 subtype were in the pathologic_M1, pathologic_T3, T4, pathologic_N1, N2, and pathologic_stage II, III (Figs. 2A–2D). In addition, 369 somatic mutation genes with significantly high mutation frequency between C1 and C2 subtypes were obtained and the mutation distribution of the top 20 genes was visually displayed. The top three somatic genes with the highest mutation frequencies were TTN (45.25%), FLG (25.25%), and ZNF536 (19.8%) (Fig. 2E).

DEGs were identified and GSEA was performed between two subtypes

Firstly, DEG analysis between C1 and C2 in TCGA-LUAD cohort screened 178 upregulated DEGs and 61 downregulated DEGs (Fig. 3A). Furthermore, GSEA showed that C1 subtype was specifically enriched in multiple pathways, including the P53_PATHWAY, GLYCOLYSIS, TGF_BETA_SIGNALING, KRAS_SIGNALING_UP, HYPOXIA, which are associated with the progression and deterioration of tumors (Fig. 3B). However, C2 subtype was enriched in the pathways such as MTORC1_SIGNALING, FATTY_ACID_METABOLISM, MYC_TARGETS_V2, and REACTIVE_OXYGEN_SPECIES_PATHWAY (Fig. 3B). These results showed that C1 subtype was more associated with glycolytic metabolism-related pathways, while C2 subtype was more related to fatty acid metabolism pathway.

RiskScore model with a high prediction performance was constructed

Univariate Cox regression analysis screened 88 genes significantly related to LUAD prognosis between C1 and C2 (p < 0.05). To refine the risk model, LASSO Cox regression analysis with 10-fold cross validation was performed to reduce gene number (Figs. 4A and 4B). Subsequently, four independent prognostic NMRGs, including 1 protective gene (CPA3) and three risk genes (DKK1, GJB3, KRT6A), were screened by multivariate Cox regression analysis (Fig. 4C) and used to establish a RiskScore model:

$$\mathrm{R}\mathrm{i}\mathrm{s}\mathrm{k}\mathrm{S}\mathrm{c}\mathrm{o}\mathrm{r}\mathrm{e}=\left(0.095\ast GJB3\right)+\left(-0.141\ast CPA3\right)+\left(0.155\ast DKK1\right)+\left(0.046\ast KRT6A\right).$$

A total of 500 samples in TCGA-LUAD were divided by the median RiskScore value into low-risk and high-risk groups (Fig. 4D). Furthermore, the AUC displayed that the RiskScore model was credible with 1-, 3-, and 5- year AUC of 0.74, 0.68, and 0.62, respectively (Fig. 4E). In TCGA-LUAD cohort, Kaplan-Meier curves demonstrated that low-risk group had better survivals in terms of progression-free interval (PFI), OS, disease-specific survival (DSS), and disease-free interval (DFI) (Figs. 4F–4I), which suggested that patients with a high RiskScore had a worse prognosis.

The RiskScore model was validated in GSE31210 dataset

To confirm the stability and reliability of RiskScore model, the GSE31210 dataset containing 226 tumor samples served as an independent validation set. Similarly, the results in validation cohort were in line with those in TCGA-LUAD training set. In detail, 226 tumor samples in GSE31210 dataset were divided into low- and high-risk groups (Fig. 5A), with low-risk patients having a higher OS probability (Fig. 5C). This suggested that the prognostic outcomes of low-risk patients were more favorable than the high-risk patients. The ROC curve also validated the reliability of the RiskScore model, with 1-, 3-, and 5-year AUC of 0.74, 0.69, and 0.59, respectively (Fig. 5B). Moreover, to further investigate the roles of the four model genes (GJB3, CPA3, DKK1, KRT6A) in LUAD prognosis, 226 patients in GSE31210 dataset were divided by the median expression of these genes into high- and low- expression groups. Except for DKK1, the expressions of GJB3, CPA3, and KRT6A were closely linked to the OS of LUAD patients. The low expression groups of GJB3 and KRT6A exhibited better OS than high expression groups, while the high expression group of CPA3 had higher OS rate than low expression group (Fig. 5D), indicating that the expressions of these model genes could affect LUAD prognosis.

Immune cell infiltration analysis and GSEA between low- and high-risk groups

ESTIMATE, TIMER, ssGSEA, and MCP-counter methods were utilized to compare the degree of immune infiltration between the two risk groups in TCGA-LUAD. Firstly, ESTIMATE analysis revealed that the high-risk group had significantly lower ImmuneScore, StromalScore, and ESTIMATEScore than low-risk group (Fig. 6A), indicating lower immune cell infiltration in the TME of LUAD patients with a high risk. TIMER analysis demonstrated that in comparison to low-risk group, the score of Macrophage was lower in high-risk group (Fig. 6B), implying that the function of Macrophage was inhibited in high-risk patients. The infiltration of 28 types of immune cells was further explored through ssGSEA, and it was observed that compared to low-risk group, the infiltration of most immune cells (such as plasmacytoid dendritic cell, eosinophil, immature dendritic cell, activated CD4 T cell, mast cell, immature B cell, monocyte, natural killer cell, activated B cell) was remarkably lower in high-risk group (Fig. 6C), suggesting an immunosuppressive environment in the high-risk group. Furthermore, MCP-counter analysis demonstrated that the infiltration of neutrophil, T cell, endothelial cell, myeloid dendritic cell and B lineage were significantly reduced in high-risk group (Fig. 6D). The GSEA results showed that the low-risk group was primarily enriched in allograft rejection and bile acid metabolism, whereas the high-risk group was more enriched in numerous HALLMARK pathways (Fig. 7).

High-risk group exhibited poor immunotherapy response and drug sensitivity

With the purpose to further investigate the potential value of the NMRG signature in the precision treatment of LUAD, the efficacy of immunotherapy and 10 common chemotherapy drugs for different risk groups of TCGA-LUAD patients was analyzed. It was found that the TIDE score of high-risk patients was notably higher (Fig. 8A) and the ratio of responders to immunotherapy in this group was lower (Fig. 8B), showing that high-risk LUAD patients had greater immune escape capabilities and less active immunotherapy responses. Additionally, the IC₅₀ values of Erlotinib, Paclitaxel, Cisplatin, Saracatinib, and CGP_082996 in high-risk group were lower than that in low-risk group (Fig. 8C), which indicated that these drugs were more effective in treating high-risk LUAD patients. However, high-risk LUAD patients exhibited low sensitivities to Rapamycin, PHA_665752, Sorafenib, Imatinib, and Crizotinib.

A nomogram with an excellent prediction performance was established

The univariate and multivariate Cox regression analysis verified that N.stage, RiskScore, and T.stage were independent factors for LUAD prognosis (Figs. 9A and 9B). Then, a nomogram was established, and the RiskScore exhibited the strongest influence on the prediction of survival probability (Fig. 9C). In addition, calibration curve showed that 1-year, 3-year, and 5-year calibration points were close to the ideal curves (Fig. 9D). According to DCA, the benefits of the nomogram were significantly higher than the baseline model (Fig. 9E), emphasizing the potential value of the nomogram in clinical decision-making. Moreover, the C-index of the nomogram was also higher than that of clinical features and RiskScore (Fig. 9F), which further confirmed its accurate prediction ability. The ROC curve of the nomogram was plotted, with an AUC value of 0.924 (Fig. 9G), which reflected a strong diagnostic performance of the nomogram.

Silencing GJB3 inhibited the migration and invasion of LUAD cells

The data from qRT-PCR revealed that the mRNA expressions of GJB3, DKK1, CPA3, and KRT6A were all notably upregulated in A549 cells than in BEAS-2B cells (Figs.10A–10D). GJB3 is a member of a group of proteins known as connexins, which fulfill crucial functions in the formation of channels and connections and in turn support the communication between cells (Zeng et al., 2024). However, the role of GJB3 in cancers has been less studied, particularly LUAD. Here, the wound healing and Transwell assays demonstrated that the silencing of GJB3 reduced the numbers of migrated and invaded A549 cells (Figs.10E and 10F).