Advertisement

Identification and Characterization of Senescence Phenotype in Lung Adenocarcinoma with High Drug Sensitivity

  • Seema Khadirnaikar
    Affiliations
    Department of Biosciences and Bioengineering, Indian Institute of Technology Dharwad, Dharwad, India

    Department of Electrical Engineering, Indian Institute of Technology Dharwad, Dharwad, India
    Search for articles by this author
  • Annesha Chatterjee
    Affiliations
    Department of Biosciences and Bioengineering, Indian Institute of Technology Dharwad, Dharwad, India
    Search for articles by this author
  • Sudhanshu Shukla
    Correspondence
    Address correspondence to Sudhanshu Shukla, Ph.D., Department of Biosciences and Bioengineering, Indian Institute of Technology–Dharwad, Dharwad, Karnataka 580011, India.
    Affiliations
    Department of Biosciences and Bioengineering, Indian Institute of Technology Dharwad, Dharwad, India
    Search for articles by this author
      Lung adenocarcinoma (LUAD) is a major health problem and has poor prognosis. Heterogeneity is a central determinant of the treatment outcome, requiring identification of new subclasses of LUAD. Senescence has emerged as a crucial regulator of metastasis and drug response. Ionizing radiation– and doxorubicin-induced senescence-associated genes in lung fibroblasts were used in K-means clustering to identify high- and low-senescence (HS and LS) classes among The Cancer Genome Atlas- LUAD (TCGA-LUAD) patients. The LS group showed significantly poorer survival (P = 0.01) and greater activation of proliferative signaling pathways, proliferation, wound healing, and genetic aberrations (P < 0.05). The TP53 mutation rate was significantly greater in the HS group (P < 0.0001), explaining the phenotype. Also, genome-wide hypomethylation was significantly greater in the LS group than in the HS group. Interestingly, pathway analysis identified silencing of Wnt signaling in the HS group. The machine learning–based recursive feature elimination technique was used to identify a 20-gene senescence signature in TCGA-LUAD samples. The presence of a senescence phenotype with poor survival was validated in an independent patient cohort and a cell-line cohort using unsupervised clustering of samples based on a 20-gene signature. On further analysis, HS cells were more resistant to drugs, particularly histone deacetylase inhibitors. Taken together, this study identified a novel subtype of LUAD with reduced Wnt signaling and high drug resistance.
      Lung adenocarcinoma (LUAD) is a common cause of cancer-related deaths worldwide.
      • Bade B.C.
      • Cruz C.S.D.
      Lung cancer 2020 epidemiology, etiology, and prevention.
      ,
      • Duma N.
      • Santana-Davila R.
      • Molina J.R.
      Non–small cell lung cancer: epidemiology, screening, diagnosis, and treatment.
      Even with advancements in therapy protocols, the prognosis of LUAD patients remains dismal.
      • Duma N.
      • Santana-Davila R.
      • Molina J.R.
      Non–small cell lung cancer: epidemiology, screening, diagnosis, and treatment.
      Inherent heterogeneity is a common cause of the variable response to drugs. Thus, it is crucial to identify novel subclasses of LUAD to design better therapy protocols.
      • Bade B.C.
      • Cruz C.S.D.
      Lung cancer 2020 epidemiology, etiology, and prevention.
      ,
      • Duma N.
      • Santana-Davila R.
      • Molina J.R.
      Non–small cell lung cancer: epidemiology, screening, diagnosis, and treatment.
      Intertumoral heterogeneity in tumors is caused by varying genetic aberrations, immune microenvironment, metastatic capacity, and senescence.
      • Kim Y.H.
      • Choi Y.W.
      • Lee J.
      • Soh E.Y.
      • Kim J.-H.
      • Park T.J.
      Senescent tumor cells lead the collective invasion in thyroid cancer.
      Of these, senescence is a crucial regulatory mechanism. Senescence is characterized as the aging and finite proliferative capacity of normal cells.
      • He S.
      • Sharpless N.E.
      Senescence in health and disease.
      However, it is also associated with a typical cellular response to various stimuli in numerous pathologic and physiologic processes.
      • He S.
      • Sharpless N.E.
      Senescence in health and disease.
      Recent studies have recognized senescence for both its anticancer and procancer phenomena.
      • Wang B.
      • Kohli J.
      • Demaria M.
      Senescent cells in cancer therapy: friends or foes?.
      A recent study used transcriptomic analysis to show that the senescence response of cells depends on the stimulus and the cell type.
      • Casella G.
      • Munk R.
      • Kim K.M.
      • Piao Y.
      • De S.
      • Abdelmohsen K.
      • Gorospe M.
      Transcriptome signature of cellular senescence.
      These observations indicate that tumors have heterogeneity in terms of senescence status.
      In this work, transcriptomic data from senescent lung fibroblast cells were used to identify the senescence subtypes of LUAD. Low-senescence (LS) samples were characterized by significant hypomethylation, greater genetic aberration, and poor prognosis. A classification model to identify senescence levels in LUAD patients and cell lines is proposed based on a machine-learning technique.

      Materials and Methods

      Lung Fibroblast Senescence Expression Data and the Identification of Senescence-Associated Genes and Patient Data Sets

      Data on senescence from WI-38, lung fibroblast cells treated with doxorubicin and ionizing radiation (IR) from the GSE130727 data set, were used. To select senescence-associated genes, the expressions of genes in senescent cells versus proliferating control cells were compared, and a fold-cutoff of >5 was applied. To select the senescence-associated genes with a probable role in cancer, the expression of genes was compared between LUAD and normal lung samples. RNA sequencing (RNA-seq) data associated with clinical parameters, from The Cancer Genome Atlas (TCGA)–Genomic Data Commons, Michigan
      • Shukla S.
      • Evans J.R.
      • Malik R.
      • Feng F.Y.
      • Dhanasekaran S.M.
      • Cao X.
      • Chen G.
      • Beer D.G.
      • Jiang H.
      • Chinnaiyan A.M.
      Development of a RNA-seq based prognostic signature in lung adenocarcinoma.
      ,
      • Dhanasekaran S.M.
      • Balbin O.A.
      • Chen G.
      • Nadal E.
      • Kalyana-Sundaram S.
      • Pan J.
      • Veeneman B.
      • Cao X.
      • Malik R.
      • Vats P.
      • Wang R.
      • Huang S.
      • Zhong J.
      • Jing X.
      • Iyer M.
      • Wu Y.-M.
      • Harms P.W.
      • Lin J.
      • Reddy R.
      • Brennan C.
      • Palanisamy N.
      • Chang A.C.
      • Truini A.
      • Truini M.
      • Robinson D.R.
      • Beer D.G.
      • Chinnaiyan A.M.
      Transcriptome meta-analysis of lung cancer reveals recurrent aberrations in NRG1 and Hippo pathway genes.
      and East Asian data sets, were used.
      • Chen J.
      • Yang H.
      • Teo A.S.M.
      • Amer L.B.
      • Sherbaf F.G.
      • Tan C.Q.
      • et al.
      Genomic landscape of lung adenocarcinoma in East Asians.

      Differential Gene Expression and Metascape Analysis

      To compare gene expression levels between the two sample groups, t-test with false discovery rate correction was performed. Metascape (https://metascape.org, last accessed September 1, 2021) analysis was performed to identify the pathways associated with a gene group.

      K-Means Clustering and Gene Set Enrichment Analysis, Cibersort, and Kaplan-Meier Analysis

      Fragments per kilobase of transcript per million mapped reads (FPKM) of senescence genes were normalized, and consensus K-means clustering was performed. Gene set enrichment analysis (GSEA) was performed using default settings. For Kaplan-Meier analysis, Prism software version 8.3 (GraphPad Software, San Diego, CA) was used, and Mantel-Cox statistics were used to calculate hazard ratios and P values. For Cibersort analysis (AlizadehLab, Stanford, CA), RNA-seq data from LUAD samples were used as a mixture file, and the LM22 immune cell type was used as a signature gene file.

      Copy Number Variation and Mutation Analysis

      Data on copy number variations from LUAD samples were downloaded from the Genomic Data Commons data portal and analyzed using Genomic Identification of Significant Targets in Cancer (GISTIC) software version 2.0 (https://www.genepattern.org/modules/docs/GISTIC_2.0, last accessed September 1, 2021). Cytobands with false discovery rates of ≤0.1 were retained, and a t-test using FPKM values was used to compare the genes in these cytobands between the high- and low-senescence groups (HS and LS). Copy number variations at these focal points were visualized using Circos plotting. To compare broad events in each chromosome arm, the Fisher exact test was performed. Values of <−0.3 and >0.3 were considered as deleted and amplified, respectively. Mutation data in mutation annotation format were downloaded from Genomic Data Commons, and all of the nonsilent variants of mutations were retained. The Fisher exact test was performed to compare all of the genes between the HS and LS groups, and a waterfall plot was plotted for selected genes.

      DNA Methylation

      For methylation, TCGA-LUAD Infinium HumanMethylation450 data was used and preprocessed using the method decribed earlier.
      • Maros M.E.
      • Capper D.
      • Jones D.T.W.
      • Hovestadt V.
      • Deimling Avon
      • Pfister S.M.
      • Benner A.
      • Zucknick M.
      • Sill M.
      Machine learning workflows to estimate class probabilities for precision cancer diagnostics on DNA methylation microarray data.
      Differential methylation analysis was performed using a linear model for microarray and RNA-seq data (Limma software version 3.46.0; https://bioconductor.org/packages/release/bioc/html/limma.html, last accessed September 1, 2021)
      • Ritchie M.E.
      • Phipson B.
      • Wu D.
      • Hu Y.
      • Law C.W.
      • Shi W.
      • Smyth G.K.
      limma powers differential expression analyses for RNA-sequencing and microarray studies.
      on quantile-normalized and logarithm-transformed β values.
      • Maros M.E.
      • Capper D.
      • Jones D.T.W.
      • Hovestadt V.
      • Deimling Avon
      • Pfister S.M.
      • Benner A.
      • Zucknick M.
      • Sill M.
      Machine learning workflows to estimate class probabilities for precision cancer diagnostics on DNA methylation microarray data.

      Identification of Senescence Gene Signature Using Machine Learning, and Drug Sensitivity Analysis

      Given that it is vital to retain the features that contribute the most in the distinction of HS samples from LS samples and to eliminate redundant features, feature selection by recursive feature elimination was performed. Caret software version 6.0 (http://topepo.github.io/caret/index.html, last accessed September 1, 2021) library in statistical software package R version 4.0.3 (https://www.R-project.org, last accessed September 1, 2021) was used to perform recursive feature elimination on a random forest model with fivefold cross-validation repeated 100 times. The top 20 features, based on variable importance returned by the random forest model, were considered for further processing. To determine the predictive ability of the chosen features, two random forest models, one with all features (418) and the other with 20 selected features, were built. A train-test split of 80% to 20% (train: C1 = 309, C2 = 106; test: C1 = 77, C2 = 26) was used in both cases. The selected features were validated by the survival analysis in the TCGA, Michigan, and East Asian cohorts. To obtain the senescence groups, K-means clustering was performed in testing cohorts using the expression of the 20 selected genes. Senescence groups were also identified and verified in cell lines using the expression data from the cancer cell line encyclopedia (CCLE). IC50 values for different cell lines from genomics of drug sensitivity in cancer (GDSC) for different drugs were used to compare the sensitivity of drugs.

      Results

      Radiation and Chemotherapy Drugs Induce Different Senescence Genes in Lung Fibroblasts

      RNA-seq data from WI-38 cells (human lung fibroblasts), associated with treatment with doxorubicin and ionizing radiation (IR), were analyzed to identify the changes in gene expression levels during senescence in lung tissue. As expected, both the IR- and doxorubicin-treated cells showed increased expression of senescence marker CDKN1A compared with that in control cells (Figure 1A). The expression changes with IR and doxorubicin treatment were checked. The expression levels of a higher number of genes were changed in IR-treated cells compared with doxorubicin-treated cells (Figure 1B). Metascape analysis showed that genes up-regulated after IR- and doxorubicin-induced senescence modulated the immune system (Supplemental Figure S1A). However, down-regulated genes in IR-induced senescence regulated the cell cycle, whereas down-regulated genes after doxorubicin-induced senescence did not show specific pathway enrichment (Supplemental Figure S1A). A Venn diagram was used to identify 315 commonly up-regulated and 200 commonly down-regulated genes after IR- and doxorubicin-induced senescence (Figure 1C). These genes were identified as senescence-related genes. Next, the expression levels of senescence genes were compared between the TCGA normal and LUAD samples, and differentially expressed genes were used for network analysis (Supplemental Figure S1, B and C). As seen earlier, senescence-induced genes with differential expression in LUAD compared with normal samples showed immune-modulatory function (Supplemental Figure S1D). In comparison, senescence-inhibited genes with differential expression in LUAD compared with normal samples were associated with cell cycle regulation (Figure 1D), suggesting inhibition of the cell cycle after senescence induction.
      Figure thumbnail gr1
      Figure 1Identification of senescence subtype of LUAD. A: The expression of CDKN1A after ionizing radiation (IR) and doxorubicin (Dox) treatment in WI-38 cells. B: Changes in expression levels of protein-coding genes (PcG) and long noncoding RNAs (LncRNA) in WI-38 cells. C: Overlap of up- and down-regulated genes. D: The senescence inhibited and differentially expressed genes in LUAD compared with normal are used for the Metascape analysis, and significant terms are used to make the Gene Ontology network. E: Heatmap of the K-means clustering results obtained using LUAD-associated senescence genes. F: Kaplan-Meier plot of the significant difference in survival between C1 and C2 patients. Data are expressed as means (SD) (A). HR, hazard ratio; MS, median survival; Trea, treated; Unt, untreated.

      Senescence-Associated Genes Identify Two Novel Subclasses of LUAD with the Difference in Senescence Level

      To identify tumors with different senescence statuses, K-means consensus clustering was performed using senescence genes differentially expressed in LUAD compared with normal samples. The analysis identified two clusters with variable expression of LUAD-associated senescence genes (Figure 1E and Supplemental Figure S1E). Heatmapping of senescence genes showed that senescence cluster C2 had higher expression of genes repressed during senescence and a lower expression of genes induced during senescence (Supplemental Figure S1E). In contrast, senescence cluster C1 showed a lower expression of senescence-repressed genes and a higher expression of senescence-induced genes (Supplemental Figure S1E). Hence, clusters C1 and C2 were named as the HS and LS groups, respectively (Supplemental Figure S1E). Additionally, the HS patients showed a significantly poor survival rate compared with that in the LS patients (Figure 1F).

      High- and Low-Senescence Groups Have Different Immune and Genetic Landscapes

      The senescence nature in the HS and LS groups was confirmed using senescence-associated gene sets in GSEA (Figure 2A). GSEA with oncogenic signature identified positive enrichment of various oncogenic pathways in the LS group (Figure 2B). Immune microenvironment plays a key role in cancer development and progression. Also, senescent cells induce the release of a variety of factors, called senescent-associated secretory phenotypes, to attract immune cells to remove senescent cells from the tissue environment.
      • Zhang B.
      • Fu D.
      • Xu Q.
      • Cong X.
      • Wu C.
      • Zhong X.
      • Ma Y.
      • Lv Z.
      • Chen F.
      • Han L.
      • Qian M.
      • Chin Y.E.
      • Lam E.W.-F.
      • Chiao P.
      • Sun Y.
      The senescence-associated secretory phenotype is potentiated by feedforward regulatory mechanisms involving Zscan4 and TAK1.
      ,
      • Coppé J.-P.
      • Desprez P.-Y.
      • Krtolica A.
      • Campisi J.
      The senescence-associated secretory phenotype: the dark side of tumor suppression.
      Therefore, the enrichment of various immune cells in the HS and LS groups was estimated. Cibersort analysis showed that T cells (CD4, CD8, and T-regulatory cells) and neutrophils were significantly more enriched in the HS group, whereas B cells, resting mast cells, and monocytes were differentially infiltrated in the LS group (Supplemental Figure S2). In a recent landmark paper, Thorsson et al
      • Thorsson V.
      • Gibbs D.L.
      • Brown S.D.
      • Wolf D.
      • Bortone D.S.
      • Yang T.-H.O.
      • et al.
      The immune landscape of cancer.
      identified the immune subtype of various cancers, including LUAD. In keeping with that, LS patients had significantly enriched C1 and C2 classes of the immune subtypes (Figure 2C). In contrast, HS patients had significantly enriched C3 subtype (Figure 2C) and higher lymphocyte infiltration (Figure 2D). As expected, the LS group showed significantly higher proliferation, wound healing, single-nucleotide variants, indel variations, and nonsilent mutation rate (Figure 2, E–H).
      Figure thumbnail gr2
      Figure 2Characterization of senescence subtypes of LUAD. A: Enrichment plots of the senescence-associated gene sets. B: Significant gene sets identified in gene set enrichment analysis using oncogenic signature gene sets, showing the enrichment of higher activation of oncogenic signaling in the LS group. C: Percents of HS and LS groups in immune subgroups. D: Lymphocyte fractions in the HS and LS groups. E and F: Proliferation and wound-healing scores in the HS and LS groups. G and H: Copy number variations, indels, and nonsilent mutations in the LS and HS groups. Data are expressed as first quartile, median, third quartile and interquartile ranges (D–H). HS versus LS: ∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001, and ∗∗∗∗P < 0.0001. E2F, E2 factor; FDR, false discovery rate; FRC_EZH2, fibroblastic reticular cell–encoding enhancer of zeste homolog 2; GO, Gene Ontology; mTOR, mammalian target of rapamycin; NES, normalized enrichment score; SNV, single-nucleotide variant; VEGF, vascular endothelial growth factor.
      Given that the HS group had overall less genetic variation, the specific genetic differences between the two groups were compared. The TP53 was the most significantly differentially mutated gene between the HS and LS groups (Figure 3A). Intriguingly, the LS group had higher mutation rates of CAD, GRIA4, and SERPINA3 (Figure 3A).
      Figure thumbnail gr3
      Figure 3Changes in DNA methylation in the HS and LS groups: A: Waterfall plot showing mutations of various genes in the HS and LS groups. B: Circos plot showing copy number variations (CNVs) in HS (outer circle) and LS (inner circle) groups. C: β values in patients in the normal lung, HS, and LS groups. D: CpGs subclassified based on the location and number of CpGs having negative and positive correlation with expression in the LS and HS groups. E: The genes silenced due to hypermethylation in the HS group are used as input for the Metascape analysis, and significant Gene Ontology (GO) terms are plotted. Data are expressed as first quartile, median, third quartile, and interquartile ranges (C). ∗P < 0.05 and ∗∗∗∗P < 0.0001. NS, not significant; UTR, untranslated region.
      Copy number variations were also less common in the LS group. Broad significant differences in copy number variation were seen in chromosomes 3p, 4p, 4q, 5p, 5q, 11p, 11q, 13q, 15q, 16p, 16q, 19p, 19q, 21q, and 22q (Figure 3B). As expected, the number of focal variations was higher in the LS group (Figure 3B). The LS group also had higher amplification of many important loci harboring genes such as CDC45, SOX2, TERT, and MET. The genes amplified and overexpressed in LS were associated with cell cycle promotion and proliferation (Supplemental Figure S3A and Supplemental Table S1). In comparison, the HS group showed amplification and overexpression of genes associated with cellular transport (Supplemental Figure S3A and Supplemental Table S1). Similarly, the LS group showed higher focal deletion and under-expression of genes associated with cellular transport and cellular movements (Supplemental Figure S3B and Supplemental Table S1). In contrast, the HS group showed deletion and down-regulation of genes associated with cell cycle and proliferation (Supplemental Figure S3B and Supplemental Table S1). These findings may explain the high-proliferation phenotype and poor survival in LS patients.

      High- and Low-Senescence Groups Have Significantly Different Epigenetic Alterations

      Epigenetic changes play a crucial role in age-related phenomena such as senescence and cancer. To understand the DNA methylation changes in the HS and LS groups, TCGA Infinium Human Methylation 450 data were analyzed. Interestingly, a significantly lower overall methylation level was observed in the LS group compared with that in the HS group (Figure 3C). Also, the overall methylation levels in the normal lung and HS groups were similar (Figure 3C), suggesting the role of DNA methylation in senescence phenotype in LUAD. Next, differential methylation analysis identified 596 hypermethylated and 685 hypomethylated CpGs (Supplemental Figure S4). The correlation analysis identified 66 CpGs (61 genes) with a significant negative and 48 CpGs (41 genes) with a significant positive correlation between methylation and expression (Supplemental Table S2). The CpGs present on CpG islands at the promoters (promoter, 5′-untranslated region, first exon) generally repress the expression, and CpGs from the gene body are associated with activation.
      • Jjingo D.
      • Conley A.B.
      • Yi S.V.
      • Lunyak V.V.
      • Jordan I.K.
      On the presence and role of human gene-body DNA methylation.
      ,
      • Arechederra M.
      • Daian F.
      • Yim A.
      • Bazai S.K.
      • Richelme S.
      • Dono R.
      • Saurin A.J.
      • Habermann B.H.
      • Maina F.
      Hypermethylation of gene body CpG islands predicts high dosage of functional oncogenes in liver cancer.
      Hence, the CpGs with significant correlation between methylation and expression were classified as island and nonisland CpGs. As expected, CpGs present on promoter islands were primarily associated with repression, and nonisland CpGs were associated with the activation of gene expression (Figure 3D). Furthermore, Metascape analysis showed that the genes suppressed by DNA methylation in HS were enriched in the Wnt pathway, development-related pathways, and senescence-regulating pathways (Figure 3E). This finding validates previous reports suggesting the role of Wnt pathway inhibition in senescence.
      • Jjingo D.
      • Conley A.B.
      • Yi S.V.
      • Lunyak V.V.
      • Jordan I.K.
      On the presence and role of human gene-body DNA methylation.
      ,
      • Arechederra M.
      • Daian F.
      • Yim A.
      • Bazai S.K.
      • Richelme S.
      • Dono R.
      • Saurin A.J.
      • Habermann B.H.
      • Maina F.
      Hypermethylation of gene body CpG islands predicts high dosage of functional oncogenes in liver cancer.

      Machine Learning–Based Validation of Senescence Phenotype with Drug Resistance in Independent Sample Sets

      To identify a senescence signature, feature selection was performed by recursive feature elimination. The analysis identified 20 genes the expression of which could be used for identifying the HS and LS groups without significant accuracy loss (Supplemental Figure S5, A–C). To validate that these 20 genes were sufficient for identifying the senescence groups in LUAD, K-means consensus clustering was performed in an independent cohort (East Asian cohort) of LUAD patients. This analysis identified two clusters, named C1 and C2 (Figure 4A). The C1 group had a high expression levels of 20 senescence genes compared with those in the C2 group (Supplemental Figure S5D). Hence, C1 was labeled as the HS group and C2 as the LS group. Subsequent Kaplan-Meier analysis showed that survival was significantly greater in the HS patients, validating the findings from the TCGA data set (Figure 4B). Furthermore, the HS status of the C1 group was confirmed using GSEA (Supplemental Figure S5E). The HS group also showed lower activation of the retinoblastoma tumor suppressor (RB), polycomb repressive complex (PRC)-2–enhancer of zeste homolog (EZH)-2, and mammalian target of rapamycin (mTOR) signaling pathways (Supplemental Figure S5F). As in the East Asian cohort, the presence of the senescence phenotype was validated in the Michigan cohort. Although the difference in survival between the HS and LS groups was nonsignificant, the trend remained the same (Supplemental Figure S6, A–C).
      Figure thumbnail gr4
      Figure 4Machine-learning model of senescence signature generation. A: The heatmap showing the K-means clustering results identified using 20 senescence signature genes. B: Kaplan-Meier plot of the difference in survival between the HS and LS groups from the East Asian-lung adenocarcinoma data set. C: K-means consensus cluster of Cancer Cell Line Encyclopedia (CCLE)-Lung cells data set, as identified using 20 senescence-associated genes. D: Heatmap of the expression of 20 selected genes in CCLE-lung cells. E: IC50 values of all of the drugs in lung cancer cells from the Genomics of Drug Sensitivity in Cancer database were used for classifying HS and LS cells. F and G: IC50 values of histone deacetylase (HDAC) inhibitors (belinostat, trichostatin A, CAY10603, and AR42) and cyclin-dependent kinase (CDK9) inhibitors (pictilisib and lestaurtinib) in lung cells classified as HS (yellow) or LS (blue). Data are expressed as first quartile, median, third quartile and interquartile ranges (E–G). HS versus LS: ∗P < 0.05, ∗∗P < 0.01, and ∗∗∗∗P < 0.0001. HR, hazard ratio; MS, median survival; NS, not significant.
      Next, cell lines from the Cancer Cell Line Encyclopedia (CCLE) data set were classified into HS (C1) and LS (C2) groups using the expression of 20-gene senescence signatures (Figure 4, C and D). The GSEA validated the senescence status of the cell lines in the HS and LS groups (Supplemental Figure S6D). Chemosensitivity data from the genomics of the Genomics of Drug Sensitivity in Cancer database were analyzed to understand drug sensitivity in senescence groups. As expected, HS cells showed significantly higher IC50 values for all drugs combined (Figure 4E). Further investigation of the IC50 values of individual drugs revealed that the histone deacetylase (HDAC) inhibitors belinostat, trichostatin A, CAY10603, and AR42 showed significantly higher IC50 in the HS group compared with the LS group (Figure 4F). Similarly, the LS group was considerably more sensitive to various cyclin-dependent kinase (CDK)-9 inhibitors (Figure 4G). These observations suggest that the senescence nature of tumor is an essential determinant of therapy decision and disease outcome.

      Discussion

      Cancer, a leading cause of fatalities across the world, is characterized by its strong proliferative signal. A reduction in tumor growth requires multiple doses of chemotherapy. Over the years, cellular senescence has been identified as another way to push cancer cells into permanent proliferation arrest.
      • He S.
      • Sharpless N.E.
      Senescence in health and disease.
      Early reports suggest that cell senescence varies depending on the cell origin.
      • Zhang B.
      • Fu D.
      • Xu Q.
      • Cong X.
      • Wu C.
      • Zhong X.
      • Ma Y.
      • Lv Z.
      • Chen F.
      • Han L.
      • Qian M.
      • Chin Y.E.
      • Lam E.W.-F.
      • Chiao P.
      • Sun Y.
      The senescence-associated secretory phenotype is potentiated by feedforward regulatory mechanisms involving Zscan4 and TAK1.
      Therefore, lung fibroblast cells treated with IR or doxorubicin were utilized here to identify the senescence genes in lung tissue.
      • Zhang B.
      • Fu D.
      • Xu Q.
      • Cong X.
      • Wu C.
      • Zhong X.
      • Ma Y.
      • Lv Z.
      • Chen F.
      • Han L.
      • Qian M.
      • Chin Y.E.
      • Lam E.W.-F.
      • Chiao P.
      • Sun Y.
      The senescence-associated secretory phenotype is potentiated by feedforward regulatory mechanisms involving Zscan4 and TAK1.
      As compared with doxorubicin-regulated genes, IR-regulated genes were highly associated with cell cycle regulation, suggesting a greater impact of IR than of chemotherapy on molecular outcomes. IR- and doxorubicin-regulated genes were used to identify HS and LS groups of patients. HS patients showed better survival and lower activation of proliferative signaling, suggesting that tumors with an HS nature are inherently less aggressive compared with tumors with an LS nature. A lower rate of mutation of the TP53 gene could be involved in the HS nature of HS-group tumors. Additionally, HS patients showed lower genetic aberration, which also supports the lesser aggressiveness of these tumors.
      Many reports have suggested that DNA methylation patterns of senescence are a crucial determinant of tumor behavior and hence are associated with cancer and aging.
      • Xie W.
      • Kagiampakis I.
      • Pan L.
      • Zhang Y.W.
      • Murphy L.
      • Tao Y.
      • Kong X.
      • Kang B.
      • Xia L.
      • Carvalho F.L.F.
      • Sen S.
      • Yen R.-W.C.
      • Zahnow C.A.
      • Ahuja N.
      • Baylin S.B.
      • Easwaran H.
      DNA methylation patterns separate senescence from transformation potential and indicate cancer risk.
      ,
      • Cruickshanks H.A.
      • McBryan T.
      • Nelson D.M.
      • VanderKraats N.D.
      • Shah P.P.
      • van Tuyn J.
      • Rai T.S.
      • Brock C.
      • Donahue G.
      • Dunican D.S.
      • Drotar M.E.
      • Meehan R.R.
      • Edwards J.R.
      • Berger S.L.
      • Adams P.D.
      Senescent cells harbour features of the cancer epigenome.
      In the current study, LS patients had significantly greater genome-wide hypomethylation than that of HS patients and the methylation-silenced genes in the HS group were enriched in Wnt signaling. The down-regulation of Wnt signaling is crucial for the onset of senescence.
      • Adams P.D.
      • Enders G.H.
      Wnt signaling and senescence: a tug of war in early neoplasia?.
      • Jeoung J.Y.
      • Nam H.Y.
      • Kwak J.
      • Jin H.J.
      • Lee H.J.
      • Lee B.W.
      • Baek J.H.
      • Eom J.S.
      • Chang E.-J.
      • Shin D.-M.
      • Choi S.J.
      • Kim S.W.
      A decline in Wnt3a signaling is necessary for mesenchymal stem cells to proceed to replicative senescence.
      • Ye X.
      • Zerlanko B.
      • Kennedy A.
      • Banumathy G.
      • Zhang R.
      • Adams P.D.
      Downregulation of Wnt signaling is a trigger for formation of facultative heterochromatin and onset of cell senescence in primary human cells.
      The data suggest that epigenetic silencing of Wnt signaling is an important factor in the determination of HS phenotype in lung cancer. The presence of senescence phenotype was also validated in two independent cohorts using the senescence gene signature as identified by machine learning–based methods. The patients with HS nature showed better survival and lower activation of proliferative pathways as found in the TCGA data set.
      To summarize, this study used RNA-seq data from senescent lung fibroblasts to identify a novel senescence phenotype in LUAD patients. The patients with HS showed a better prognosis and poor activation of mTOR. Additionally, The study used a machine-learning model to retain the minimal number of genes required to identify the senescent phenotype, and further validated it in two independent data sets of LUAD patients and CCLE cell lines. The study also showed that cell lines from the HS group were significantly more resistant to drugs.

      Acknowledgment

      We thank The Cancer Genome Atlas for providing patient data.

      Supplemental Data

      • Supplemental Figure S1

        A: Pathway analysis using the up- and down-regulated genes after ionizing radiation (IR)- and doxorubicin (Dox)-induced senescence, and significantly enriched Gene Ontology terms. The dashed line indicates the significance cutoff [false discovery rate (FDR) of <0.05]. B and C: Expression of genes commonly up-regulated (B) and down-regulated (C) after IR- and Dox-induced senescence, compared in LUAD and normal lung, and FDR and logarithm twofold changes. The dashed lines indicate the significance cutoff. D: Senescence-induced and differentially expressed genes in LUAD compared to normal tissue in the Metascape analysis, and significant terms. E: Heatmap of LUAD-associated senescence genes in C1 and C2 groups as identified by K-means clustering.

      • Supplemental Figure S3

        A: Genes differentially amplified in HS and LS groups identified by analysis of copy number variations. The genes amplified and overexpressed in the LS group are used as input for Metascape analysis, with the top 25 enriched Gene Ontology (GO) terms. Similarly, genes amplified and overexpressed in the HS group are used as input for Metascape analysis, with enriched GO terms. B: Genes differentially deleted in HS and LS groups identified by copy number variation analysis. The genes deleted and under-expressed in the LS group are used as input for Metascape analysis, with the top 25 enriched GO terms. Similarly, genes deleted and under-expressed in the HS group are used as input for Metascape analysis, with the top 25 enriched GO terms.

      • Supplemental Figure S5

        A: Receiver operating characteristic analysis comparing the sensitivity and specificity of gene sets in identifying HS and LS groups in The Cancer Genome Atlas (TCGA) data set, showing the areas under the curve (AUC). B: K-means consensus clustering of TCGA patients based on the expression of 20 genes as identified by random forest analysis (RFA). C and D: Heatmaps of the expression of 20 genes as identified by RFA in TCGA samples (C) and in EAS-LUAD (D) HS and LS groups. E and F: The significant gene sets identified in gene set enrichment analysis (GSEA) (E) and EAS-LUAD (F) data using oncogenic signature gene sets are plotted to show the enrichment of higher activation of oncogenic signaling in the HS and LS groups. EZH, enhancer of zeste homolog; FDR, false discovery rate; GO, Gene Ontology; mTOR, mammalian target of rapamycin; NES, normalized enrichment score; PRC, polycomb repressive complex; RB, retinoblastoma tumor suppressor; RF, random forest.

      • Supplemental Figure S6

        A: K-means consensus clustering of the Michigan cohort based on the expression of 20 genes as identified by random forest analysis (RFA). B: Heatmap showing the expression of 20 genes as identified by RFA in Michigan samples. C: K-means plot of the difference in survival between the HS and LS groups as identified in K-means clustering. The dashed lines indicate median survival. D: Enrichment of the senescence-associated gene sets, showing the enrichment of senescence-associated genes in the HS group. GO, Gene Ontology.

      References

        • Bade B.C.
        • Cruz C.S.D.
        Lung cancer 2020 epidemiology, etiology, and prevention.
        Clin Chest Med. 2020; 41: 1-24
        • Duma N.
        • Santana-Davila R.
        • Molina J.R.
        Non–small cell lung cancer: epidemiology, screening, diagnosis, and treatment.
        Mayo Clin Proc. 2019; 94: 1623-1640
        • Kim Y.H.
        • Choi Y.W.
        • Lee J.
        • Soh E.Y.
        • Kim J.-H.
        • Park T.J.
        Senescent tumor cells lead the collective invasion in thyroid cancer.
        Nat Commun. 2017; 8: 15208
        • He S.
        • Sharpless N.E.
        Senescence in health and disease.
        Cell. 2017; 169: 1000-1011
        • Wang B.
        • Kohli J.
        • Demaria M.
        Senescent cells in cancer therapy: friends or foes?.
        Trends Cancer. 2020; 6: 838-857
        • Casella G.
        • Munk R.
        • Kim K.M.
        • Piao Y.
        • De S.
        • Abdelmohsen K.
        • Gorospe M.
        Transcriptome signature of cellular senescence.
        Nucleic Acids Res. 2019; 47: 7294-7305
        • Shukla S.
        • Evans J.R.
        • Malik R.
        • Feng F.Y.
        • Dhanasekaran S.M.
        • Cao X.
        • Chen G.
        • Beer D.G.
        • Jiang H.
        • Chinnaiyan A.M.
        Development of a RNA-seq based prognostic signature in lung adenocarcinoma.
        J Natl Cancer Inst. 2017; 109: djw200
        • Dhanasekaran S.M.
        • Balbin O.A.
        • Chen G.
        • Nadal E.
        • Kalyana-Sundaram S.
        • Pan J.
        • Veeneman B.
        • Cao X.
        • Malik R.
        • Vats P.
        • Wang R.
        • Huang S.
        • Zhong J.
        • Jing X.
        • Iyer M.
        • Wu Y.-M.
        • Harms P.W.
        • Lin J.
        • Reddy R.
        • Brennan C.
        • Palanisamy N.
        • Chang A.C.
        • Truini A.
        • Truini M.
        • Robinson D.R.
        • Beer D.G.
        • Chinnaiyan A.M.
        Transcriptome meta-analysis of lung cancer reveals recurrent aberrations in NRG1 and Hippo pathway genes.
        Nat Commun. 2014; 5: 5893
        • Chen J.
        • Yang H.
        • Teo A.S.M.
        • Amer L.B.
        • Sherbaf F.G.
        • Tan C.Q.
        • et al.
        Genomic landscape of lung adenocarcinoma in East Asians.
        Nat Genet. 2020; 52: 177-186
        • Maros M.E.
        • Capper D.
        • Jones D.T.W.
        • Hovestadt V.
        • Deimling Avon
        • Pfister S.M.
        • Benner A.
        • Zucknick M.
        • Sill M.
        Machine learning workflows to estimate class probabilities for precision cancer diagnostics on DNA methylation microarray data.
        Nat Protoc. 2020; 15: 479-512
        • Ritchie M.E.
        • Phipson B.
        • Wu D.
        • Hu Y.
        • Law C.W.
        • Shi W.
        • Smyth G.K.
        limma powers differential expression analyses for RNA-sequencing and microarray studies.
        Nucleic Acids Res. 2015; 43: e47
        • Zhang B.
        • Fu D.
        • Xu Q.
        • Cong X.
        • Wu C.
        • Zhong X.
        • Ma Y.
        • Lv Z.
        • Chen F.
        • Han L.
        • Qian M.
        • Chin Y.E.
        • Lam E.W.-F.
        • Chiao P.
        • Sun Y.
        The senescence-associated secretory phenotype is potentiated by feedforward regulatory mechanisms involving Zscan4 and TAK1.
        Nat Commun. 2018; 9: 1723
        • Coppé J.-P.
        • Desprez P.-Y.
        • Krtolica A.
        • Campisi J.
        The senescence-associated secretory phenotype: the dark side of tumor suppression.
        Pathol Mech Dis. 2010; 5: 99-118
        • Thorsson V.
        • Gibbs D.L.
        • Brown S.D.
        • Wolf D.
        • Bortone D.S.
        • Yang T.-H.O.
        • et al.
        The immune landscape of cancer.
        Immunity. 2018; 48: 812-830.e14
        • Jjingo D.
        • Conley A.B.
        • Yi S.V.
        • Lunyak V.V.
        • Jordan I.K.
        On the presence and role of human gene-body DNA methylation.
        Oncotarget. 2012; 3: 462-474
        • Arechederra M.
        • Daian F.
        • Yim A.
        • Bazai S.K.
        • Richelme S.
        • Dono R.
        • Saurin A.J.
        • Habermann B.H.
        • Maina F.
        Hypermethylation of gene body CpG islands predicts high dosage of functional oncogenes in liver cancer.
        Nat Commun. 2018; 9: 3164
        • Xie W.
        • Kagiampakis I.
        • Pan L.
        • Zhang Y.W.
        • Murphy L.
        • Tao Y.
        • Kong X.
        • Kang B.
        • Xia L.
        • Carvalho F.L.F.
        • Sen S.
        • Yen R.-W.C.
        • Zahnow C.A.
        • Ahuja N.
        • Baylin S.B.
        • Easwaran H.
        DNA methylation patterns separate senescence from transformation potential and indicate cancer risk.
        Cancer Cell. 2018; 33: 309-321.e5
        • Cruickshanks H.A.
        • McBryan T.
        • Nelson D.M.
        • VanderKraats N.D.
        • Shah P.P.
        • van Tuyn J.
        • Rai T.S.
        • Brock C.
        • Donahue G.
        • Dunican D.S.
        • Drotar M.E.
        • Meehan R.R.
        • Edwards J.R.
        • Berger S.L.
        • Adams P.D.
        Senescent cells harbour features of the cancer epigenome.
        Nat Cell Biol. 2013; 15: 1495-1506
        • Adams P.D.
        • Enders G.H.
        Wnt signaling and senescence: a tug of war in early neoplasia?.
        Cancer Biol Ther. 2008; 7: 1706-1711
        • Jeoung J.Y.
        • Nam H.Y.
        • Kwak J.
        • Jin H.J.
        • Lee H.J.
        • Lee B.W.
        • Baek J.H.
        • Eom J.S.
        • Chang E.-J.
        • Shin D.-M.
        • Choi S.J.
        • Kim S.W.
        A decline in Wnt3a signaling is necessary for mesenchymal stem cells to proceed to replicative senescence.
        Stem Cells Dev. 2015; 24: 973-982
        • Ye X.
        • Zerlanko B.
        • Kennedy A.
        • Banumathy G.
        • Zhang R.
        • Adams P.D.
        Downregulation of Wnt signaling is a trigger for formation of facultative heterochromatin and onset of cell senescence in primary human cells.
        Mol Cell. 2007; 27: 183-196