Differential Expression and Alternative Splicing of Transcripts Associated With Cisplatin-Induced Chemoresistance in Nasopharyngeal Carcinoma

Radiotherapy and adjuvant cisplatin (DDP) chemotherapy are standard administrations applied to treat nasopharyngeal carcinoma (NPC). However, the molecular changes and functions of DDP in NPC chemo-resistance remain poorly understood. In the present study, transcriptomic sequencing between 5-8F and 5-8F/DDP cells was performed to identify differential expression and alternative splicing (AS) characteristics in DDP-resistant NPC cells. Transcriptomic profiling identified 1,757 upregulated genes and 1,473 downregulated differentially expressed genes (DEGs). Bioinformatic analysis revealed that these DEGs were associated with or participated in important biological regulatory functions in NPC. Validation of 20 signiﬁcant DEGs using quantitative real-time reverse transcription PCR showed that the expression patterns of 17 mRNAs were in accordance with the sequencing data. Intron retention was identified as the major AS event in chemoresistant cells. Furthermore, the expression level of matrix metalloproteinase 1 (MMP1), which was one of the most upregulated mRNAs in the chemoresistant cell lines, was signiﬁcantly associated with the migration, invasion, and proliferation of NPC cells in vitro. Our study revealed that dysregulated genes and AS-mediated DDP chemoresistance might play important roles in NPC development and progression. Targeting aberrantly expressed genes might clarify the pathogenesis of NPC and contribute to developing new therapeutic strategies for NPC.


INTRODUCTION
Nasopharyngeal carcinoma (NPC) is one of the most highly invasive and metastatic head and neck cancers in southern China and Southeast Asia (Cao et al., 2011;Torre et al., 2015). The vague symptoms and aggressiveness of NPC mean that more than 70% of patients with NPC are diagnosed with locoregionally advanced disease and often have an unfavorable prognosis (Lai et al., 2011;Chua et al., 2016). According to the National Comprehensive Cancer Network (NCCN) guidelines, platinumbased concurrent chemoradiotherapy is the standard of care for patients with locoregionally advanced nasopharyngeal carcinoma . However, many patients develop recurrence, distant metastasis, or both because of chemoresistance (Brabec and Kasparkova, 2005;Xie et al., 2010). Thus, the chemoresistance of NPC to platinum-based therapy not only affects treatment efficacy, but also impacts on the prognosis of patients with NPC.
Cisplatin (DDP), a DNA damaging agent, exhibits its cytotoxicity and apoptosis-inducing activities by forming DNA adducts or by targeting cancer signaling pathways (Chen et al., 2015;Amable, 2016;Wang et al., 2018). Combinations of DDP-based chemotherapies have been widely used to treat various cancers, including NPC (Agbale et al., 2016;Dugbartey et al., 2016;Sun et al., 2016). However, high-dose DDP is frequently associated with severe vomiting, ototoxicity, and hematotoxicity (Forastiere et al., 2003;Ang et al., 2014;Tang et al., 2018). In addition, clinical studies indicate that many patients acquire DDP resistance during cancer chemotherapy. Understanding the molecular mechanism underlying DDP chemoresistance is of crucial importance to develop novel therapeutic strategies to treat NPC.
In the present study, we established 5-8F and SUNE-1 DDP chemoresistance models to investigate the differential gene expression and alternative splicing (AS) events associated with NPC chemoresistance. First, we compared the chemoresistance of 5-8F DDP and 5-8F cell lines using transcriptomic sequencing and quantitative real-time reverse transcription PCR (qRT-PCR) to identify dysregulated mRNAs that participate in chemoresistance. Then, Gene Ontology (GO) and pathway analysis were performed to better understand the differentially expressed mRNAs. Next, we selected one of the most upregulated genes, MMP1 (matrix metalloproteinase 1), to validate migration, invasion, and proliferation functions in vitro. Subsequently, aberrant alternative splicing events were further explored, which had not been previously reported in NPC chemoresistance. Our research provides new insights into chemoresistance of nasopharyngeal carcinoma.

Cell Counting Kit 8 (CCK8) Assay
Cells (1 × 10 3 ) were seeded in each well of 96-well plates for overnight incubation. The cells were then treated with different concentrations of DDP for 72 h. A CCK8 assay was then performed to assess cell viability, as described previously  Cell Apoptosis Assays The cell death of 5-8F, SUNE1, 5-8F/DDP, and SUNE1/DDP cell lines induced by DDP (2.5 mg/ml) was analyzed using flow cytometry with Annexin-V/propidium iodide (PI) assays according to the manufacturer's instructions (BD Biosciences, Bedford, MA, USA).

RNA Extraction and Quality Control
Total RNA was extracted from the NPC cell lines using the TRIzol reagent (Invitrogen, Grand Island, NY, USA) according to the manufacturer's instructions. RNA integrity was assessed by standard denaturing agarose gel electrophoresis. The quality and amount of RNA was assessed using a NanoDrop ND-1000 spectrophotometer (Thermo Scientific, Rockford, IL, USA) and the isolated RNAs were stored at −80°C before transcriptomic sequencing analysis.

Construction cDNA Library
RNA (2 mg) was utilized in the RNA sample preparation, strictly according to the manufacturer's protocol. First, ribosomal RNA was removed using a RiboZero Magnetic Gold Kit (Epicentre, Madison, WI, USA), and residual RNAs were cleaned using ethanol precipitation. Sequencing libraries were generated using the rRNA-depleted RNA with KAPA Stranded RNA-Seq Library Prep Kit (Illumina, San Diego, CA, USA). RNA integrity was evaluated by using an RNA Nano 6000 Assay Kit from the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). The libraries were sequenced on an Illumina Hiseq 4000 platform.

Identification of Differently Expressed Genes
The analysis of differences in mRNA expression between the two groups was performed using the DEGseq (2010) R package. The P-value was adjusted using the q-value. The threshold for significant differential expression was set as a q-value < 0.05 and |log2(fold change)| > 1.5 by default.
Quantitative Real-Time Reverse Transcription PCR qRT-PCR was performed using SYBR Green PCR Master Mix (Applied Biosystems,Foster City, CA, USA) on a CFX96 Touch Sequence Detection System (Bio-Rad, Hercules, CA, USA). All samples were normalized to internal controls and fold changes were calculated using relative quantification (2 −DDCt ). All experiments were performed in triplicate. The primer sequences are shown in Supplementary Table S1.

GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment Analysis
The differentially expressed mRNAs were selected and subjected to GO and KEGG pathway analysis. For GO analysis (http:// geneontology.org/), the corresponding genes were annotated and classified according to biological process (BP), cellular component (CC), and molecular function (MF). For KEGG analysis (http://www.genome.jp/kegg/), the differentially pathways were ranked by their enrichment scores.

Western Blotting Analysis
Cells were lysed on ice in Radioimmunoprecipitation assay (RIPA) buffer. The protein concentration was determined using the Bradford method. Protein lysates (20 mg) were separated by sodium dodecyl sulfate polyacrylamide gel electrophoresis (SDS-PAGE) and electrophoretically transferred to a polyvinylidene fluoride membrane (Millipore, Billerica, MA, USA). The membranes were then blocked with 5% skim milk and incubated with antibodies against MMP1 (Cell Signaling Technology,Danvers, MA, USA, 1:10,000) and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) (Proteintech, Rosemont, IL, USA; 1:5,000) overnight at 4°C. Species-matched secondary antibodies were then added and incubated at room temperature for 1 h. The immunoreactive proteins were detected using BeyoECLPlus (Beyotime, Shanghai, China).

Transient Transfection and Stable Cell Line Establishment
Plasmids pEnter-MMP1-Flag and pEnter-vector were obtained from Vigene Bioscience (Rockville, MD, USA). Short interfering RNA (siRNA) oligonucleotides targeting MMP1 were purchased from GenePharma (Shanghai China); the siRNA sequences are shown in Supplementary Table S1. Plasmid (2 mg) and siRNA oligonucleotide (100 nmol/L) transfections were carried out using Lipofectamine 3000 (Invitrogen) according to the manufacturer's instructions. The cells were used for further studies at 48 h after transfection.

Wound Healing Assay
HONE1 and SUNE1 cells grown to near confluence in 6-well plates were incubated in serum-free medium for 24 h of starvation. Linear wounds were created in the cell monolayers using a P-200 pipette tip, followed by an additional 48 h of starvation. Images were captured and documented under a microscope at × 100 magnification at 0 and 24 h.

Migration and Invasion Assays
Cells (5 × 10 4 or 1 × 10 5 ) suspended in 200 ml of serum-free medium were plated into Transwell chambers (8-mm pores, Corning, Corning, NY, USA) precoated with (migration assay) or without (invasion assay) Matrigel (BD Biosciences) and cultured at 37°C for 12 or 24 h. The cells were fixed with methyl alcohol, stained with crystal violet, and counted under a microscope in 10 random fields of view per well.

Statistical Analysis
SPSS statistical software version 19.0 was used for data analysis (IBM Corp., Armonk, NY, USA). Chi-square and ANOVA tests were used to analyze the differences among the two groups. Statistical significance was set at a P < 0 .05. The RNA sequencing data has been deposited in the Gene Expression Omnibus (accession number: GSE135083).

Establishment of DDP-Resistant NPC Cell Lines
To establish the DDP-resistant cells, SUNE1 and 5-8F cell lines were exposed to increasing concentrations of DDP for more than 6 months. Compared with their parental cells, chemoresistance to DDP was observed in 5-8F/DDP and SUNE1/DDP cells. As shown in Figures 1A, B, 2.5 mg/ml DDP led to about 50% cell growth inhibition in both SUNE1 and 5-8F cells, respectively. However, SUNE1/DDP and 5-8F/DDP cells produced resistance to the growth inhibitory properties of 2.5 mg/ml DDP. Thus, the resistant NPC cells were continuously maintained in RPMI 1600 medium containing 2.5 mg/ml DDP. To further validate the DDP-resistant NPC cell lines, flow cytometry assay was used to detect the apoptosis ability of 5-8F/DDP and SUNE1/DDP cells. As shown in Figures 1C, D, under 2.5 mg/ml DDP treatment, 5-8F/DDP and SUNE1/DDP cells showed a lower apoptosis rate than 5-8F and SUNE1 cells.

RNA Sequencing of Chemoresistance Cells
5-8F and 5-8F/DDP cells were collected to perform a standard RNA sequencing analysis. The outcome of sequencing demonstrated that 11,798 mRNAs were detected in total. The mRNA expression patterns of the samples are presented as heat maps (Figures 2A, B). To clarify the expression signatures of the dysregulated mRNAs, we analyzed upregulated or downregulated mRNAs identified in 5-8F/DDP cell lines according to their classification and chromosome distribution ( Figure 2C). Compared with the 5-8F cells, 3,230 mRNAs were differentially expressed in 5-8F/DDP cells |log2(fold change)| > 1.5), among which 1,757 were upregulated and 1,473 were downregulated ( Supplementary Tables 2 and 3). To verify the transcriptomic data, we selected the 20 most significantly dysregulated mRNAs, including 10 upregulated mRNAs and 10 downregulated mRNAs, and then validated their expression levels using qRT-PCR. The results showed that the expression patterns of 17 mRNAs were consistent with the sequencing data ( Figures 2D, E), which demonstrated the reliability of the RNAseq data.

GO and KEGG Pathway Analysis
To explore the potential function of the differentially expressed mRNAs in chemoresistance, GO analysis was performed to describe biological process (BP), cellular component (CC) and molecular function (MF; Supplementary Tables 4-6). The GO terms were determined by calculating the Enrichment Score (P < 0.05). The aberrantly expressed genes were mainly enriched for GO terms related to regulation of anatomical structure morphogenesis, regulation of signal transduction, and response to organic substance involved (BP); endomembrane system, cytoplasm region, and vesicle (CC); and protein binding, peptidase binding, and amide binding activity (MF). The top ten highest and most significant GO terms are shown in Figures 3A-C.
Pathway analysis based on the KEGG database was performed, which identified 20 pathways with significant differences (P < 0.05) in gene expression (Supplementary Table 7). The pathway terms of top ten highest Enrichment Scores are shown in Figures 4A, B. The pathway analysis results suggested that the upregulated mRNAs were part of several signaling pathways, including pathways in cancer (hsa05200), fluid shear stress and atherosclerosis (hsa05418), focal adhesion (hsa04510), and the apoptosis pathway (hsa04210). However, the downregulated mRNAs participated in metabolic pathways (hsa01100), citrate cycle (TCA cycle) (hsa00020), propanoate metabolism (hsa00640), and steroid biosynthesis pathways (hsa00100), which may be associated with NPC chemoresistance.

Alternative Splicing Analysis
To clarify the potential AS in 5-8F/DDP cells, A3SS, A5SS, MXE, RI, and SE events were analyzed. As shown in Figure 5A, 6310 AS events (2,524 upregulated and 3,786 downregulated; Supplementary Table 8) were identified as differentially expressed in 5-8F/DDP cells compared with 5-8F cells. The percentage of A3SS, A5SS, MXE, RI, and SE AS events were 9.6%, 7.2%, 8.2%, 16.6%, and 58.4%. Thus, SE was one of the most frequent AS events. The results showed that AS events were extensive and complicated in chemoresistance. To determine the signatures of the dysregulated AS events, the top ten upregulated or downregulated AS events in A3SS, A5SS, MXE, RI, and SE were validated using rMATS ( Figures 5B-N). The results indicated that AS is an important molecular event and might play a vital role in chemoresistance in NPC.

Silencing MMP1 Promotes NPC Cell Migration and Invasion In Vitro
To evaluate whether aberrant expression of MMP1, which was one of the most upregulated mRNAs in the chemoresistant cell lines, affects the metastasis ability of NPC cells, HONE1, and SUNE1 cells were transiently transfected with MMP1 plasmid and siRNAs targeting MMP1. As shown in Figure 6A, western blotting validated that MMP1 protein level was obviously elevated, meanwhile, increased vimentin expression and decreased E-cadherin expression after overexpression of MMP1 in NPC cells. The wound healing and Transwell migration and invasion assays showed that the cell migration, and invasion abilities of HONE1 and SUNE1 cells stably overexpressing MMP1 were remarkably reduced compared with that of cells transfected with the vector plasmid ( Figures  6B, C). Conversely, silencing MMP1 increased the migratory and  invasive abilities of NPC cells (Figures 6D-I). These results suggest that MMP1 inhibits NPC cell migration and invasion in vitro.

DISCUSSION
NPC is a malignancy arising from the nasopharyngeal epithelium, and it is mainly endemic in Southeast Asia and China (Tsao et al., 2014). NPC is frequently diagnosed with locoregionally advanced disease because of its anatomic location and ambiguous symptoms. Radiotherapy combined with chemotherapy is the standard therapeutic approach to treat NPC. However, because many patients develop either intrinsic or acquired drug resistance to DDP, clinical treatment failure rates remain high. The remarkable difference in treatment benefits might be associated with genetic factors. In the present study, we compared the mRNA expression differences between 5-8F/DDP and 5-8F cells, and further explored , were consistent with the RNA sequencing data. All experiments were performed at least three times; data are mean ± SD. *P < 0.05, **P < 0.01 vs. control, Student's t-test.

Zhang et al. Transcriptomic Profiling in NPC Chemoresistance
Frontiers in Genetics | www.frontiersin.org February 2020 | Volume 11 | Article 52 the expression and function of a dysregulated gene and AS, which might regulate chemoresistance in NPC. Accumulating evidence demonstrates the close association between chemotherapy resistance and dysregulated gene expression. Several studies have confirmed that the aberrant expression of certain mRNAs is involved in NPC che more sistance. Peng et al. (2019) showed that hypermethylation of ARNTL (encoding aryl hydrocarbon receptor nuclear translocator like) promotes NPC tumorigenesis and inhibits cisplatin sensitivity by activating CDK5 (cyclin   (Zhang B. et al., 2019). However, few studies have systematically explored the relationships between dysregulated expression of mRNAs and chemoresistance in stable DDP-resistant NPC cells.
results (Liu et al., 2013). Our functional assays validated that MMP1 could increase NPC proliferation, migration, and invasion, which indicated that MMPs play an important role in NPC chemoresistance.
GO analysis predicated that the aberrant mRNAs were associated with BP, CC, and MF in NPC. The GO terms such as regulation of signaling and regulation of cell migration and motility in BP predicted that the associated gene products led to the chemoresistance of NPC. The results of pathway analysis showed that the dysregulated mRNAs were associated with 20 signaling pathways in NPC. The correlation between these pathway and tumor progression, including NPC, has been proven in previous studies. Especially, focal adhesion (Liu et al., 2018;Yuan et al., 2019;Yu et al., 2019), apoptosis (Li M. et al., 2019;Liang et al., 2019;Wang et al., 2019;Xie et al., 2019), and the tumor necrosis factor (TNF) signaling pathway (Lu et al., 2011;Ou et al., 2015;Huang et al., 2017;Deng et al., 2018) have been reported to be closely related to NPC development.
AS changes are frequently observed in cancer, and AS is starting to be recognized as an important signature in tumor progression and therapy (Climente-Gonzalez et al., 2017). However, the impact and relevance of AS on tumorigenesis All experiments were performed at least three times; data are mean ± SD. *P < 0.05, **P < 0.01 vs. control, Student's t-test. and chemoresistance remain mostly unknown. In the present study, we carried out a systematic analysis to characterize the AS events and found that AS (A3SS, A5SS, MXE, RI, SE) were frequently dysregulated in DDP chemoresistant cells, especially RI (13.7%). Recently, several studies confirmed that AS is also involved in tumor malignant features. Downregulation of QK1 (QKI, KH domain containing RNA binding) led to the alternative splicing change in NUMB (NUMB endocytic adaptor protein), which promoted cell proliferation (Zong et al., 2014). An SE event in MST1R (macrophage stimulating 1 receptor) has been related to the acquisition of cell motility during cancer cell invasion (Ghigna et al., 2005). The smallmolecule modulators of pre-mRNA splicing are capable of restoring the original BRAF (B-raf proto-oncogene, serine/ threonine kinase) splicing and rescue growth of therapyresistant cells (Salton et al., 2015). The modulation of these events can recapitulate the tumor phenotype or revert to a normal phenotype (Ghigna et al., 2010;Bechara et al., 2013). Thus, determining alterations in AS is essential to understand the functional chemoresistance in cancer and targeting AS might provide a new way to alleviate chemoresistance.
In conclusion, dysregulated gene-and AS-mediated DDP chemoresistance might play important roles in NPC development and progression. These findings indicated that targeting aberrantly expressed genes and AS might provide novel insights into the pathogenesis of NPC and contribute to the development of new therapeutic strategies to treat NPC.

DATA AVAILABILITY STATEMENT
All datasets generated and analyzed for this study have been deposited in the Gene Expression Omnibus (accession number: GSE135083).

AUTHOR CONTRIBUTIONS
JZha, HJ, and TX designed the research. TX, JZhe, YT, RL, BW, JL, and AX acquired and analyzed the data. JZha, HJ, and YY wrote the manuscript.