Screening common signaling pathways associated with drug resistance in non‐small cell lung cancer via gene expression profile analysis

Abstract Lung cancer is the leading cause of cancer‐related deaths worldwide. Although several therapeutic strategies have been employed to curb lung cancer, the survival rate is still poor owing to the development of drug resistance. The mechanisms underlying drug resistance development are incompletely understood. Here, we aimed to identify the common signaling pathways involved in drug resistance in non‐small cell lung cancer (NSCLC). Three published transcriptome microarray data were downloaded from the Gene Expression Omnibus (GEO) database comprising different drug‐resistant cell lines and their parental cell lines. Differentially expressed genes (DEGs) were identified and used to perform Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. An overlapping analysis was performed for KEGG pathways enriched from all the three datasets to identify the common signaling pathways. As a result, we found that metabolic pathways, ubiquitin‐mediated proteolysis, and mitogen‐activated protein kinase (MAPK) signaling were the most aberrantly expressed signaling pathways. The knockdown of nicotinamide phosphoribosyltransferase (NAMPT), the gene involved in metabolic pathways and known to be upregulated in drug‐resistant tumor cells, was shown to increase the apoptosis of cisplatin‐resistant A549 cells following cisplatin treatment. Thus, our results provide an in‐depth analysis of the signaling pathways that are commonly altered in drug‐resistant NSCLC cell lines and highlight the potential strategy that facilitates the development of interventions to interfere with upregulated signaling pathways as well as to boost downregulated signaling pathways in drug‐resistant tumors for the elimination of multiple resistance of NSCLC.


| INTRODUCTION
Lung cancer is the most common and lethal malignancy, with a 5-year survival rate of only 16%. 1 Non-small cell lung cancer (NSCLC) is the most common type of lung cancer that accounts for more than 80% of all lung cancer cases. 2,3 Chemotherapeutic drugs such as gemcitabine, taxane, and platin, alone or in combination, are the first line of treatment for patients with NSCLC. In addition, new therapeutic strategies have been applied for lung cancer treatment, including immunotherapy. [4][5][6] Immunotherapy using antibodies that block immune checkpoints, programmed cell death 1 (PD-1)/programmed cell death 1 ligand (PD-L1), has shown impressive antitumor effects and clinical benefits for the treatment of NSCLC. The Food and Drug Administration (FDA) has approved anti-PD-1 antibodies, nivolumab, and pembrolizumab, for treatment of solid tumors including advanced NSCLC. [7][8][9] The China Food and Drug Administration (CFDA) has recently approved nivolumab for the treatment of advanced NSCLC that lacks epidermal growth factor receptor (EGFR) mutations and anaplastic lymphoma kinase (ALK) mutations. However, most patients with late-stage NSCLC show resistance to these drugs that has led to an increase in mortality rates. 10 Drug resistance is a gene-driven and pathway-mediated process. Gene heterogeneities and mutations have been shown to play an important role in influencing drug efficacy and resistance in lung cancer; these phenomenon include mutations in tumor suppressor protein 53 (TP53), 6,11 kirsten rat sarcoma viral oncogene (KRAS), 12,13 PI3-kinase subunit alpha (PIK3CA), PI3-kinase subunit delta (PIK3CD), 14,15 and EGFR. 16 In addition, noncoding RNAs (eg, microRNAs [miRNAs]) are known to alter drug resistance or sensitivity in lung cancer under chemotherapy treatment. 1 These molecular mechanisms cause drug efflux, drug metabolism/inactivation, drug target alteration, DNA repair, or apoptosis deficiency. Over the past decades, therapies targeting either gene mutation or amplification have offered promising results, but the acquired resistance following treatment has demanded further investigations. It is acknowledged that genetic alterations in signaling pathways controlling apoptosis, cell cycle, and cell growth are common hallmarks of drug-resistant tumors. Alterations in signaling pathways vary with different drug treatments. The drugs that target key molecules or block the downstream signaling pathways may cause feedback regulations and activate alternative signaling pathways to resist drug-mediated destruction, as observed in tumor cells; this phenomenon has become a major obstacle for the elimination of drug resistance. Hence, it is essential to identify the common signaling pathways that mediate drug resistance.
Microarray technology has recently gained popularity for the investigation of gene alterations in tumorigenesis, metastasis, cancer recurrence, and drug resistance as well as to identify biomarkers for tumor diagnosis, prognosis, and therapy. [17][18][19][20][21][22] Through RNA-sequencing analysis, many genes, RNAs, including messenger RNAs (mRNAs), longnoncoding RNAs (lncRNAs), and miRNAs, and proteins have been reported to play a vital role in lung cancer initiation, progression, and recurrence. Larsen et al 23 identified a distinct gene expression profile for recurrent lung squamous cell carcinoma through genome-wide profiling. Sun et al 24 reported the decrease in the expression of insulin-like growth factor-binding protein 3 (IGFBP3) in cisplatin-resistant lung cancer cells, resulting in an increase in the activation of insulin-like growth factor 1-receptor (IGF-1R) signaling and the subsequent resistance to cisplatin and radiation. Therefore, advances in RNA technology have made it possible to systematically study the genetic changes and commonly altered signaling pathways in drug-resistant tumor cells.
In this study, we screened the Gene Expression Omnibus (GEO) database and analyzed the selected three datasets with a multistep strategy ( Figure 1). We first compared the differences at the gene level between drug-resistant tumor cell lines and parental cell lines and identified several differentially expressed genes (DEGs). We performed Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of the DEGs for each dataset. An overlapping analysis of the KEGG pathways enriched from the three datasets was subsequently performed and 17 common pathways were identified. We detected the common signaling pathways associated with drug resistance and provided potential therapeutic options and possibilities to overcome drug resistance in patients with lung cancer.

| Microarray data
The transcriptome profiles of GSE6914 and GSE77209 were downloaded from GEO (https ://www.ncbi.nlm.nih.gov/geo/) database. A total of 24 samples were used, including 10 drugresistant NSCLC cell lines and 14 parental cell lines. GSE6914 comprised one dataset, Calu3, containing four gemcitabineresistant Calu3 NSCLC cell lines and four parental cell lines. These eight samples were accessed using GPL96 Affymetrix Human Genome U133A Array. GS E77209 comprised two datasets, H1299T18 and H1355T16. The H1299T18 dataset had three taxane-platin-resistant H1299T NSCLC cell lines and five parental cell lines, while the H1355T16 dataset was composed of three taxane-pla tin-resistant H1355 NSCLC cell lines and five parental cell lines. These 16 samples were accessed using GPL10058 Illu mina HumanHT-12 v4.0 expression BeadChip.

| Data preprocessing
Normalization is very crucial for the comparison of microarray datasets. We used the GEOquery package to download the raw probe-level data (CEL. files) through the GEO database. The Z-score method was used for the normalization of the three microarray data. Values of multiple probes that corresponded to the same gene were averaged as the final expression level for the specific gene.

| Screening and clustering of the DEGs
The limma package was used for the identification of the DEGs between drug-resistant cell lines and parental cell lines. The random variance model (RVM) t test was applied to filter the DEGs between the two groups, and the cut-off value was set as |log (fold change)| >1.2 and false discovery rate (FDR) <0.05. Hierarchical clustering of the DEGs was based on the Euclidean distance, and was performed with EPCLUST. [25][26][27] Venn diagram package was used to perform Venn analysis of the DEGs in three datasets. Unique DEGs were selected.

| Enrichment analysis of unique DEGs
The GO analysis was used to analyze the biological functions of the genes, while KEGG pathway enrichment analysis was performed to investigate the signaling pathways that were related to the unique DEGs. The bioconductor package was used to perform GO and KEGG enrichment analyses. In particular, two-sided Fisher's exact and chi-squared tests were used to classify the GO category, FDR and q values were calculated to correct the P-values. [28][29][30][31][32] For the overlapping KEGG pathways, we overlapped the enriched KEGG pathways from the three datasets and identified the common KEGG pathways that were consistently altered in three datasets.

| RNA extraction and real-time polymerase chain reaction (RT-PCR)
The total RNA from cisplatin-resistant A549 cells was extracted with Trizol (Invitrogen Corporation, Carlsbad, CA, #A33250). NanoDrop 2000 (Thermo Scientific, Waltham, MA) was used to detect the concentration and purity of total RNA. The 1 μg of total RNA was reversely transcribed into complementary DNA (cDNA) with Prime Script RT reagent kit with genomic DNA eraser (TaKaRa, Tokyo, Japan, #RR037A). The primers used for cDNA preparation were oligo(dT). Real-time quantitative PCR was performed on Agilent Mx3005P (Santa Clara, CA) consisting of specific primers and SYBR Premix Ex Taq II (TaKaRa, Tokyo, Japan, #RR8202). Primers used for this experiment are listed in Supplementary Table S10. All primers were purchased from Sangon Biotech (Shanghai, China), and the PCR program was as follows: 95°C for 30 s, 95°C for 5 s, and 60°C for 30 s. The range of CT values in this experiment was 21.25-30.53. We used the 2 −ΔΔCt method for the analysis of the relative gene expression level. Each experiment was independently performed in triplicates, and glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used for the normalization of data.

| RNA interference
For transient expression knockdown, small-interfering RNA (siRNA)-NAMPT and negative control were designed and purchased from Gene Pharma Company (Shanghai, China); the siRNA sequences are shown in Supplementary Table  S11. The siRNAs (si-1-NAMPT-349, si-2-NAMPT-1757, and si-control-NAMPT) were diluted in diethyl pyrocarbonate water at a final concentration of 20 μM according to the manufacturer's description. A549 and cisR-A549 cells were plated in six-well plates at a density of 5 × 10 5 cells per well. After reaching approximately 60% confluency, the cells were transfected with the siRNA or negative control (2 μM/ well) using Lipofectamine 3000 (Invitrogen, Carlsbad, CA, #L3000001) for 48 hours. After transfection, the cells were collected for RT-PCR analysis or treated with cisplatin.

| Flow cytometry for apoptosis analysis
After 24 hours of treatment with cisplatin, the supernatant was harvested and tumor cells were washed twice with ice-cold phosphate-buffered saline (PBS). The cells were suspended in Annexin-V-binding buffer (BioLegend, San Diego, CA, #422201) at a final concentration of 10 6 cells/mL, and 100 μL of cell suspension was transferred into a 1.5 mL centrifuge tube and incubated with 5 μL of Alexa Fluor 647 Annexin-V fluorescein isothiocyanate (FITC; BioLegend, San Diego, CA, #640906) and RNase (Thermo Fisher Scientific, USA, #AM2286) for 15 min at 4°C in the dark. After incubation, the samples were treated with propidium iodide (5 μL) (Sigma, Santa Clara, CA, USA, #P4170) and immediately analyzed with flow cytometry (BD, San Diego, CA, USA, FACSCanto II).

| Statistical analysis
Data were analyzed with SPSS software (version 16.0; SPSS Inc, Chicago, IL) and GraphPad Prism 7 (version 5.0; San Diego, CA), and the results were presented as mean ± SD. The ratio data were subjected to log2 transformation and the proportion data were subjected to logit transformation. The two-tailed unpaired t test was used to compare the difference in the relative gene expression and apoptosis ratios between experimental and control groups. The results were presented as histograms with overlaid dot plots; the whiskers represented error bars, and the upper box boundaries represented an average value. The dots represented the mean values of two technical repetitions. Each experiment has at least three biological replicates. P-value less than 0.05 was considered statistically significant.

| DEGs identified between drugresistant cell lines and control parental cell lines
The transcriptome profiles of GSE6914 and GSE77209 were downloaded from the GEO database. The two profiles comprised 10 drug-resistant cell lines and 14 paired parental cell lines. GSE6914 carried one dataset, Calu3, consisting of four gemcitabine-resistant NSCLC cell lines, and four parental cell lines, while GSE77209 included two datasets, H1299T18 and H1355T16, comprising three taxane-platinresistant NSCLC cell lines (H1299 or H1355) and five of their parental cell lines. The characteristics of the three individual datasets are shown in Table 1. We first identified the DEGs from the above three datasets and listed the top 10 DEGs (Supplementary Tables S1-S3). We performed a two-dimensional hierarchical clustering analysis of DEGs for each dataset. The results revealed the significant differences in the drug-resistant cells as compared to their parental cells, as observed from the clustering of the DEGs. These observations indicate the high quality of these datasets and the reliability of the subsequent analysis ( Figure 2).

| GO analysis of DEGs
To investigate the functional role of the DEGs from the three independent datasets, we performed GO enrichment analysis. We observed 572 significantly enriched GO terms in Calu3 dataset, including 174 upregulated and 378 downregulated terms. In addition, we observed 307 GO terms in H1299T18 dataset, including 107 upregulated and 200 downregulated terms, and 343 GO terms in H1355T16 dataset, including 181 upregulated and 162 downregulated terms. The top 15 GO terms are shown in Figure 3, and the detailed information of the top 10 significantly upregulated or downregulated GO terms is listed in Supplementary Tables S4-S6. The results showed the alterations in the small molecule metabolic process in drug-resistant tumor cells.  Figure 4, and the detailed information of the top 10 significantly upregulated or downregulated KEGG pathways is listed in Supplementary Tables S7-S9. A great difference was found in the enriched upregulated and downregulated pathways among the three datasets. Moreover, we found that the treatment of H1355 and H1299 with the same drugs taxane and platin resulted in variations in the alteration patterns, indicating that the cell lineage may affect the activation of signaling pathways. We also found that some signaling pathways were aberrantly altered in the three datasets (eg, metabolic pathways), indicating that these pathways were central for acquiring drug resistance.

| Overlapping analysis of KEGG pathways
To further identify the KEGG pathways that were common in drug resistance, we performed an overlapping analysis of KEGG pathways that were enriched in all three datasets ( Figure 5). We reported 17 overlapping KEGG pathways, including 10 upregulated and seven downregulated pathways (Figure 6), the detailed information is listed in Table 2. The results showed that metabolic pathways, upregulated ubiquitin-mediated proteolysis, and MAPK signaling pathways were the major pathways mediating drug resistance, as evident from the most common significant alterations observed.

| Genes and common signaling pathways altered in cisplatin drug-resistant tumor cells
We next used cisplatin-resistant A549 tumor cells to validate the findings observed in our study. Some well-known drug resistance-related genes were detected with quantitative PCR, including multidrug resistance 1 (MDR-1), ATPbinding cassette subfamily C member 1 (ABCC1), ABCC2, ABCC3, and ABCC4. The results revealed that the expression of these genes was upregulated in cisplatin-resistant tumor cells, suggesting that these cell lines are resistant to cisplatin ( Figure 7A). We detected the mRNA expression levels of the genes that contributed to drug resistance identified in our study. As a result, we found that most of these genes were  Figure 7B). To further validate these results, we knockdown the expression of NAMPT with RNA interference, the efficiency of knockdown is shown in Figure 7C. The apoptosis of tumor cells was analyzed with flow cytometry, and the gating strategy is shown in Supplementary Figure S1. The knockdown of NAMPT expression dramatically increased the apoptosis of cisplatin-resistant tumor cells following cisplatin treatment ( Figure 7D), while minimum effects were observed on A549 cells. Western blot analysis showed that the expression of phosphorylated p44/42 MAPK(Thr202/Tyr204) was upregulated in cisplatin-resistant tumor cells ( Figure 7E) and unprocessed picture is shown in Supplementary Figure S2. These results suggest that the MAPK signaling pathways may be the common link in drug-resistant tumor cells.

| DISCUSSION
Somatic alterations in signaling pathways are common at varying frequencies and combinations in tumor cells and seem crucial in the development of resistance to various drugs. Therefore, the identification of the commonly altered signaling pathways in drug-resistant tumor cells is essential for the development of effective therapeutic strategies.
In this study, we compared the gene expression profiles of 24 samples comprising gemcitabine-resistant and taxane-platin-resistant NSCLC cell lines and their parental cell lines. We integrated three microarray datasets and identified the common signaling pathways associated with drug resistance. DEGs were identified for each dataset, and GO and KEGG enrichment pathway analysis for DEGs were performed to explore the molecular mechanisms underlying drug resistance development for each dataset. The functional enrichment analysis of GO terms and KEGG pathways showed striking differences between three drugresistant cell lines, indicating that the selective activation of signaling pathways is crucial for mediating drug resistance in tumor cells.
Drug resistance is a major obstacle observed during chemotherapy treatment, and different pathways are activated in the tumor cells in response to different drug treatments. Therefore, the identification of the common signaling pathways that are important to mediate drug resistance in NSCLC is desirable to eliminate drug resistance. We performed an overlapping analysis of three KEGG pathways for each dataset and found most significant alterations in metabolic pathways. Metabolic reprogramming is a hallmark of cancer development. Many studies have confirmed increased aerobic glycolysis, fatty acid synthesis, and glutamine metabolism to be associated with therapeutic resistance in cancer. 33 In breast cancer, fatty acid synthase (FASN) induces docetaxel/ trastuzumab/adriamycin resistance and lactate dehydrogenase A (LDHA) contributes to paclitaxel/trastuzumab resistance. 34,35 Aberrant metabolism has been thought to induce drug resistance in cancer cells; thus, the strategies targeting metabolism, for instance, glucose transporters (GLUTs), hexokinase (HK), pyruvate kinase M2 (PKM2), LDHA, pyruvate dehydrogenase kinase (PDK), and glutaminase (GLS), were shown to achieve a remarkable progress in reducing drug resistance in experimental and clinical studies. [36][37][38][39][40][41][42][43] Tavassoly et al 44 reported different metabolic pathways in allopurinolsensitive and -resistant cell lines. Fatty acid catabolic process and triglyceride process were enriched in the resistant cells, while the pathways and processes related to oxidative stress were likely to be dominant in sensitive cells.
In this study, we found that NAMPT upregulation was the most common. NAMPT is a rate-limiting enzyme in the salvage pathway for the biosynthesis of nicotinamide adenine dinucleotide (NAD + ) from nicotinamide. NAMPT has been reported to activate MAPK signaling. As per our results, MAPK signaling was synergistically upregulated with NAMPT expression. 45 In addition, NAMPT has also been associated with chronic inflammation in pancreatic cancer and is thought to contribute to drug resistance. 46 The knockdown of NAMPT expression significantly increased the apoptosis of cisplatin-resistant A549 tumor cells following cisplatin treatment. In contrast, the gluconeogenic enzyme F I G U R E 6 Common signaling pathways in three datasets. A, Upregulated common KEGG pathway in three datasets. B, Downregulated common KEGG pathways in three datasets fructose-1,6-bisphosphatase 1 (FBP1) is encoded by a wellknown tumor suppressor gene 47 and known to be downregulated in drug-resistant cell lines. These results indicate that the alternative metabolic reprogramming was essential in mediating resistance to drug killing. We similarly identified certain molecules such as alpha-aminoadipic semialdehyde synthase (AASS) and coenzyme Q 7 homolog ubiquinone (COQ7) whose functions were unknown in tumor cells.
Ubiquitin-mediated proteolysis pathway was the second most significantly upregulated pathway shared by the three datasets. Ubiquitin-mediated proteolysis is involved in most cellular processes, including cell cycle, DNA repair, transcription, apoptosis, angiogenesis, protein synthesis, polyamine biosynthesis, and antigen presentation. 48 The abnormal activation of ubiquitination has been linked to many diseases. Accumulating evidence has proved that the dysregulation in ubiquitin-mediated proteolysis may play a crucial role in tumorigenesis and progression. The mutation in Clb gene, which is required for EGFR internalization and lysosomal degradation, results in the inhibition of the ubiquitin-mediated degradation and has been linked to gastrointestinal tumor formation. 49 However, in our results, we identified that Cullin 4B (CUL4B), a member of . E, Representative western blot for total and phosphorylated P44/42 MAPK (Thr202/Tyr204) protein expression in A549 and CisR-A549 cells. β-actin was used as loading control. The number represents the protein size. Graphic display method refers to this articles. 62,63 The data in (A-C) were made log2 transformation, and analyzed by unpaired t tests. The data in (D) were made logit transformation and analyzed by unpaired t tests. The dots represent the mean value of the two technical repetitions, results are representative of three independent experiments E3 ubiquitin ligases, was overexpressed in drug-resistant cell lines. CUL4B promoted the degradation of P53 and inhibited the expression of phosphatase and tension homology (PTEN) through posttranscriptional modifications. 50 In addition, recent studies have reported that CUL4B increases the expression of human EGFR2 (HER2) in gastric cancer cells to promote tumor invasion. Zhang et al showed that ubiquitination regulates the stability of MDR1 gene product, P-glycoprotein, and affects the functions of this membrane transporter, resulting in multidrug resistance. This phenomenon was confirmed by Liu et al These authors found that the seven-in-absentia homologue 1 (Siah1), an E3 ubiquitin ligase that regulates ubiquitination, was downregulated in P-glycoprotein-mediated multiple drugresistant cancer cells. Siah1 exhibits the function of decreasing P-glycoprotein level, and low expression level of Siah1 was shown to induce multidrug resistance in cancer cells. 51,52 Therefore, the treatment of drug-resistant cells with β-elemene was shown to increase the expression of E3 ubiquitin ligases to enhance the efficacy of DOX treatment. 53 This evidence, as well as our study, indicated the important role of ubiquitination in tumor progression. MAPK signaling pathway is well-known to participate in mediating drug resistance and is crucial for the damageinduced DNA replication. Therefore, aberrant activation of MAPK signaling may prevent the death of cancer cells from drug treatment. 54 Numerous studies have confirmed that MAPK signaling is involved in drug resistance induced by various chemotherapy drugs such as gemcitabine, platinum, 5-fluorouracil, and tamoxifen. [55][56][57][58] Ercan et al 59 revealed the amplification of MAPK1 gene in EGFR kinase inhibitorresistant NSCLC cell lines. Inhibition of MAPK signaling using MEK inhibitor, CI-1040, resulted in the restoration of sensitivity to WZ4002. Furthermore, MAPK feedback activation was another mechanism underlying drug resistance observed following EGFR inhibitor treatment. Dysregulation of MAPK signaling pathway was shown to correlate with the progression of pancreatic ductal adenocarcinoma. The blockade of MAPK signaling using FBP1-derived small peptide could overcome the gemcitabine-induced ERK activation and increase the efficacy of anticancer treatment. 60,61 In our study, upregulated MAPK signaling level was reported with a concomitant overexpression of activating transcription factor 2 (ATF2), an isoform of P38 MAPK family. Furthermore, dual-specificity phosphatase 10 (DUSP10, MKP5) and downregulated FBP1 expression in drug-resistant NSCLC tumor cell lines were observed among the three datasets. Fang et al also revealed the upregulated expression of MAPK signaling pathway in cisplatin-resistant A549 cells, thereby supporting our analysis and validation (western blot analysis). These authors also found that the phosphoinositide 3-kinase (PI3K)/protein kinase B (AKT) signaling pathway was upregulated in cisplatin-resistant A549 cells. However, we failed to observe any consistent activation of PI3K/AKT signaling in the three datasets, as PI3K/AKT signaling was upregulated in gemcitabine-resistant Calu3 cell line but downregulated in taxane-platin-resistant NSCLC cell line (H1299 or H1355). 30 Janus kinase/ signal transducer and activator of transcription (JAK-STAT) signaling pathway have been reported to be upregulated in allopurinol-resistant cell lines, but we failed to report similar observation in either taxane-platin-resistant cells or gemcitabine-resistant cells. These differences may be related to the different killing mechanisms of the chemotherapy drugs.
To our knowledge, this is the first study to focus on the common signaling pathways involved in mediating drug resistance in NSCLC. Although many pathways identified in our study have already been known, our results prove that these pathways play an important role in mediating resistance to different drugs. In addition, we found that the genes in the same pathways may show different levels of alterations in response to different drugs, and revealed the functions of several genes involved in drug resistance. However, our study has a few limitations. For instance, only two types of drug-resistant cell lines (gemcitabine and taxane-platin) were used, only three cell lines were used, and the sample size was small. In addition, the backgrounds of tumor cell lines used were different, and our results lack the validation of mouse experiments and tumor specimen analysis. Further investigations in suitable mouse models and tumor samples from patients with NSCLC that are confirmed as drug resistant are warranted to validate the common pathways identified in this study.
In summary, our study provides a deep understanding of the signaling pathways that contribute to drug resistance in NSCLC and highlight some key molecules involved in these altered pathways. Notably, the most common signaling pathways identified in our study may facilitate the development of therapeutic interventions for drug-resistant NSCLC.

ACKNOWLEDGMENT
This work was supported by National Natural Sciences Foundation of China (No. 81771781).