Gene expression profiling and functional analysis reveals that p53 pathway-related gene expression is highly activated in cancer cells treated by cold atmospheric plasma-activated medium

Materials and Methods

Plasma device and PAM generation

The CAP device was a typical needle-ring plasma jet (Liu et al., 2016). The plasma jet is made of a polytetrafluoroethylene (PTFE) tube, a tungsten steel needle and a stainless steel ring. The PTFE tube has an inner diameter of 3 mm and the tungsten steel needle is coaxially placed in the center of the tube. The inner diameter of the stainless steel ring is 2 mm and is placed at the end of the PTFE tube. The distance between the tip of the needle and the ring electrode is 1 mm. The needle serves as high voltage (HV) electrode and the ring serves as the ground electrode. The HV electrode is connected to a DC power (purchased from Zhaofu Corp, Xian, China) via a 500 kΩ ballast resistor. When the air flows through the PTFE tube, a cold plasma plume is generated in the surrounding air. The plasma jet was used to vertically irradiate the media in a 6 cm plate. A total of 5 mL media were irradiated with the glow region of the plasma plume for 5 min (6 mm distance from the jet tip). Then, the media were balanced in cell incubator and used to stimulate SCC15 cells. The plume core (3 mm distance from the jet tip) temperature was detected under different conditions by Liu et al. (2016).

Cells and PAM treatment

Human oral cancer cell line SCC15 cells (Cat.3111C0001CCC000524) were purchased from National Infrastructure of Cell Line Resource (http://www.cellresource.cn) and kept in our laboratory. SCC15 cells were generally cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Life technologies, NY, USA) supplemented with 10% FBS, 50 µM β-ME, 100 U/ml penicillin, and 100 μg/ml streptomycin under the standard cell culture conditions (37 °C, 5% CO₂ environment). 1.5 × 10⁶ SCC15 cells were seeded in 6 cm plates (Corning, NY, USA) the day before treatment. After 24 h culture, cells were washed with PBS 3 times. PAM was generated as mentioned above and added to the culture at a final volume of 5 mL. Then cells were incubated for 1 h at 37 °C. The medium was used as a control. The cells were washed and harvested with 1 mL Trizol reagent (Life technologies, NY, USA) by using a cell scraper. Then, the samples were quick-frozen in liquid nitrogen for RNA isolation.

Cell viability assay

Cell Viability was measured with the Cell Counting Kit-8 (CCK-8; Transgen Biotech, Beijing, China) under the manufacturer’s instruction. In brief, 1.0 × 10⁴ SCC15 cells were seeded in 96-well plates the day before treatment. After 24 h culture, PAM was generated and added to each well at a final volume of 50 µL for the indicated times. Here we did a time course stimulation (0, 15, 30, 60, 120 min). To facilitate detection and promote consistency, we added the PAM in reverse order. After treatment, cells were washed and incubated with the CCK-8 reagent at the same time (90 µL fresh media plus 10 µL CCK-8). Then cells were incubated for 1 h and the absorbance of OD 450 nm was measured by using a Microplate Reader (SynergyHTX, BioTek, Winooski, VT, USA). The relative viability was calculated according to the manufacturer’s instructions.

RNA isolation, cDNA library construction and RNA sequencing

Three independent experiments were performed for each group to obtain biological replicates. For sample preparation, 3 µg RNA per sample was used as input material. DNA contaminations were removed from the samples with DNase I. The purity, concentration and integrity of the RNA were assessed using a Nanodrop Spectrophotometer (Thermo Scientific, MA, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, CA, USA). The sample libraries were constructed with Illumina TruSeq RNA Sample Preparation Kit (Illumina, CA, USA) according to the manufacturer’s instructions. The quantity and quality of the libraries were assessed using the Qubit Fluorometer (Life Technologies, CA, USA) and the Agilent Bioanalyzer 2100 System. The libraries were then sequenced using the Illumina HiSeq platform.

Differential gene expression analysis

To obtain high-quality clean reads for subsequent analysis, adaptor reads, low-quality reads and reads with more than 10% unknown bases (poly-N) were removed from raw data. The Q30 and GC content of the clean data were also calculated. The clean reads were mapped to the human reference genome (GRCh37) using TopHat2 (Kim et al., 2013). The expression levels were calculated by using the fragments per kilobase per million reads method (FPKM) (Mortazavi et al., 2008). Differential gene expressions were analyzed by using the DEseq R package (Anders & Huber, 2010). The absolute value of Fold Change ≥2 and False Discovery Rate (FDR) ≤0.01 was adopted as the standard for determining the significance of differential gene expression.

Functional analysis of DEGs

To identify functional categories of differentially expressed genes, Gene Oncology (GO) enrichment analysis was performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) (Huang da, Sherman & Lempicki, 2009b; Huang da, Sherman & Lempicki, 2009a). Following the instructions of the DAVID manual, DEGs were uploaded and the function charts were generated. The groups with a P-value <0.05 and gene counts more than two were examined. KEGG pathway enrichment analysis of DEGs was performed using KOBAS software (Xie et al., 2011). The rich factor is calculated as the ratio of the numbers of DEGs enriched in this pathway, to the numbers of all genes annotated in the same pathway. The Q-value is the corrected P-value with threshold <0.05. Gene Set Enrichment Analysis (GSEA) was performed using the Java GSEA implementation in our study (http://software.broadinstitute.org/gsea/index.jsp) (Subramanian et al., 2005). The gene lists of KEGG or hallmark gene signature were adopted from The Molecular Signatures Database (MSigDB). As a metric for ranking genes in GSEA, the absolute signal to noise value of gene expression was used, and other parameters were set to default values.

Quantitative RT-PCR

Cells were processed and harvested under the same conditions as the sequencing sample. Total RNA was isolated by using Trizol reagent and was reverse-transcribed with the PrimeScript RT Kit (TAKARA Biotech, Dalian, China). The qRT-PCR assays were conducted on the Bio-Rad CFX96 System using the TransStart Tip reagent (TransGen Biotech). Primers used in the present study were designed by using Primer premier software (Table S1). Relative mRNA levels were calculated by the comparative cycle threshold method and GAPDH was used as the internal control for each sample.

CAP device and PAM treatment

The CAP device is a typical needle-ring plasma jet designed by Liu et al. (2016). Ambient air was used as the working gas, which has been used by several groups in anti-cancer research (Kim et al., 2010a; Kim et al., 2010b; Arndt et al., 2013). A schematic diagram of the system is presented in Fig. 1A. Cold plasma was generated when air flowed through the PTFE tube. The air was supplied by a pump and the flow rate controlled by gas flowmeter at a rate of 3 L/min, plus the input voltage of DC power was −5 KV. PAM was generated by using a general strategy. The plasma jet was used to vertically irradiate the media in a 6 cm plate. After irradiation, the media were used to stimulate cancer cells.

To characterize the kinetics of PAM-treated cancer cells, relative viabilities were detected by the CCK-8 method. For 15 min treatment, PAM was able to stimulate a decrease of cell viability of target SCC15 cells (85.22 ± 2.8 %) (Fig. 1B). Relative viability was kinetically augmented with increasing incubation time of PAM (30 min, 85.0 ± 4.9%; 60 min, 71.8 ± 6.6%; 120 min, 51.3 ± 5.9%). The results suggested that PAM showed a noticeable anti-cancer ability on SCC15 cancer cells. These data were representatives of at least three separate experiments. To elucidate the early transcriptional changes, sequencing samples were prepared in 6 cm plates for 60 min incubation with PAM.

RNA sequencing and read mapping

We established 6 libraries with the following designations. PAM1, PAM2 and PAM3 were from PAM-treated SCC15 cells while Ctr1, Ctr2 and Ctr3 were mock-treated with three biological replicates respectively. The cDNA libraries were prepared and sequenced by Illumina HiSeq platform, which generated a total of 175,673,694 raw reads (Table S2A). After removing the adaptors and low-quality reads, we obtained 172,801,869 clean reads, with a high quality of Q30 ≥ 88.98%. We then mapped the trimmed clean reads onto the human genome and 75.21% to 79.59% of the clean reads were mapped uniquely to the genome, whereas a small percentage of them were mapped multiple times onto the genome (Table S2B). The whole subsequent analysis was based on the uniquely mapped reads.

Differentially expressed genes (DEGs) in response to PAM

The levels of gene expression were expressed as FPKM values and displayed a comparable distribution among 6 parallel libraries. A total of 15,446 genes were annotated in this study and are available in Table S3. The biological replicates of both control and PAM-treated samples exhibited high correlation (Pearson’s Correlation Coefficient R2 >0.81) in FPKM values (Fig. S1). DEseq was employed to screen for DEGs and a total of 934 annotated genes were identified to be differentially expressed when considering exclusively a stringent threshold of FDR ≤ 0.01 and fold-change ≥ 2 (Table S4). To obtain a global view of these DEGs, hierarchical clustering was constructed with normalized FPKM values to characterize changes across six samples (Fig. 2). The most up-regulated gene was Glutathione-specific gamma glutamylcyclotransferase 1 (CHAC1, 16-fold), a pro-apoptotic component of the unfolded protein response, suggesting the involvement of apoptotic pathway after PAM stimulation.

Gene ontology (GO) classifications

GO is a universally acknowledged gene functional enrichment database and is generally applied to search for enriched GO terms in DEGs (Ashburner et al., 2000). We performed GO analysis using web-based DAVID tool (Huang da, Sherman & Lempicki, 2009b; Huang da, Sherman & Lempicki, 2009a). 845 out of 934 profiled DEGs were assigned to 1,211 GO terms, including 1,101 biological processes, 16 cellular components and 94 molecular function terms. Typical enriched GO terms are shown in Fig. 3. The GO terms of molecular function category were concentrated in “binding” (704 DEGs, 75.4%) and “protein binding” (549 DEGs, 58.8%) The highest percentages of GO terms under cellular component class were “cell” (785 DEGs, 84.0%) and “cell part” (783 DEGs, 83.8%). We noted that GO terms were mainly categorized into “biological processes”, with wide distributions and extensive assignments. The most prevalent “biological processes” assignment was “cellular”, including 759 DEGs and accounting for 81.3% of all DEGs. Besides, some important assignments, such as “biological regulation” (582 DEGs, 62.3%), “metabolic” (573 DEGs, 61.3%) and “signaling” (392 DEGs, 42.0%), were highly enriched, suggesting that the biological processes of cancer cells were widely changed after PAM treatment. The data of GO classifications are available in Table S5.

“P53 pathway” was highly enriched by KEGG mapping

KEGG database is a collection of various pathways, representing the molecular interactions and reaction networks (Kanehisa et al., 2004). To identify signaling pathways involved in PMA-treated cells, we had mapped the KEGG database and found that identified DEGs were significantly enriched in 40 KEGG pathways (Q value ≤ 0.01, Table S6). The top 20 enriched pathways are shown in Fig. 4. DEGs were highly clustered in several signaling pathways, such as “p53 signaling”, “Hippo signaling”, “TNF signaling”, “AGE-RAGE signaling” and “FoxO signaling”, suggesting that PAM may perform its function through these pathways.

The “p53 pathway” was identified as the most significantly enriched pathway in KEGG analysis (Rich Factor = 0.23 and Q value < 0.001). A total of 16 DEGs were enriched in this pathway (Fig. 5A). Moreover, 12 genes (GADD45A, SERPINE1, SERPINB5, PPM1D, IGFBP3, SESN2, CCNE2, GADD45B, THBS1, GADD45G, CDKN1A, PMAIP1) were up-regulated and 4 genes (MDM-4, CCNB2, TP73, ADGRB1) were down-regulated (Fig. 5B). To further confirm these observations, eight genes were selected for qRT-PCR validation (Fig. 5C). The results of the qRT-PCR analysis were also consistent with our RNA-seq data and validate its reliability. In addition, we detected CCNB2 and Apaf1 expression by Western blot (Fig. S2). CCNB2 is a member of the cyclin family and is essential for cell cycle regulation. Apaf1 is the key component of apoptosome and acts in both p53 and apoptosis pathways. The mRNA expression levels of CCNB2 and Apaf1 were regulated in SCC15 cells after PAM stimulation (Figs. 5 and 6). We incubated SCC15 cells with PAM for different times and validated that the protein levels of CCNB2 and Apaf1 were changed, which were in agreement with our previous results.

Gene set enrichment analysis (GSEA)

GSEA is a useful computational algorithm that determines whether a pre-defined set of genes (signature) shows differences between two groups (Mootha et al., 2003; Subramanian et al., 2005). To gain further insights at the gene-set level, GSEA was performed against KEGG or hallmark gene-set signature. Consistently, “p53 pathway” was found to be enriched in PAM phenotype (Fig. 6A). Importantly, instead of 16 DEGs identified by KEGG analysis, 19 more genes with fewer changes (Fold Change ≤ 2) were found in GSEA (Figs. 6B, 5A). In addition, apoptosis and hypoxia pathways were also identified as the principal involved pathways in this assay, although they were not enriched in KEGG analysis (Fig. S3 and Table S7). Our results not only establish the valuable application of GSEA in expression analysis but also provide novel insights into cancer cell signaling induced by PAM.