Evaluation of the Interplay between the ADAR Editome and Immunotherapy in Melanoma

Abstract

Background: RNA editing is a highly conserved posttranscriptional mechanism that contributes to transcriptome diversity. In mammals, it includes nucleobase deaminations that convert cytidine (C) into uridine (U) and adenosine (A) into inosine (I). Evidence from cancer studies indicates that RNA-editing enzymes promote certain mechanisms of tumorigenesis. On the other hand, recoding editing in mRNA can generate mutations in proteins that can participate in the Major Histocompatibility Complex (MHC) ligandome and can therefore be recognized by the adaptive immune system. Anti-cancer treatment based on the administration of immune checkpoint inhibitors enhance these natural anti-cancer immune responses. Results: Based on RNA-Seq datasets, we evaluated the editome of melanoma cell lines generated from patients pre- and post-immunotherapy with immune checkpoint inhibitors. Our results reveal a differential editing in Arthrobacter luteus (Alu) sequences between samples pre-therapy and relapses during therapy with immune checkpoint inhibitors. Conclusion: These data pave the way towards the development of new diagnostics and therapies targeted to editing that could help in preventing relapses during immunotherapies.

1. Introduction

The diversification of messenger RNA (mRNA) sequences from genomic DNA relies on posttranscriptional mechanisms, such as alternative splicing and RNA editing [1,2]. In mammals, RNA editing involves the deamination of adenosines (A) or cytidines (C) generating in situ inosine (I) or uridine (U), respectively. While adenosine deaminase acting on RNA (ADAR) [3] and adenosine deaminase acting on tRNA (ADAT) [4] deaminate A residues, enzymes from the APOBEC family [5] deaminate C residues. In the following, we focus on ADAR activities, as A-to-I modifications globally represent 97% of RNA editing events [6], in mammals and especially in humans. A-to-I editing events can occur in coding and non-coding sequences of mRNA [7]. In general, they are prominent in double-stranded RNA stretches formed by inverted non-coding repeats such as Arthrobacter luteus (Alu) and long interspersed element (LINE) located in mRNA untranslated regions (UTRs) and introns [8,9,10]. RNA editing in non-coding sequences can influence alternative splicing, nuclear retention, and transcript degradation (e.g., recognition by miRNA) [11]. In addition, RNA editing can affect sites in the protein coding region of mRNAs, leading to potential amino acid changes, known as recoding editing [12]. In humans, recoding editing is rare while A-to-I in Alu elements is abundant and accounts for 97% of all available events [13]. Both ADAR1 and ADAR2 can perform recoding or non-recoding editing. A third ADAR protein called ADAR3 is expressed only in the brain and its deaminase activity has not been yet proven [7].

A-to-I RNA editing is essential in maintaining cellular homeostasis [14,15,16] and has been implicated in several diseases ([17] reviewed in [18,19,20]). In particular, the disruption of the controlled expression of ADAR1 and ADAR2 has been shown to contribute to cancer pathogenesis. Based on current reports, recoding and non-coding editomes with hypo- or hyper-editing levels appear to be dependent on cancer-types and genes [21,22]. For example, the DNA base excision repair glycosylase enzyme NEI-like protein 1 (NEIL1), encoding antizyme inhibitor 1 (AZIN1), Ras homologue family member Q (RHOQ), and protein tyrosine phosphatase non-receptor type 6 (PTPN6) are found hyper-edited in cancer compared to healthy tissues, while gamma-aminobutyric acid type A receptor alpha 3 subunit (Gabra3), glutamate receptor subunit B Glur-B (also known as GRIA2), and insulin-like growth factor-binding protein 7 (IGFBP7) are found to be hypo-edited in cancer compared to normal tissue. These editing events have been reported to impact protein functions. Other recoding editing, such as filamin B (FLNB), cyclin I (CCNI), coatomer protein complex subunit α (COPA), and the component of oligomeric Golgi complex 3 (COG3) have not been functionally characterized. Non-coding editomes in cancer can also be different when compared to the corresponding healthy tissues, such as with hypo-editing for melanoma and glioma or hyper-editing (of microRNA) in non-small cell lung carcinoma (NSCLC). This editome is involved in cancer development, progression, invasion, metastatic potential, recurrence, and resistance (reviewed by [18]).

ADARs are also associated with immunity in several ways. First, the elongated form of ADAR1 is located in the cytosol and is induced by type I interferon, while constitutive ADAR2 and the short form of ADAR1 are nuclear. Second, by editing dsRNA, ADARs avoid the stimulation of innate immune responses by endogenous dsRNA [23] and thereby, lowering ADAR1 increases tumor inflammation [24,25]. Third, adaptive immunity is impacted by ADARs; a recent report shows that peptides containing amino acids generated through ADAR-recoding events are human leukocyte antigen (HLA) ligands. In particular, CCNI-edited peptides act as cancer antigens capable of activating tumor infiltrating lymphocytes (TILs) and thereby mediate cancer cell death in melanoma [26]. Thus, ADAR-recoding can impact the Major Histocompatibility Complex (MHC) ligandome and thereby the specific anti-cancer T-cell response.

Although ADARs affect immunity, studies investigating the potential relationship between editome and cancer immune evasion are lacking. Melanoma is the most aggressive form of skin cancer and is currently best treated through the administration of immune checkpoint inhibitors: monoclonal antibodies blocking negative signals generated by cytotoxic T-lymphocyte-associated protein 4 (Ipilimumab) or programmed cell death-1 (e.g., Pembrolizumab or Nivolumab) in activated T lymphocytes [27]. Since such therapies rely on the potential for the immune system to fight cancer cells, and since ADARs are involved in the immune response in multiple ways (see above), we used RNA-Seq datasets to analyze the editome of melanoma cell lines made from tumors obtained before immunotherapy and from Ipilimumab and Pembrolizumab-resistant (relapsing) metastasis.

2. Materials and Methods

2.1. Patients and Cell Lines

The study was approved by the Kantonale Ethikkommission under approval number KEK-ZH Nr. 2014-0425, EK647 and EK800. After informed consent was given by the patients, biopsies were used to generate melanoma cell lines as described [28]. DNA was used to prepare a customized target library using the Nimblegen SeqCap EZ kit (Nimblegen). Sequencing was performed on an Illumina Hiseq 4000 machine at the Functional Genomics Center of the University of Zürich (FGCZ). For targeted sequencing, we sequenced 0.125 lanes per sample, paired-end, with 150 bp reads [28].

Cell line purity was estimated based on the mutant allele frequency (MAF) calculated as follows: MAF = mutant copies/(wild-type copies + mutant copies). Mutations in BRAF (V600E) and NRAS (Q61R, Q61K, Q61L) in melanoma are reported to be heterozygous, and thus cell line purity was calculated as 2 × MAF. RNA was extracted with the QIAGEN RNeasy kit. RNA capture was performed with the TruSeq RNA library Prep Kit v2 (Illumina). RNA sequencing was sequenced at 125 bp paired end on a HiSeq4000 at the FGCZ.

Datasets for 14 patients were available for analysis. The samples were grouped by treatment: 12 samples were assigned to the group not yet subjected to Ipilimumab and Pembrolizumab therapy (12 cell lines taken from 10 patients) and 11 samples were assigned to the relapsed group (11 cell lines taken from 4 relapsing patients during treatment by Ipilimumab and Pembrolizumab). Those relapses happened during immunotherapy as they appeared within less than one month after the last antibody injection. For two patients, we collected cell lines before immunotherapy and at the point of relapse.

Raw sequencing reads and RNA counts can be retrieved with GEO accession number GSE162095 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE162095).

2.2. Bioinformatics and Statistics

Gene-level counts were quantified from pair-end reads aligned to the GRCh38 genome using the quantMode option implemented in STAR (it enables the aligner to count number reads per gene while mapping) [29]. A differential gene expression analysis was performed using Qlucore software v 3.5. and using the R script for limma-voom [30]. Significant differential expression was set to p < 0.05 and the false discovery rate q < 0.05. Gene expression levels for the ADAR genes were obtained from our USZ melanoma website http://phil.shinyapps.io/ZMCE. Differences between the two groups of patient samples were analyzed through an unpaired t-test conducted in GraphPad where p < 0.05 values were considered significant. RNA editing events were quantified with ExpEdit [31,32,33,34], a web interface based on REDItools [31]. The program explores known human RNA editing positions annotated in DARNED (approximately 333,215 edited sites, of which 221,595 are A-to-G edited sites) from FASTQ RNA-seq input files. The default settings displayed by the interface were selected for analysis: the qPhred score was set to ≥10 (base quality score), and the mapping quality score was set to ≥20 [26]. The obtained files were adapted for further analysis with Qlucore software v 3.5 focusing on editing events of ≥10%. Recoding and non-recoding editing was analyzed separately with the exclusion of known Single Nucleotide Polymorphism (SNPs) from dbSNP and of any nucleotide changes that did not resemble A to I (G) [9]. Principal component analysis (PCA) plots and heat maps applying hierarchical clustering were generated using Qlucore software. The Alu editing events were analyzed with an unpaired t-test with the Benjamini and Hochberg false discovery rate multiple testing correction where significantly different Alu editing was measured with a p value of <0.05 and a q value of <0.1. Alu editing indexes were calculated as the ratio of the sum of the reads with edited sites (I = G in DARNED) to the total number of reads [22]. Correlation plots were analyzed in GraphPad applying Pearson correlation, where a p value of <0.05 was considered significant.

3. Results

We analyzed 23 RNA seq; 12 datasets were collected from cell lines made from biopsies before the administration of immune checkpoint inhibitors (“before IT”), and 11 datasets were collected from cell lines made from biopsies of relapsing melanoma in patients treated with Ipilimumab and Pembrolizumab (“relapsed”) (

Table 1. Patients and corresponding cell lines. The table presents the patients’ demographics and treatments as well as the characteristics of the cell lines.

	Gender	Tumor	Excision Place	Cell Lines	Treatment	Immunotherapy	Duration	Cell Line Purity
Patient 1	M	Metastasis	Lymph node	MM130820	Before IT			87%
Patient 2	F	Metastasis	Trunk	MM150325	Before IT			98%
	F	Metastasis	Trunk	MM150604	Relapsed	Ipilimumab	3 weeks	74%
Patient 3	F	Metastasis	Lymph node	MM140905	Before IT			99%
Patient 4	M	Metastasis	Back	M130107	Before IT			59%
	M	Metastasis	Lymph node	M121102	Before IT			65%
Patient 5	F	Metastasis	Lung	M121008	Before IT			20%
Patient 6	F	Metastasis	Lymph node	MM130626	Before IT			60%
Patient 7	F	Metastasis	Leg	MM130106	Before IT			57%
	F	Metastasis	Leg	MM131205	Relapsed	Ipilimumab	3 weeks	78%
	F	Metastasis	Leg	MM131206	Relapsed	Ipilimumab	3 weeks	65%
Patient 8	M	Metastasis	Trunk	MM150405	Before IT			58%
	M	Metastasis	Lymph node	M150404	Before IT
Patient 9	F	Metastasis	Skin	MM130926	Before IT			50%
Patient 10	F	Metastasis	Lymph node	M980513	Before IT			83%
Patient 11	M	Metastasis	Arm	MM140902	Relapsed	Ipilimumab	3 months	100%
	M	Metastasis	Arm	M130830	Relapsed	Ipilimumab	3 months	100%
Patient 12	M	Metastasis	Back	MM130434	Relapsed	Ipilimumab	3 months	77%
Patient 13	M	Metastasis	Trunk	MM121224	Relapsed	Ipilimumab	3 months	75%
Patient 14	M	Metastasis	Skin	M130420	Relapsed	Ipilimumab	10 days	100%
	M	Metastasis	Liver	M130421	Relapsed	Ipilimumab	10 days	100%
	M	Metastasis	Testis	M130425	Relapsed	Ipilimumab	10 days	100%
	M	Metastasis	Lung	MM130427	Relapsed	Ipilimumab	10 days	97%

). Clinical features of the melanoma patients from which cell lines were derived are indicated for each patient in Table 1 and for both groups in

Table 2. Demographics and characteristics of the cell lines for both groups.

Characteristic	Relapsed N = 11 1	Before IT N = 12 1
Age at sampling date	48 (18–60)	61 (32–86)
Sex
F	3 (27%)	7 (58%)
M	8 (73%)	5 (42%)
Mutation
BRAF	2 (18%)	7 (58%)
BRAF and NRAS	1 (9.1%)	0 (0%)
NRAS	8 (73%)	5 (42%)
Excision place
Arm	2 (18%)	0 (0%)
Back	1 (9.1%)	1 (8.3%)
Leg	2 (18%)	1 (8.3%)
Liver	1 (9.1%)	0 (0%)
Lung	1 (9.1%)	1 (8.3%)
Lymph node	0 (0%)	5 (42%)
Skin	1 (9.1%)	2 (17%)
Testis	1 (9.1%)	0 (0%)
Trunk	2 (18%)	2 (17%)

¹ Statistics presented: median (minimum-maximum); n (%).

. Our methods to establish cell lines were optimized so that melanoma cells were more effectively retrieved and that the cell culture represented the full range of in vivo tumor heterogeneity [28]. It is not possible to verify that the cell lines will have similar gene expressions or editing profiles than tumor cells in the biopsies, however as all cell lines analyzed here were produced in the exact same conditions, we would expect that if the cell culture derivation would affect gene expression or editing, it would be similar for all, thus not affecting the comparison between the two groups. From a multivariant analysis of limma and voom, no gene was found to be significantly up- or downregulated in one group or the other when p values of less than 0.05 and q values of less than 0.05 were used (Supplementary Figure S1A and Supplementary Table S1). In addition, no clear clustering of the two groups can be evidenced using PCA (Supplementary Figure S1B). Thus, we extended our transcriptome analysis by exploring the role of ADARs in checkpoint inhibitor immunotherapies against melanoma (Supplementary Table S2).

3.1. Expression of the ADAR Genes before IT Versus in Relapsing Tumours

We evaluated changes in gene expression between the two groups for the Adar1 and Adarb1 genes corresponding to the ADAR1 and ADAR2 enzymes, respectively. The analysis shows that ADAR1 gene expression (long and short forms of ADAR1 combined) remained unchanged between the groups (Figure 1A), while ADAR2 was significantly enhanced in the “relapsed” group (Figure 1B) based on an unpaired t-test. This result was unexpected, as ADAR2 was previously reported to act as a tumor suppressor [35,36] while ADAR1 has been preferably deemed a tumor oncogene (reviewed by [18]). Thus, the role of ADARs in relapse during therapy by immune checkpoint inhibitors may not be the same as the role of such enzymes in the natural course of melanoma development.

3.2. Recoding Editing

Clinically relevant recoding editing in cancer has been reported. Corresponding levels of editing are shown in Supplementary Figure S2 and Supplementary Table S3. ADAR2 is responsible for most recoding editing events including COG3_I635V and COPA_I164V [37]. While ADAR2 is upregulated in “relapsed”, the two forms of recoding editing were not significantly different between “before IT” and “relapsed” (Figure 2A). These forms of recoding editing show a strong positive correlation with the expression of ADAR2, confirming their ADAR2 dependence (Figure 2B). The recoding editing of COG3 and COPA was first identified as clinically relevant in more than one form of cancer in Han et al. [21]. The editing of COG3 and COPA has been shown to correlate with shortened progression-free survival in renal clear cell carcinoma, and high levels of COG3 editing have been associated with resistance to fluorouracil and austocystin D while high levels of COPA editing have been significantly associated with resistance to austocystin D and lapatinib [37]. Meanwhile, it has been reported that RNA editing at both CCNI R75G and COPA I164V generates MHC-associated epitopes [26]. Our data suggest that in patients treated with immune checkpoint inhibitors, MHC epitopes derived from edited COPA are not strongly immunogenic, as otherwise they would be downregulated when the immune response is enhanced by blockade of immune checkpoint inhibitors.

3.3. Non-Recoding Editing in Alu Elements

The global noncoding editome represented by the Alu editing indexes shows no differences between the two groups (Figure 3A and Supplementary Table S4). Additionally, in line with previous reports, we find a positive correlation between ADAR1 expression and Alu editing indexes (Figure 3B), confirming ADAR1’s major contribution to the noncoding editome. However, a focus on non-recoding single editing sites based on a qPhred score of ≥10 indicates Alu hyper-editing in relapsed samples of Gap junction gamma-1 protein (GJC1) at site 42877641 (Figure 3C and Supplementary Tables S5 and S6). Interestingly, GJC1 Alu editing frequency at this site shows no correlation with ADAR1 expression, but it correlates positively with ADAR2 expression (Figure 3D). This editing of Alu elements has been observed for nuclear-retained Cat2 transcribed nuclear RNA [38]. Furthermore, Alu editing in GJC1 is reported as clinically relevant, showing differential editing levels across tumor subtypes and tumor stages and correlating with patients’ overall survival rates [21]. However, further studies must be conducted to determine the contributions of this Alu editing to cancer progression and to responses to therapy.

3.4. Principal Component Analysis

To obtain an overview of potential differences in editing patterns observed between “Before IT” tumors and “Relapsed” during treatment with checkpoint inhibitors, we applied a principal component analysis. It was done using the two group comparison statistical analysis (t-test) provided by the Qlucore software and the following parameters: filtering by variance 0.15, p value was 0.05, and q value was 0.775. The outcome by the software was 57 Alu editing sites in 57 genes for the sample representation in a 3D PCA plot and heatmap representation of Alu editing sites (Figure 4 and Supplementary Table S7). Editing events were not normalized as a single editing site was used as a unique identifier for variable analysis and not the complete gene with all the editing events. From the PCA plot including the 57 Alu editing events, we can clearly separate the two groups of patients based on the corresponding cell lines. This separation is confirmed in a hierarchically clustered heatmap showing that Alu signatures can differentiate cell lines from patients who relapse from cell lines taken from patients before IT (Figure 4A,B). Conversely, no recoding editing signatures could be identified to separate the two groups (Supplementary Figure S2).

Collectively, our study reveals an editing signature in Alu elements that characterizes tumors relapsing during treatment with immune checkpoint inhibitors (Figure 3C and Figure 4). Surprisingly, this special signature may be associated with the higher expression of ADAR2 (at least for GJC1, Figure 3D). While ADAR2 has been envisioned as a tumor suppressor, its increased expression in “relapses” would indicate that under treatment with immune checkpoint inhibitors, inhibiting ADAR2 could help prevent tumor recurrence. Some targets of ADAR2, mostly recoding editing events may be of relevance, but they could not be identified in the present study (non-significant changes in the recoding editing of COPA and COG3 and no grouping of recoding events, Supplementary Figure S2). Meanwhile, the known recoding editing that generates MHC epitopes that are recognized by anti-cancer T-cells was not downregulated, indicating the weak significance of these epitopes for the immune control of cancer. The combination of several recoding and non-coding editing events (Alu editing of GJC1) eventually mediated by ADAR2 may be responsible for the potential advantages that the overexpression of this gene provides for recurrence during treatment with immune checkpoint inhibitors. Although RNAseq data from whole tumor tissues are publicly available (for example from Liu et al. [39]), an analysis of A to G editing and its eventual correlation with clinical outcome would require the generation of tumor cell lines. Indeed, in whole tumor tissues, RNA sequences are originating not solely from tumor cells but also from non-tumor stroma cells, immune cells, etc. Thus, editing analysis of these sequencing files would not provide information exclusively on melanoma cells but would primarily reflect the heterogeneity of the tumor samples (percentages of immune cells, non-tumor cells, blood vessels, etc.) that varies from one sample to another. We foresee that our results will galvanize the analysis of A to G editing in tumor cell lines made from patients with a precisely known cancer history and will enable the identification of further correlations between editing and cancer outcome. Our study does not address the biological or biochemical phenomenon that connects the Alu-editing signature to the immunotherapeutic treatment. It however is the very first study that shows a correlation between RNA editing in melanoma and clinical outcome. We foresee that based on those observational results, further studies can be undertaken to more precisely decipher the mechanisms leading to the differential RNA editing in tumor cells during immunotherapy of cancer.

4. Conclusions

Our study points to differential editing in Alu sequences between melanoma cell lines obtained before therapy and melanoma cell lines made from relapses during treatment with immune checkpoint inhibitors. Those findings may be of relevance for diagnostic and prognostic tools as well as for the development of drugs or treatments that may lower the risks of relapses during therapy with immune checkpoint inhibitors.