Characterization of molecular subtypes based on chromatin regulators and identification of the role of NPAS2 in lung adenocarcinoma

Abstract

Background

Chromatin regulators (CRs) are critical epigenetic modifiers and have been reported to play critical roles during the progression of various tumors, but their role in lung adenocarcinoma (LUAD) has not been comprehensively studied.

Methods

Differential expression and univariate Cox regression analyses were conducted to identify the prognostic CRs. Consensus clustering was applied to classify the subtypes of LUAD based on prognostic CRs. LASSO-multivariate Cox regression method was used for construction of a prognostic signature and development of chromatin regulator-related gene index (CRGI). The capacity of CRGI to distinguish survival was evaluated via Kaplan–Meier method in multiple datasets. Relationship between CRGI and tumor microenvironment (TME) was evaluated. Additionally, clinical variables and CRGI were incorporated to create a nomogram. The role of the prognostic gene NPAS2 in LUAD was elucidated via clinical samples validation and a series of in vitro and in vivo experiments.

Results

Two subtypes of LUAD were classified based on 46 prognostic CRs via consensus clustering which had significantly different survival and TME. A prognostic signature consisting of six CRs (MOCS, PBK, CBX3, A1CF, NPAS2, and CTCFL) was developed and proved to be an effective survival predictor in multiple independent datasets. The prognostic signature was also demonstrated to be an indicator of TME and sensitivity to immunotherapy and chemotherapy. The nomogram was suggested to be a simple tool that can predict survival accurately. Clinical samples show that NPAS2 is highly expressed in LUAD tissues, and in vitro and in vivo experiments demonstrated that inhibition of NPAS2 impeded malignant progression of LUAD cells.

Conclusions

Our study comprehensively unveiled the functions of CRs in LUAD, developed a classifier to predict survival and response to treatments, and suggested that NPAS2 promoted LUAD progression for the first time.

Introduction

Lung cancer is one of the most frequently diagnosed cancers which is made up of diverse histological subtypes, with lung adenocarcinoma representing approximately 54.7% of all lung cancer cases [1]. Despite the encouraging progresses achieved in screening and treatment, the prognosis of patient with LUAD is not optimistic as only 26.4% of LUAD individuals can survive 5 years or longer [1]. Application of predictive biomarkers can aid in stratifying patients with varying degrees of mortality risk. Poor survival of LUAD makes it necessary to discover more robust biomarkers to achieve individual treatment.

CRs consisting of DNA methylators, histone modifiers, and chromatin remodelers can control epigenetic alteration which is one of the cancer hallmarks [2]. CRs are crucial modulators in a variety of biological processes such as energy metabolism and activation of T cells [3, 4]. Dysregulation of CRs affects chromatin structure and expression of cancer-related genes, leading to changes in biological processes and malignant behaviors [2, 5]. Perturbation of CRs is common in various cancers, and mutation or aberrant expression of CRs was associated with outcome of malignancies [6]. For example, DNMT1 which belongs to the family of DNA methyltransferases not only regulated stemness and tumorigenicity of liver cancer and breast cancer, but also predicted poor prognosis of cervical cancer and breast cancer [7,8,9]. Another study has revealed that CRs served as a predictor of survival and therapeutic response in bladder cancer [10]. The above studies suggest that CRs possess great potential to be reliable predictors of prognosis.

In this study, we comprehensively analyzed expression of CRs in LUAD and identified the prognostic CRs to characterize molecular subtypes of LUAD. Additionally, we developed a prognostic predictor based on prognostic CRs to predict survival and TME for LUAD individuals. We also investigated the role of NPAS2 in progression of LUAD via a series of in vitro and in vivo experiments.

Materials and methods

Identification of prognostic CRs

Transcriptome, methylation data, and genome files of the LUAD samples were collected from The Cancer Genome Atlas (TCGA) database (https://gdc-portal.nci.nih.gov/). Chromatin regulators (CRs) were collected from previous studies [2] (Additional file 1: Table S1). After removal of genes with low expression and data normalization, the differentially expressed CRs were identified by “edgeR” package based on the screening criteria of adjusted p value (FDR) < 0.05 and |log2 fold change (FC)|> 1.0. The differentially expressed CRs which had statistically significant effect on survival in in univariate Cox regression analysis were considered prognostic CRs. Mutation landscape and copy number variations (CNV) of the CRs were assessed based on genome data of the TCGA-LUAD cohort.

Subtype characterization of LUAD

Expression data of GSE37745, GSE31210, and GSE50081 were collected from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). After those without complete clinicopathological information or followed up less than a month were removed, LUAD individuals from TCGA, GSE37745, and GSE31210 were merged into a large LUAD cohort which consisted of 800 samples. “SVA” package which used the ComBat function to remove both known batch effects and other potential latent sources of variation was applied in our study to remove the batch effect between the LUAD datasets [11]. We summarized the clinicopathological variables of all the LUAD samples in Additional file 2: Table S2. The consensus clustering algorithm, an unsupervised clustering method, is able to detect potential groups of a dataset based on intrinsic characteristics, making it a popular data mining technique in cancer research [12]. In this study, consensus clustering was executed by “ConsensusClusterPlus” package to characterize the subtypes of LUAD based on expression of prognostic CRs [13]. Principal component analysis (PCA) was carried out to confirm the capability of the prognostic CRs in distinguishing the clusters identified in consensus clustering. Kaplan–Meier survival curves were created to test the prognostic value of the CR subtypes.

The relationship of CR-related subtypes with methylation and TME

The chromatin regulators control epigenetic alteration including DNA methylation. The CpG sites with significantly different methylation levels between the CR-related subtypes were screened according to β value using “limma” package, and the cutoff criteria were FDR < 0.05. TME exerts tremendous influence in tumorigenesis and is associated with prognosis of patients. The TME is mainly composed of stromal cells and immune cells. The abundance of immune cells, stromal cells, and tumor cells was quantified via Estimation of STromal and Immune cells in MAlignant Tumors using Expression data (ESTIMATE) method by “estimate” package. There are various subtypes of immune cells which fulfill different functions. Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) method was utilized to calculate the abundance of the major subtypes of immune cells. Gene set variation analysis (GSVA) was executed to investigate which pathways differed between the two subtypes with FDR < 0.05 as screening criteria.

Construction of a CR-associated signature

Further selection of prognostic CRs was carried out by Least Absolute Shrinkage and Selection Operator (LASSO) Cox regression analysis. Then, multivariate Cox regression analysis was used to construct a signature based on the CRs selected in LASSO method. CRGI was calculated via mRNA levels of CRs and corresponding coefficients identified in the multivariate Cox regression analysis, and the median value of CRGI was set as the division criterion to classify patients into high- and low-CRGI groups. The relationship between the subtypes and CRGI groups was presented by Sankey diagram.

Validation of the prognostic signature

PCA was conducted in the TCGA cohort to confirm the capability of the CRs in the signature to discriminate between those with different CRGIs. The survival difference between the high- and low-CRGI groups was displayed by Kaplan–Meier survival curves. We also calculated the survival percentage of the high- and low-CRGI groups and the distribution of CRGI between survivors and nonsurvivors. The exact effect of CRGI on survival was assessed by univariate and multivariate Cox regression analyses. The predictive capability of CRGI was inferred by receiver operating characteristic (ROC) curves. To further validate the prognostic value of CRGI, PCA and Kaplan–Meier survival analysis were carried out in GSE37745 and GSE31210.

TME patterns of the high- and low-CRGI groups

As has been described above, ESTIMATE and CIBERSORT algorithms were run to investigate the patterns of TME in the high- and low-CRGI groups. Additionally, the correlation between CRGI and immune cells was also evaluated. Expression of immune checkpoint genes is associated with function of immune cells and therapeutic response of immunotherapy. Thus, the relationship between the expression of immune checkpoint genes and CRGI was assessed.

Construction and assessment of a nomogram

Considering that clinicopathological variables are contributing factors of survival, we developed a nomogram consisting of clinicopathological parameters and CRGI to improve the predictive precision of CRGI. Clinicians can use the points obtained from the nomogram to predict survival probability of patients. Calibration curves and ROC curves were generated, and C-index and area under the curve (AUC) were calculated to evaluate predictive capability of the nomogram, CRGI, and other clinical variables.

Relationship between CRGI and sensitivity to treatment

The sensitivity of the high- and low-CRGI groups to immunotherapy and chemotherapy which are effective treatments for LUAD individuals was also predicted. Response of the high- and low-CRGI groups to immunotherapy and chemotherapy was compared to figure out whether CRGI could be a tool for stratification of therapeutic benefits. Immunophenoscore (IPS) was developed based on expression of immune-related genes, which can predict response to CTLA-4 and PD-1 antibodies [14]. IPS was positively associated with efficacy to immunotherapy [14]. IPS of 20 solid cancers in the TCGA database is provided in The Cancer Immunome Atlas (TCIA) (https://tcia.at/home). We collected IPS of the TCGA-LUAD cohort and compared IPS between the high- and low-CRGI groups. Chemotherapeutic agents such as cisplatin and paclitaxel can achieve satisfactory efficacy in patients with LUAD. Half-maximal inhibitory concentration (IC50) of chemotherapeutic agents was acquired using “pRRophetic” package to evaluate therapeutic benefits of chemotherapy in the high- and low-CRGI groups.

Clinical specimens and immunohistochemistry (IHC)

We also collected lung adenocarcinoma tissues and paired normal lung tissues of 11 patients from Wuhan Tongji Hospital for IHC analysis. According to the manufacturer’s instructions, immunohistochemical staining of NPAS2 was conducted. The specific primary antibodies anti-NPAS2 (Invitrogen, PA5-98824) was used for IHC. The IHC scores of each patient were calculated to evaluate the NPAS2 protein expression levels based on the staining intensity and quantity.

Cell culture and siRNA transfection

The human lung adenocarcinoma cell lines including PC-9, A549, HCC827, H1975 and the normal lung epithelial cell line BEAS-2b were stored in the oncology laboratory of Tongji Hospital, Wuhan, China. All these cells were grown in RPMI-1640 medium containing 10% FBS (Gibco, USA) and maintained in the 37 °C incubator with 5% CO₂. The siRNA oligonucleotides (negative control siRNA and NPAS2 siRNA) were synthesized by GENERAL BIOL (Anhui, China), and siRNA transfection was carried out with Lipofectamine 3000 (Invitrogen, USA). The siRNA sequences were: si-NC, 5′-UUCUCCGAACGUGUCACGUTT-3′ (Sense) and 5′-ACGUGACACGUUCGGAGAATT-3′ (Antisense); si-NPAS2#1, 5′-CCUCAGCACUAAAGGACAATT-3′ (Sense) and 5′-UUGUCCUUUAGUGCUGAGGTT-3′ (Antisense); si-NPAS2#2, 5′-CCACCAAGCUGAUGGCAGATT-3′ (Sense) and 5′-UCUGCCAUCAGCUUGGUGGTT-3′ (Antisense).

Quantitative real-time PCR (qRT-PCR)

Total cell RNA was isolated using TRIzol reagent (TaKaRa, Japan) and transcribed to cDNA using the reverse transcription reagents (Vazyme, China). Then, the qRT-PCR was conducted with a Real-Time PCR System (Applied Biosystems 7900HT, USA). The relative expression levels of NPAS2 were calculated using the 2^−ΔΔCT method. The primer sequences were listed as follows: NPAS2-forward, 5′-CACAGAGCACCTCCAATCATAG‑3′ and NPAS2-reverse, 5′-GTAGCAACACGACTTCCCTT‑3′; GAPDH-forward, 5′-GACAGTCAGCCGCATCTTCT‑3′ and GAPDH-reverse, 5′-GCGCCCAATACGACCAAATC‑3′.

Western blot

Western blot assay was conducted as previously described [15]. Briefly, the extracted protein samples were separated by SDS-PAGE and transferred to PVDF membranes. Next, the membranes were blocked using 5% nonfat milk and incubated with specific primary antibodies (GAPDH, 1:5000, 10494-1-AP, Proteintech and NPAS2, 1:1000, PA5-98824, Invitrogen). After incubation overnight at 4 °C, the membranes were washed and incubated with secondary antibodies (1:5000, SA00001-2, Proteintech). Finally, the protein bands were detected with the ECL reagent by the G: BOX Chemi X system (Syngene, UK).

CCK8 assay

Two thousand to three thousand cells/well were plated into 96-well plates. After adherent, the cells were incubated with 100 μl medium containing 10% CCK8 (MedChem Express, USA) each well for 1–2 h at 37 °C. The absorbance (450 nm wavelength) was detected with a microplate reader (BioTek, Winooski, VT, USA). After culturing cells for another 24, 48, and 72 h, these operations were repeated.

EdU assay

Five thousand cells/well were plated into 96-well plates. After overnight attachment, different reagents were added successively according to the instructions of EdU kit (Meilunbio, China). The EdU-positive cells were detected using a fluorescence microscope (Leica).

Colony formation assay

Cells were seeded in 6-well plates (500 cells/well) and cultured for 10 days. Then, colonies were fixed (4% formaldehyde), stained (0.1% crystal violet), photographed, and counted.

Wound healing assay

Cells were plated into 12-well plates and cultured to confluence. Then, cells were scratched with 200 μl pipette tips and cultured in medium containing 1% FBS. At 0 h and 48 h, wound closure was observed and micrographs were captured with an optical microscope (Nikon).

Transwell assays

Transwell assays were performed using the 24-well transwell chambers (8-μm pores, NEST Biotechnology). For invasion assay, the upper chambers were precoated with 50 μl diluted Matrigel (1:8) and maintained in 37 °C incubator overnight. 50,000 cells in 200 μl serum-free medium were placed into the upper chambers, and 600 μl medium containing 20% FBS was placed into the lower chambers. After culture for another 24 h or 48 h, the nonmigrated cells were removed using cotton swabs and the migratory cells across the aperture were fixed (4% paraformaldehyde), stained (0.1% crystal violet), photographed, and counted.

Animal experiments

According to the target sequences of si-NC and si-NPAS2#1, we constructed the shRNA lentivirus vectors and transfected into A549 cells with polybrene. The male BALB/c nude mice (five weeks old) were used to construct xenograft model by subcutaneously injecting NC and sh-NPAS2 A549 cells (5 × 10⁶) into the right axilla of each mouse. The tumor volumes were measured with the formula: tumor volume (mm³) = 1/2 (length × width²). At the experimental endpoint, the mice were killed, and the xenograft tumors were separated and fixed for further IHC analyses. The primary antibodies were listed as follows: NPAS2 (Invitrogen, PA5-98824), Ki-67 (Abcam, ab15580).

Statistical analysis

The analyses performed in our study were executed via R software (version 4.1.0) and GraphPad 8.0. Assessment of differences between two groups was made by the Wilcoxon test or Student’s t test. In Kaplan–Meier survival analysis, log-rank test was utilized to compare the survival time of different groups. Correlation between numeric variables was identified by Spearman correlation analysis. p value < 0.05 was considered statistically significant.

Results

Identification of prognostic CRs

The workflow of this study is presented in Fig. 1. A total of 535 LUAD and 59 para-tumoral lung samples from TCGA were included in the differential expression analysis which identified 117 dysregulated CRs in LUAD, including 21 downregulated CRs and 96 upregulated CRs (Fig. 2A). Among all the dysregulated CRs, 46 CRs were identified as contributing factors of overall survival (Fig. 2B). Somatic mutations of the prognostic CRs occurred in 37.79% of the LUAD samples with PRDM9 had the highest mutation frequencies (Fig. 2C). Additionally, CNV of the prognostic CRs can be observed in LUAD samples (Fig. 2D).

Fig. 1

Workflow of this study

Fig. 2

Expression and genome landscape of prognostic CRs. A Expression of differentially expressed CRs between LUAD and peritumoral tissues. B Forest plot showing effect of 46 prognostic CRs on survival. C Somatic mutation landscape of 46 prognostic CRs. D CNV frequencies of 46 prognostic CRs

Characterization of LUAD subtypes based on prognostic CRs

A total of 800 LUAD samples were subjected to consensus clustering based on expression of 46 CRs, and 2 discrete subtypes were identified in consensus clustering (Fig. 3A). There are 318 cases in C1 and 482 cases in C2. PCA confirmed that prognostic CRs can separate C1 from C2 (Fig. 3B). Cluster 1 showed significant survival advantage over cluster 2 and had lower expression of most of the prognostic CRs than cluster 2 (Fig. 3C,  D). A total of 3827 differentially methylated sites were identified between the high- and low-CRGI groups, suggesting that patients with different CRGI had different DNA methylation levels. The most significant sites in differential methylation analysis were displayed in heatmap (Additional file 3: Fig. S1). Cluster 1 had higher abundance of stromal cells and immune cells and lower tumor purity than cluster 2 (Fig. 4A–D). Immune cell including plasma cells, monocytes, dendritic cells, and mast cells are more abundant in cluster 1 (Fig. 4E). Pathways such as bladder cancer, pancreatic cancer, and p53 signaling pathway are activated in cluster 2, while asthma, cell adhesion molecules cams, complement, and coagulation cascades are activated in cluster 1 (Fig. 4F).

Fig. 3

Identification of two subtypes of LUAD. A Two subtypes identified by consensus clustering. B PCA plot showing distribution of the two subtypes based on 46 prognostic CRs. C Survival difference between the two subtypes. D Differential expression of 46 prognostic CRs between two subtypes

Fig. 4

TME patterns and related-pathways of the two subtypes. A–D Differences in ESTIMATE score, tumor purity, immune score, and stromal score between the two subtypes. E Differences in immune cell subtypes between the two subtypes. F Differences in pathways between the two subtypes

Construction of a signature based on CRs

Given that the prognostic CRs can distinguish LUAD individuals with favorable survival from those with poor survival, we constructed a signature based on prognostic CRs to predict survival for patients with LUAD. A signature made up of 6 prognostic CRs (MOCS, PBK, CBX3, A1CF, NPAS2, and CTCFL) was developed by stepwise LASSO and multivariate Cox regression analyses (Additional file 4: Fig. S2A–C). Among the CRs in the signature, MOCS1 was tumor suppressor gene affecting survival positively, while the other genes were oncogenes whose high expression was associated with poor survival (Additional file 4: Fig. S2C). CRGI was computed using expression level of CRs and coefficients in the signature, and LUAD samples were classified into high-CRGI group and low-CRGI group (Additional file 5: Table S3). Most of the patients in the cluster 1 were classified into the low-CRGI group, and the low-CRGI group accounted for a high percentage of the living patients (Additional file 4: Fig. S2D), corresponding to the results that cluster 1 had lower CRGI than cluster 2 (Additional file 4: Fig. S2E).

Validation of the prognostic signature

Multiple analyses were conducted in the TCGA-LUAD dataset to show the prognostic value of CRGI. The high-CRGI group can be separated from the low-CRGI group by the expression of the CRs in the signature (Fig. 5A). Increased CRGI can be observed in nonsurvivors compared with survivors (p < 0.001) (Fig. 5B). Furthermore, CRGI is a risk factor of survival whose effect on survival is independent of the clinicopathological variables [hazard ratio (HR) = 1.641, 95% confidence interval (CI): 1.467–1.836, p < 0.001] (Fig. 5C, D). The AUCs at 1, 3, and 5 years were 0.740, 0.692, and 0.690, respectively (Fig. 5E). The prognostic value of CRGI was further examined in GSE37745 and GSE31210. Low CRGI was associated with survival advantage in both GSE37745 (p < 0.001) and GSE31210 (p = 0.0076) (Additional file 6: Fig. S3A, B). Different expression patterns of the CRs in the signature can be observed between the high- and low-CRGI groups in both GSE37745 and GSE31210 (Additional file 6: Fig. S3C, D). Cell cycle is fundamental biological process, and cell cycle-related genes are significant contributing factor of prognosis of patients with cancer [16, 17]. Some CRs in our model have been reported to be involved in cell cycle including PBK [18], CBX3 [19,20,21], and NPAS2 [22]. To clarify whether the prognostic value of CRGI depended on cell cycle-related genes, we removed these genes and re-validate the prognostic value of CRGI. Analyses of survival and TME suggested that CRGI was still a predictor of outcome and immune infiltration (Additional files 7 and 8: Figs. S4, S5).

Fig. 5

The prognostic value of CRGI in the TCGA dataset. A PCA plot showing distribution of the high- and low-CRGI groups. B Survival difference between the high- and low-CRGI groups. C Univariate Cox regression of CRGI and clinicopathological variables. D Multivariate Cox regression of CRGI and clinicopathological variables. E ROC curves of CRGI in survival prediction

Differences in TME between the high- and low-CRGI groups

Compared those with the high CRGI, LUAD samples with low CRGI had higher abundance of immune cells and stromal cells (Fig. 6A). Additionally, CRGI was negatively associated with immune cells such as the mast cells, eosinophils, activated NK cells, monocytes, and CD4⁺ memory cells and positively associated with neutrophils and macrophages (Fig. 6B). The CRs in the prognostic signature were also correlated with immune infiltration (Fig. 6B). Differential expression of immune checkpoint genes such as PD-L1, LAMA3, IFNG, GZMB, and CD28 can be observed between the high- and low-CRGI groups (Fig. 6C).

Fig. 6

Correlation between CRGI and TME. A Differences in ESTIMATE score, tumor purity, immune score, and stromal score between high- and low-CRGI groups. B Correlation between CRGI and infiltration of immune cells. C Differences in expression of immune checkpoint between high- and low-CRGI groups

Construction and validation of a nomogram

A nomogram which consisted of CRGI and clinicopathological variables was constructed and validated in TCGA, GSE37745, GSE31210, and GSE50081 (Fig. 7A). The high consistency between predicted survival and actual survival was presented via the calibration curves (Fig. 7B–E). The results of C-index illustrated that the predictive accuracy of CRGI and clinicopathological parameters was enhanced by incorporating them into a nomogram (Fig. 7B–E). High accuracy of the nomogram was also observed from the ROC curves, indicating that the nomogram is a reliable tool for survival prediction (Fig. 7B–E).

Fig. 7

Development and validation of a nomogram. A Nomogram based on CRGI and clinicopathological parameters. B Calibration curves, C-index, and ROC curves of the nomogram in survival prediction in the TCGA. C Calibration curves, C-index, and ROC curves of the nomogram in survival prediction in GSE37745. D Calibration curves, C-index, and ROC curves of the nomogram in survival prediction in GSE31210. E Calibration curves, C-index, and ROC curves of the nomogram in survival prediction in GSE50081

CRGI was associated with response to antitumor treatment

The low-CRGI group showed higher IPS which represents immunogenicity in “CTLA4⁻ PD1⁻” (p < 0.001), “CTLA4⁻ PD1⁺” (p = 0.025), “CTLA4⁺ PD1⁻” (p < 0.001), and “CTLA4⁺ PD1⁺” (p = 0.025) groups (Additional file 9: Fig. S6A), indicating that better efficacy of anti-CTLA-4 treatment and anti-PD-1 treatment may be achieved in the low-CRGI group. IC50 values of docetaxel, paclitaxel, and cisplatin are lower in the high-CRGI group compared with the low-CRGI group (p < 0.001) (Additional file 9: Fig. S6B), which suggests that the high-CRGI group are more sensitive to these chemotherapeutic agents. CRGI is not only a survival predictor but also a stratification tool for making optimal treatment.

NPAS2 is highly expressed and predicts poor prognosis in LUAD

NPAS2 is the most significant contributing factor of survival in multivariate Cox regression analysis among all the signature genes. Therefore, we further investigated the mRNA expression levels and prognostic values of NPAS2 and found that NPAS2 is upregulated in LUAD tissues compared to normal lung tissues and predicts shorter OS (Fig. 8A–E). Furthermore, we also validated the high expression of NPAS2 in LUAD at protein levels by using own clinical samples based on IHC staining (Fig. 8F, G).

Fig. 8

NPAS2 is overexpressed and associated with poor prognosis in LUAD. A–B The mRNA expression of NPAS2 in TCGA-LUAD and GSE31210 datasets. C–E Kaplan–Meier curves of LUAD patients based on the optimal cutoff in TCGA-LUAD, GSE37745, and GSE31210 datasets. F–G Representative images of IHC staining of NPAS2 protein and the comparison of IHC score in normal lung and LUAD tissues

Knockdown of NPAS2 dampened malignant behaviors of LUAD cells

The oncogenic properties of NPAS2 have been demonstrated in various cancers like hepatocellular carcinoma and thyroid carcinoma, but the role of NPAS2 in LUAD remains unexplored. In this study, we first investigated the mRNA expression of NPAS2 in BEAS-2b, a normal lung epithelial cell line, and four LUAD cell lines (A549, PC-9, HCC827, and H1975) via qRT-PCR, indicating that NPAS2 expression was significantly elevated in all LUAD cells compared with BEAS-2b (Fig. 9A). Due to the highest expression of NPAS2 in A549 and PC-9 cells, they were selected for further experiments with siRNA knocking down. The mRNA and protein level of NPAS2 was remarkably suppressed by si-NPAS2#1 and si-NPAS2#2 in both A549 and PC-9 cells (Fig. 9B–D). The effect of NPAS2 on LUAD cells proliferation was evaluated by CCK8, EdU, and colony formation assays. The CCK8 assays indicated that NPAS2 knockdown impaired cell viability of A549 and PC-9 cells (Fig. 9E, F), and the results of EdU assay showed that inhibiting NPAS2 dramatically reduced the proportion of EdU-positive cells (Fig. 9G, H). Additionally, colony formation ability of LUAD cells was remarkably attenuated after NPAS2 loss (Fig. 9I, J). Collectively, these results suggested that NPAS2 knockdown dramatically suppressed LUAD cells proliferation in vitro.

Fig. 9

Knockdown of NPAS2 impeded proliferative activity of LUAD cells. A Expression of NPAS2 in BEAS-2b, A549, PC-9, HCC827, and H1975 cells. B–D Efficacy of siRNAs in suppressing expression of NPAS2 in A549 and PC-9 cells. E–F OD values of A549 and PC-9 cells in CCK8 assay after knockdown of NPAS2. G–H Proliferative activity of A549 and PC-9 cells in EdU assay after knockdown of NPAS2. I–J Colony formation assay for A549 and PC-9 cells after knockdown of NPAS2. *p < 0.05, ** p < 0.01, *** p < 0.001

Next, to assess the effect of inhibiting NPAS2 on migration and invasion, wound healing and transwell assays were conducted in vitro. In wound healing assays, A549 and PC-9 cells with NPAS2 suppressed migrated more slowly compared to the control group (Fig. 10A, B). Additionally, the A549 and PC-9 cells with NPAS2 knockdown exhibited weaker migration and invasive abilities in transwell assays (Fig. 10C–F). All the above results suggest that NPAS2 is a cancer driver gene in LUAD, and that inhibition of NPAS2 suppresses malignant progression of LUAD cells.

Fig. 10

Knockdown of NPAS2 impeded migration and invasion capacities of LUAD cells. A–B Migration ability of A549 and PC9 cells in wound healing assay after knockdown of NPAS2. C–F Migration and invasion abilities of A549 and PC-9 cells in transwell assays after knockdown of NPAS2. *p < 0.05, ** p < 0.01, *** p < 0.001

NPAS2 knockdown suppresses tumor growth of LUAD in vivo

To further determine the oncogenic role of NPAS2 in LUAD in vivo, we construct xenograft model by subcutaneously injecting NC and sh-NPAS2 A549 cells in BALB/c nude mice. As we expected, we found that NPAS2 knockdown significantly inhibited tumor growth, and the tumor volume and weight of the NPAS2 knockdown group were notably smaller than the NC group (Fig. 11A–C). Further IHC staining confirmed that NPAS2 was down-expressing in the tumor tissues of NPAS2 knockdown group. Besides, the Ki-67 expression level was also significantly less than NC group (Fig. 11D, E). Collectively, these results demonstrated that NPAS2 could promote tumor growth of LUAD in vivo and function as an oncogene of LUAD.

Fig. 11

Knockdown of NPAS2 inhibited the growth of A549 cells in vivo. A Gross appearance of tumor xenografts in NC and sh-NPAS2 groups. n = 5. B Tumor growth curve of two groups. C The tumor weights on day 24 of each group. D–E HE staining and IHC staining of NPAS2 and Ki-67 in the tumor tissues. *p < 0.05, **p < 0.01, ***p < 0.001

Discussion

High malignancy and poor outcome of LUAD make it necessary to explore more robust prognostic indicators. Chromatin regulators are key molecules that affect numerous biological processes and malignant progression. In this study, we characterized two subtypes of LUAD based on expression of CRs which had significantly different prognosis and TME. Compared with the cluster 2, the cluster 1 had better survival, lower tumor purity, and higher infiltration degree of immune cells. Considering the prognostic value of CR subtypes, we construct a CR-based indicator to predict survival and TME for LUAD individuals.

Six key CRs were selected to develop a signature which was used to calculate CRGI. The survival probability of LUAD individuals decreased with the increase in CRGI in multiple LUAD datasets, suggesting that CRGI is a remarkable risk factor for patients with LUAD. The negative effect of CRGI on survival is independent of clinicopathological parameters, and the effect of CRGI on survival is more significant than the clinicopathological parameters. Therefore, CRGI had the potential to be alternative or complement to clinical parameters.

All the genes in the signature were oncogenes except MOCS1 whose increased expression was associated with favorable prognosis. NPAS2 is a circadian clock gene and has been reported to be a crucial regulator of tumorigenesis and immune infiltration. NPAS2 was identified as an oncogene in hepatocellular carcinoma and was related to infiltration of immune cells [23]. In accordance with our findings, increased expression of NPAS2 has been suggested to indicate poor survival in lung cancer [24, 25]. NPAS2 was overexpressed in anaplastic thyroid carcinoma and facilitated malignant progression [26]. The oncogenic properties of NPAS2 were associated with its role in promoting the Warburg effect [27]. In our study, NPAS2 was demonstrated to prompt malignant progression of LUAD via a series of in vitro experiments, which was consistent with the bioinformatic findings and precious studies. MOCS1 was essential for synthesis of molybdenum cofactor, and gene mutation of MOCS1 was associated with molybdenum cofactor deficiency [28,29,30]. Upregulation of MOCS1 was reported in patients with higher arterial stiffness [31]. However, the role of MOCS1 in oncogenesis has not been studied. Our study identified the prognostic role of MOCS1 in lung cancer for the first time. A1CF, also known as APOBEC1 complementation factor, regulates various cellular processes and exerts a carcinogenic role in renal cell carcinoma, endometrial cancer, and glioma through RNA editing function [32,33,34]. Inhibition of A1CF dampens the malignant behaviors of cancer cell [32]. Although CTCFL is not expressed in normal tissues except testis and embryonic stem cells [35, 36], increased expression of CTCFL could be observed in various carcinomas and was correlated with malignant behaviors and drug resistance [37,38,39,40]. CTCFL can induce expression of oncogenes and promote oncogenic properties of carcinomas such as gastric cancer and ovarian cancer [41,42,43]. Upregulation of CTCFL is required for resistance to ALK inhibition in cancer cells via promoting chromatin interactions [37]. CTCFL can also serve as a promising target of immunotherapy to improve prognosis of breast cancer patients [44]. PBK is upregulated in hepatocellular carcinoma which is a risk factor of survival and associated with antitumor immunity in hepatocellular carcinoma [45, 46]. Additionally, PBK promotes metastasis and resistance to oxaliplatin by regulating pathways such as ETV4-uPAR and PTEN in hepatocellular carcinoma [47, 48]. The prognostic value of PBK in prediction of survival and antitumor immunity was also proved in colorectal cancer [49,50,51]. Therefore, the prognostic CRs identified in our study also exert fundamental role in various cancers and have potential to be reliable survival predictor and therapeutic targets.

Clinicopathological variables are frequently used and effective parameters to predict prognosis of cancer in clinical practice. Therefore, CRGI and clinicopathological variables were incorporated to establish a nomogram, and survival probability of patients with LUAD could be estimated accurately according to the points obtained in nomogram. The nomogram had the highest C-index among all the predictors, which demonstrates that combining CRGI and clinicopathological variable can improve predictive ability of single parameter and that the nomogram we established is a reliable and simple tool to predict survival.

CRGI was associated with components of TME, including immune cell and stromal cell. Additionally, expression of immune checkpoint differs between LUAD samples with different CRGIs. The high-CRGI group had higher IPS than the low-CRGI group, suggesting that those with high CRGI are more suitable for immune checkpoint inhibitors. Although the low-CRGI group are more likely to respond to immunotherapy, the high-CRGI group are more sensitive to chemotherapy such as docetaxel, paclitaxel, and cisplatin since the high-CRGI group had higher IC50 values than the low-CRGI group. Thus, CRGI is not only a predictor of survival but also a useful tool for treatment selection for patients with LUAD.

In conclusion, this study characterized two CR-related subtypes of LUAD which had varying prognosis and TME. We also constructed a signature based on CRs which was used to calculate an index termed CRGI to predict survival, TME, and sensitivity to immunotherapy and chemotherapy for LUAD, contributing to making personal treatment plans for LUAD individuals. Additionally, a nomogram with high predictive accuracy was established as a simple tool for survival prediction in LUAD. Our study aids in understanding CRs in LUAD and exploring novel prognostic predictor for patients with LUAD.