HeLa (immortalized human cervical cancer) cells were purchased from the Institute of Biochemistry and Cell Biology, Chinese Academy of Sciences, HepG2 (immortalized human liver cancer) cells were purchased from the American Type Culture Collection (ATCC; Manassas, VA, USA). The cell lines were identified by short terminal repeat (STR) genotyping, which revealed a correspondence of >80% of the markers tested (32). Cells were cultured in Dulbecco's modified Eagle's medium (DMEM) with 10% fetal bovine serum (FBS).

To knockdown CK18, different shRNAs against human CK18 (shRNA1: 5′-GATGACACCAATATCACACGA-3′; shRNA2: 5′-CTTCATGAAGAAGAACCACGA-3′; shRNA3: 5′-CCTGCTGAACATCAAGGTCAA-3′) were designed, and scrambled shRNA (Scr-shRNA) that targeted a non-specific sequence (5′-ACTGGACCAGGCAGCAGCGTCAGAAGACT-3′) was used as the control. These shRNAs were transfected into HeLa cells using Lipofectamine 2000 transfection reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA), according to the manufacturer's protocol. Transfected cells were harvested after 48 h for analysis.

Total RNA was isolated from cells using TRIzol reagent (Ambion; Thermo Fisher Scientific, Inc) and cDNA synthesis was carried out using standard procedures. qPCR was performed on the Bio-Rad S1000 Thermal Cycler, using Bestar SYBR-Green RT-PCR Master Mix (DBI Bioscience, Shanghai, China). The PCR conditions are consisted of denaturing at 95°C for 10 min, 40 cycles of denaturing at 95°C for 15 sec, annealing and extension at 60°C for 1 min. The primers of CK18 used for quantitative real-time PCR (qPCR) were: Forward, AAAGGCCTACAAGCCCAGAT and reverse, CACTGTGGTGCTCTCCTCAA. Gene expression levels were calculated using the 2^−ΔΔCq method (33) and CT values were normalized using glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as an internal standard. The primers of GAPDH were: Forward, GGTCGGAGTCAACGGATTTG and reverse, GGAAGATGGTGATGGGATTTC.

An MTT assay was used to assess cell proliferation. Briefly, HeLa cells were cultured in 96-well plates and transfected with the vector using Lipofectamine 2000, according to the manufacturer's protocol. The cells were then incubated at 37°C for 48 h. Subsequently, MTT solution (5 mg/ml, 0.025 ml) was added to each well, and the cells were incubated for another 4 h. After centrifugation at 4,000 rpm for 15 min, the supernatant was removed from each well. The colored formazan crystals produced from MTT in each well were dissolved in DMSO (0.15 ml) and the optical density (OD) values were measured at 490 nm.

HeLa cells (5×10⁴) were seeded in 24-well culture plates. Once the cells reached 70% confluence, the cells were were transfected with the vector using Lipofectamine 2000, according to the manufacturer's protocol. The cells were then incubated at 37°C for 48 h and viable cells were harvested and washed twice with phosphate-buffered saline (PBS). Viable cells were double-stained with fluorescein isothiocyanate (FITC)-conjugated Annexin V and 7-amino actinomycin D (7-AAD) (4A Biotech Co., Ltd. Beijing, China).

The percentage of cell apoptosis was defined as the sum of the right lower quadrant and upper quadrant.

The RNA was treated with RQ1 DNase (Promega Corp., Madison, WI, USA) to remove DNA. The quality and quantity of the purified RNA were assessed by measuring the absorbance at 260 nm/280 nm (A260/A280) using SmartSpec Plus spectrophotometer (BioRad Laboratories, Inc., Hercules, CA, USA). The integrity of each RNA sample was further verified by 1.5% agarose gel electrophoresis.

For each sample, a total of 1 µg RNA was used for RNA-seq library preparation. PolyA mRNAs were purified and concentrated with oligo (dT)-conjugated magnetic beads (Invitrogen; Thermo Fisher Scientific, Inc.) before directional RNA-seq library preparation. The fragmented RNAs were end repaired and an adaptor sequence was ligated at the 5′end. Reverse transcription was then performed using an RT primer harboring a 3′adaptor sequence and randomized hexamer. The cDNAs were purified, amplified, and stored at −80°C until they were used for sequencing.

For high-throughput sequencing, the libraries were prepared following the manufacturer's instructions. An Illumina HiSeq 4000 sequencing system (ABLife, Inc., Wuhan, China) was used to collect data from 151-bp pair-end sequencing.

Raw reads were filtered to remove the adaptors, PolyN reads, and low-quality bases using FASTX-Toolkit (version 0.0.13). Short reads, less than 16 nt, were also discarded. Clean reads were then aligned to the GRCh38 genome using TopHat2 software (34), with four mismatches. Uniquely mapped reads were used to calculate the reads number and FPKM value (fragments per kilobase of transcript per million fragments mapped) for each gene.

EdgeR software (35) was used to measure FPKM values and to analyze the differential expression of genes using RNA-Seq data, in order to identify DEGs. To determine whether a gene was differentially expressed, the results were analyzed based on the fold change (fold change ≥2 or ≤0.5) and a false discovery rate (FDR<0.05).

We also analyzed the two replicates separately (named as simple pair). In detail, the DEGs were identified between shCK18-1st vs. Ctrl-1st and shCK18-2nd vs. Ctrl-2nd, respectively. We then obtained the overlapping genes from the DEGs of each replicate. The upregulated and downregulated DEGs shared by these two simple pairs were then overlapped with those identified by edgeR as biological replicates.

Gene Ontology (GO) and enriched KEGG pathway analyses were carried out using the KOBAS 2.0 server (36) to predict the gene function and calculate the frequency distribution of functional categories. The hypergeometric test and Benjamini-Hochberg FDR controlling procedure were used to define the enrichment of each pathway (corrected P-value <0.05).

The alternative splicing events (ASEs) and regulated alternative splicing events (RASEs) between the samples were defined and quantified by using the ABLas pipeline, as previously described (37). In brief, eight types of ASEs were identified, based on the splice junction reads. The eight possible types of ASE included Cassette exon (CassetteExon), Exon skipping (ES), Mutual exclusive exon skipping (MXE), A5SS, A3SS, the MXE combined with alternative 5′promoter (5pMXE) combined with alternative polyadenylation site (3pMXE) and intron retention.

Having identified the ASEs in each RNA-seq sample, Fisher's exact test was selected to determine statistical significance, using the alternative reads and model reads of the samples as input data. The altered ratio of alternatively spliced reads and constitutively spliced reads between compared samples was defined as the RASE ratio. The P-value <0.05 and RASE ratio >0.2 were set as the thresholds for detection of RASEs.

To determine the validity of the RNA-seq data, qPCR was performed for selected DEGs and normalized with the housekeeping gene GAPDH. Primers are designed in exon regions, and sequences are presented in Table I. The same RNA samples for RNA-seq and RNA samples isolated from CK18 knockdown in HepG2 (using the same shRNAs) were used for qPCR. The PCR conditions consisted of denaturing at 95°C for 10 min, 40 cycles of denaturing at 95°C for 15 sec, and annealing and extension at 60°C for 1 min. PCR amplifications were performed in triplicate for each sample.

Concurrently, a qPCR assay was used to analyze ASEs. The primers for detecting the pre-mRNA splicing are presented in Table I. To detect one of the alternative isoforms, one primer was designed in the alternative exon and an opposing primer was designed in a constitutive exon. Other alternative isoforms were detected using a boundary-spanning primer for the sequence encompassing the exon-exon junction, with the opposing primer in a constitutive exon.

CK18 knockdown and control Hela cells were lysed in RIPA buffer containing 50 mM Tris-HCl (pH 7.4), 150 mM NaCl, 1.0% deoxycholate, 1% Triton X-100, 1 mM EDTA and 0.1% SDS. Following centrifugation of the homogenate (20,000 × g, 15 min) the supernatants were used for western blotting. Protein concentrations were measured using the BCA protein assay (Beyotime Institute of Biotechnology, Nanjing, China) with bovine serum albumin (BSA) as a standard. Equal amounts (20 µg/lane) of protein samples were loaded on 12% SDS-PAGE gel for separation and then transferred to a polyvinylidene difluoride (PVDF) membrane. The membrane was incubated with blocking solution (1X TBS; 0.05% Tween-20; 5% non-fat milk) at room temperature for 1 h and incubated overnight with primary antibodies raised against CK18 (dilution 1:1,000; cat. no. A0389), CASP8 (dilution 1:1,000; cat. no. A0215) and GAPDH (dilution 1:1,000; cat. no. AC027) (all from ABclonal, Wuhan, China). Immunoreactive proteins were detected using an ECL chemiluminescence system (Clinx Science Instruments Co., Ltd., Shanghai, China) with default settings and GAPDH as the normalized control.

RNA-seq data presented in this study have been deposited in the Gene Expression Omnibus of NCBI and are accessible through GEO series accession no. GSE119255.

Experimental data are presented as the mean ± standard deviation of at least three experiments. Statistical analyses were performed with SPSS software (version 17.0; SPSS, Inc., Chicago, IL, USA). Significance of differences was evaluated with Student's t-test when only two groups were compared. P<0.05 was considered to indicate a statistically significant difference.

The expression of CK18 was examined in HeLa cells transduced with different shRNAs (shRNA1, shRNA2 and shRNA3) or an empty vector by RT-qPCR. Compared with cells treated with the empty vector mRNA expression was effectively reduced in cells treated with CK18-shRNAs. The effective shRNA (shRNA1) also led to a significant reduction in the protein level (Fig. 1A). Therefore, CK18-shRNA1 was used in subsequent experiments in order to knockdown CK18 expression in HeLa cells. Knockdown of CK18 (CK18 KD) led to a significant increase in cell proliferation and a significant decrease in cell apoptosis (Fig. 1B and C).

CK18 KD and control cells were used to construct cDNA libraries for sequencing on an Illumina HiSeq 4000 platform. Two biological replicates were used and a total of 76.7±4.9M 150 nucleotide paired-end raw reads per sample were obtained. The raw reads were filtered by removing low-quality reads and reads containing N and adaptor sequences, leaving 73.5±4.7M clean reads for downstream bioinformatics analysis. The clean reads were then mapped onto the human GRCh38 genome using TopHat2: 77.99-81.98% were aligned and ~86.77-95.05% were uniquely mapped (Table II). Quantification of genes and transcripts was reassessed using Cufflinks (38), to compare gene expression patterns across individuals. FPKM values were then calculated. There were 22,677 genes expressed in the RNA-seq data (available upon request). Effective KD of CK18 was further confirmed by in parallel RNA-seq analysis (Fig. 2A).

Using criteria of an absolute fold change ≥2 and FDR <0.05 with the edge R package (35), 294 upregulated and 299 downregulated genes related with CK18 KD (data not shown) were identified. A volcano plot was constructed to display the significantly expressed genes that were associated with CK18 KD (Fig. 2B). The heatmap demonstrated distinguishable transcription profiles between CK18 KD and the control groups (Fig. 2C).

The two replicate RNA-seq datasets were divided into two groups (case 1 with control 1; case 2 with control 2). Using the criteria of an absolute fold change ≥2 with the EdgeR software package, 174 co-upregulated and 146 co-downregulated genes overlapping in two groups associated with CK18 KD were identified. A Venn-diagram demonstrated the upregulated and downregulated genes between the two different analytical strategies (Fig. 2D).

GO enrichment analysis was performed to further explore the biological function of DEGs with two different analytical strategies. All three ontologies of GO analysis, molecular function, cellular component and biological process, were obtained. The top biological process terms of the GO analysis that involved upregulation or downregulation of genes following CK18 KD are presented in Fig. 3A (repetition) and B (simple pair). Most of the terms overlapped; the upregulated genes were mainly associated with processes involved in multicellular organismal development, signal transduction, cell differentiation and the apoptosis process. The downregulated genes, conversely, were mainly associated with negative regulation of transcription from RNA polymerase II promoter, innate immune response, signal transduction, and regulation of transcription, DNA-dependent. Based on the KEGG analysis, the top ten pathways involved in up- and downregulated DEGs were also presented in Fig. 3. Coinciding with the GO analyses, many pathways were associated with the immune response and apoptosis. The B-cell receptor signaling pathway, the T-cell receptor signaling pathway and the Jak-STAT signaling pathway were all significantly enriched in downregulated gene sets.

Regulation of alternative splicing (AS) of CK18 in the transcriptome sequencing data was also explored. Between 54.2 and 57.32% of the uniquely mapped reads from CK18 KD and control samples were junction reads (Table II). A total of 58.9% of annotated exons (216,356 out of 367,321) were detected. For splicing junctions, 147,540 annotated and 86,961 novel junctions were detected using TopHat2.

ASEs were analyzed using ABLas software to investigate global changes in the occurrence of alternative splicing. We detected 15,786 known ASEs in the model gene that we designated in the reference genome, along with 40,946 novel ASEs, excluding intron retention (Table III).

We identified 263 high confidence RASEs by using a stringent cutoff of P≤0.05, along with changed AS ratio ≥0.2. CK18-regulated ASEs, which included 35 alternative 3'splice sites (A3SS), 40 alternative 5'splice sites (A5SS), 44 examples of exon skipping (ES) and 24 cassette exons, 12 mutually exclusive 5′UTRs (5pMXE), four mutually exclusive 3′UTRs (3pMXE), six mutually exclusive exons (MXE), 84 examples of intron retention (IR), eight examples of alternative 5'splice site & exon skipping (A5SS & ES) and six examples of alternative 3'splice site & exon skipping (A3SS & ES), are summarized in Fig. 4A. These results indicated that CK18 KD had a broad influence on splicing. Coupled to the transcription data, there were no significant expression changes in alternative splicing genes (Fig. 4B). The alternative spliced genes that were identified by GO analyzed were, however, enriched in liver development, positive regulation of I-κB kinase/NF-κB cascade, glucose metabolic process, positive regulation of NF-κB transcription factor activity, mRNA processing and small molecule metabolic process. These results were similar to those obtained by the transcriptional GO analysis, further confirming the important regulatory role of CK18 in these biological functions. The detailed results of the GO and KEGG pathway analyses are presented in Fig. 4C. It was observed that alternative splicing of CASP8 and FAS, the known upstream genes regulating the production of CK18, was under the regulation of CK18.

To confirm the important regulatory function of CK18 on both gene expression levels of alternative splicing in HeLa and HepG2 cell lines, we validated certain DEGs and ASEs that were important in immunity and apoptosis by RT-qPCR. The ASEs detected by PCR primer pairs (Table I) were designed to amplify both long splicing isoforms and the short splicing isoforms in the same reaction.

Most of the DEGs associated with immune function and apoptosis that we validated were in agreement with RNA-seq results (Fig. 5A and B). For CASP8 gene expression, however, there was no significant change in RNA-seq data, whereas the qPCR and western blotting results revealed an evident decrease in CK18 KD cells (Fig. 5C). Two important ASEs were located in FAS and CTNNB1, which have been well validated as key genes in apoptotic pathways (Fig. 6).