We acquired RNA sequencing data from the TCGA database (https://tcga-data.nci.nih.gov) [2]. Bam files were analyzed using VarScan2 [16]. We compared data from normal and tumor samples with the following variant calling parameters: min-coverage, 8; min-var-freq, 0.2; mina-avg-qual, 20; and min-reads, 2. Variant call files were annotated with ANNOVAR [17]. We selected the positions that were converted from A to G, nonsynonymous single nucleotide variants (nsSNVs), after filtering somatic mutations that were enrolled in TCGA. We also removed all of the data from dbsnp137, excluding cDNA data and the positions within similar genomic regions.

To compare GLI1 RNA editing between tumor and normal tissue samples, 10 human colorectal tumor tissues and the matched normal tissues were obtained from the tissue bank of Seoul National University Hospital (Seoul, Korea).

Genomic DNA (gDNA) from cell lines was isolated using the QIAMP DNA mini kit (Qiagen, Valencia, CA) according to the manufacturer’s instructions. Total RNA was extracted from each cell line using an RNeasy mini kit (Qiagen) according to the manufacturer's instructions, after which 1 μg of total RNA was used to synthesize cDNA with random hexamers and ImProm-II reverse transcriptase (Promega, Madison, WI).

RNA editing sites were amplified by polymerase chain reaction (PCR) using target-specific primers in gDNA or cDNA templates. Information regarding the primer is described in S1 Table. Purified PCR products were sequenced with a BigDye Terminator v3.1 cycle sequencing kit (Applied Biosystems, Carlsbad, CA), then analyzed with a 3730 ABI capillary electrophoresis system (Applied Biosystems).

We designed GLI1 target-specific primers for pyrosequencing analysis. The reverse primer was biotin-labeled. Single-stranded biotinylated PCR products were processed for pyrosequencing analysis according to the manufacturer’s standard protocols (PyroMark Q96 ID, Qiagen). The primer sequences were as follows: forward, 5’-GACCGTCCTGCTCCAGCTAG-3’; reverse, 5’-CCCAACTTCTGGCTCTTCCTGT-3’; sequencing, 5’-AATGCTGCCATGGAT-3’.

The inosine chemical erasing (ICE) assay is composed of four steps: (1) RNA cyanoethylation, (2) cDNA synthesis by reverse transcription, (3) PCR amplification, and (4) direct sequencing. We cyanoethylated total RNA from cancer cells containing edited transcripts for 15 minutes at 70°C alongside controls lacking cyanoethylation. The ICE assay was performed as previously described [18].

SNU-254 and SNU-81 human CRC cell lines were obtained from the Korea Cell Line Bank [19]. The cell line was grown in RPMI-1640 with 10% fetal bovine serum and gentamicin (10 μg/mL) at 37°C in a 5% CO2-humidified atmosphere. siRNA was purchased from Mbiotech Inc. (Seoul, Korea). A total of 20 nM siRNA was applied to cells in culture using G-fectin (Genolution, Seoul, Korea) according to manufacturer's standard protocols.

Cultured cells were washed with ice-cold phosphate-buffered saline and lysed with lysis buffer (50 mM Tris-HCl, pH 7.5, 1% NP-40, 0.1% sodium deoxycholate, 150 mM NaCl, 50 mM NaF, 1 mM sodium pyrophosphate, 1 mM EDTA, and protease/phosphatase inhibitors). Lysates were transferred to new eppendorf tubes after centrifugation at 13,000 rpm for 20 minutes. Protein concentrations were quantified with a Bicinchoninic Acid Protein Assay Reagent (Pierce, Rockford, IL) according to manufacturer's instructions. Samples containing equal quantities of total proteins were resolved on sodium dodecyl sulfate–polyacrylamide denaturing gel, then transferred to nitrocellulose membranes. The membranes were subsequently incubated in blocking solution containing 1% skim milk and 1% bovine serum albumin for 1 hour at room temperature, then probed overnight at 4°C with primary antibodies (anti–adenosine deaminase acting on RNA 1 [ADAR1], Abcam, Cambridg, MA; anti-actin, Santa Cruz Biotechnology, Santa Cruz, CA).

The statistical significance of the results was calculated by unpaired Student’s tests, with p-values of < 0.05 considered to be significant.

We downloaded 39 paired tumor and nontumor RNA transcriptome bam files from the TCGA database to analyze RNA editing patterns. The amount of total reads can influence the number of editing sites [20]; therefore, we evaluated only those samples with similar reads between normal and tumor tissues to reduce bias related to total reads. We utilized strict filtering conditions to find RNA editing candidates, using the VarScan2 software to identify RNA editing sites. VarScan2 was originally developed to identify somatic mutations and copy number alterations in exome sequencing data, but it can also be applied to analyze RNA-seq data. We set the analysis standards as follows: total reads, more than 8; editing frequency, 20; quality score, 20; and edited reads, 3. We then only selected nonsynonymous A to G conversions based on changes in AA in the refGene column, which was generated by ANNOVAR. We removed the positions enrolled in dbsnp137 except from those that originated from cDNA, similar genomic regions, and somatic mutations from the TCGA database (S2 Fig.). Overall, we identified 939 editing sites from 39 paired samples, including recurrent positions, ranging from 12 to 42 sites in each sample, with an average of 24 (Fig. 1A). We also investigated the percentage of hypo- and hyper-RNA editing within the samples. We defined hyperediting as a situation in which the percentage of RNA editing was 10% greater in tumors than adjacent normal tissue. There was no trend in hypoediting or hyperediting within the samples investigated. Among the 39 samples, 26 were classified as hyperedited and 13 were hypoedited (Fig. 1B).

We then focused on the recurrent hyper- or hypo-RNA editing sites in CRC. First, we listed the positions that were detected in more than two samples. The total number of recurrent detected editing positions was 39. Although we used relatively strict filtering standards, the possibility of false-positives remained; therefore, we employed additional filtering steps using the IGV AND PLAST programs. After these rounds of filtering, we found 10 recurrent editing positions in nine genes (Table 1).

To identify RNA editing biomarkers for CRC, we investigated RNA editing by dividing edited messages into two categories: hypoediting, when the level of RNA editing is higher in normal than tumor samples; and hyperediting, when the level of RNA editing is higher in tumor samples. The results revealed four positions in three genes that were hypoedited in cancer. Hypo-RNA editing of IGFBP7 (c.A284G, p.K95R) was detected in 33 samples. The other hypo-RNA editing position of IGFBP7 (c.A232G, p.R78G) was detected in five samples (Fig. 2A and B), COPA, Coatomer Protein Complex, Subunit Alpha, was detected in four samples (Fig. 2C), and ZNF358, Zinc Finger Protein 358, was detected in two samples (Fig. 2D). These messages were significantly less frequently edited in cancer. Moreover, PDLIM, NEIL1, SRP9, and APMAP had both hyper- and hypo-RNA editing in 39 samples.

GLI1 RNA editing was detected in five samples, and was increased in colon cancer compared with normal tissue (Fig. 3A). We validated this position using Sanger sequencing (Fig. 3A) based on the pyrosequencing results of CRC cell lines (S3 Fig.). The guanosine peak was detected in cDNA, but not gDNA. After identifying GLI1 701G nsSNV, we confirmed the hyperedited pattern in CRC using Sanger sequencing in our tissue samples (Fig. 3C). There were increased guanosine peaks in tumor samples. To quantitatively analyze A-to-I RNA editing of GLI1, we conducted pyrosequencing in 10 paired samples of CRC and normal samples to measure the edited cDNA sequence. The extent of RNA editing of GLI1 was 29.3% in tumor tissues, while it was only 15.1% in normal tissues (Table 2, Fig. 3D). When compared with normal tissues, the level of RNA editing was increased by about 2-fold in tumor tissues. We also analyzed the clinical significance of GLI1 RNA editing in a TCGA cohort (n=39). More GLI1 was expressed in GLI1 RNA editing samples than in samples without GLI1 RNA editing (S4 Fig.). High expression or GLI1 RNA editing showed a significant relationship with the clinical prognosis (S5 Fig.). The clinical significance of GLI1 RNA editing merits further investigation in a larger patient population. Moreover, we found that the age, sex, primary site, histology, stage, microsatellite stable/microsatellite instable status, and KRAS mutation was not significantly related to GLI1 RNA editing (S6 Table). These data suggest that GLI1 is hyperedited in CRC.

To confirm the reliability of our RNA editing discovery methods, we validated editing events using Sanger sequencing of gDNA and cDNA from the CRC cell lines. The cell lines that contain RNA editing were selected based on our RNA sequencing data. RNA editing of SRP9 and NEIL1 was validated because they were also detected in our cell line RNA sequencing data. Sanger sequencing confirmed the presence of a G peak in the cDNA of SRP9 and NEIL1 (Fig. 4A). The ICE assay is another method to validate A-to-I RNA editing, as cyanoethylated inosine will not be converted to guanosine by reverse transcriptase. We conducted an ICE assay to validate RNA editing of SRP9 and NEIL1. From these experiments, we found that the guanosine peaks were significantly decreased in the treated versus the nontreated samples in both positions (Fig. 4B). These data confirm the existence of A-to-I RNA editing from our analysis and, therefore, exclude the possibility of gDNA contamination.

ADAR family enzymes mediate A-to-I RNA editing. To investigate whether ADAR enzymes were involved in the RNA editing sites we discovered in CRC, we conducted RNAi-mediated ADAR knockdown. The expression of ADAR was reduced to approximately 40% of the wild type mRNA and protein levels (Fig. 4D). We next performed Sanger sequencing to determine the ratio of A/G. The editing level was reduced in both SRP9 and NEIL1 editing positions (G: 45.6%→31.0% in SRP9, G:74.3%→43.9% in NEIL1) in si-RNA ADAR treated cells (Fig. 4C). Although there were limitations to these samples, the results show that the identified RNA editing candidates exist, and are affected by ADAR enzyme knockdown.