Global analysis of A-to-I RNA editing reveals association with common disease variants

Statistics

R v.3.2.2 was used for all statistical analyses (R Core Team, 2015). If not otherwise specified, P-values were calculated with Welch’s t-test. All statistical tests were two-sided.

Data sources

RNA-seq and genotype data (generated with the Illumina Infinium OmniExpressExome-8 chip) in this study originated from samples collected of up to 7 tissues (artery aorta, internal mammary artery, whole blood, subcutaneous fat, visceral fat, liver, and skeletal muscle) and 2 cell types (macrophages and foam cells) (Franzén et al., 2016) from 855 human subjects. Not all tissues/cell types were collected from all subjects. The rationale for collecting these tissues and their role in coronary artery disease have previously been described (Björkegren et al., 2015). Ethical permits, biopsy and experimental procedures have previously been described (Franzén et al., 2016) (Karolinska Institutet Dnr 154/7 and 188/M-12). The samples from aortic wall (n = 538) and internal mammary artery (n = 552) were exclusively prepared from rRNA-depleted libraries generated with the Illumina RiboZero protocol. The remaining samples were almost all poly(A)-enriched. Of the 4301 samples included in this study, 52% (2,267/4,301) were sequenced with a strand-specific protocol, and the remaining were sequenced with a non-strand-specific protocol. A complete summary of tissues, samples and protocols is shown in Tables S1 and S2. Editing calls are listed in Dataset S1.

Metabolite and protein data were generated by Olink AB (Sweden) and Nightingale (formerly Brainshake, Finland).

Genome version, annotations and scripts

The human genome GRCh38 and GENCODE (Harrow et al., 2012) v.24 annotations were used throughout this study. RNA editing sites were mapped to genomic features with ANNOVAR (Wang, Li & Hakonarson, 2010) v.2015-03-22. miRBase (Kozomara & Griffiths-Jones, 2014) v.21 was used for miRNA annotations. Scripts used to process data can be retrieved from: https://github.com/oscar-franzen/rnaed.

Computational scheme for detection of RNA editing sites

Sequencing reads were aligned with the human genome using STAR (Dobin et al., 2013) v.2.5.1b, with the maximum number of alignments set to 1, and the ratio of mismatches to sequence length set to 0.1. Only uniquely mapped reads were considered (mapping quality = 255). The strand specificity of sequencing was confirmed by counting reads mapped to sense and antisense transcripts with HTSeq (Anders, Pyl & Huber, 2014) v.0.6.0. For all samples, the predicted strandedness matched the recorded information from the sequencing facility. Only samples with more than 1,000,000 uniquely mapped reads in genes were analyzed. The first step in our pipeline was collapsing polymerase chain reaction (PCR) duplicates with the rmdup command in samtools v.1.2. This step limits the number of false positives from the library construction process (i.e., amplification of errors introduced by the reverse transcriptase). Alignments were then scanned for single nucleotide variants (SNVs) by parsing the output of the mpileup command in samtools. SNVs were removed if located within simple repeats, homopolymers ≥ 5 bp, within 5 bp of splice junctions. SNVs located on the mitochondrial genome, on unplaced contigs, and within the major histocompatibility complex were also removed. Intersections were called with bedtools (Quinlan & Hall, 2010) v.2.21.0. Furthermore, SNVs listed in any of the following databases were discarded: NCBI dbSNP b141, b146, and b147; Exome Aggregation Consortium variants v.0.3.1; NHLBI ESP6500 variants; Scripps’ Wellderly (Erikson et al., 2016); and COSMIC (Forbes et al., 2015) v.77. For the remaining SNVs, all sequencing reads overlapping the site were extracted and realigned with GSNAP (Wu & Nacu, 2010) v.2016-05-25, using settings that allow mismatches and splice junctions. For each sequencing read, the alignment coordinates (chromosome, strand, start and stop positions) from STAR and GSNAP were compared, and discordant alignments were removed. A concordant alignment was defined as >99% overlap of the start-stop interval on the same chromosome and strand. To call an RNA editing event, the base quality score had to be ≥20. To avoid errors introduced from random hexamer priming, the SNV had to be >6 bases from the 5′ start position of the sequencing read. SNVs indicating more than one type of variant per site were ignored. The final criteria required at least two non-identical sequencing reads to support the RNA editing event.

121.2 billion single-end reads were mapped to unique locations in the human genome (median  = 27,539,768 reads/sample). Collapsing PCR duplicates reduced the numbers to ∼42.9 billion aligned reads (median  = 9,527,872 reads/sample). In total, 22,124,713 RNA editing events were called, mapping to 2,631,392 unique sites. 83.3% (2,866,471/3,438,790) and 92.6% (17,304,283/18,685,923) of the identified events were A-to-G mismatches in strand-specific and non-strand-specific libraries, respectively. In the non-strand-specific data, these numbers also include the reverse complement of A-to-G (i.e., T-to-C).

Quality of detected editing events

Sequencing errors are more likely to occur toward the end of reads. We did not find 5′ or 3′ positional biases in sequencing reads for canonical events (A-to-I and C-to-U; Fig. S13). Next, assuming that G-to-A events in strand-specific samples were sequencing or mapping errors, we estimated that the false discovery rate (FDR) was 3.1% (89,490/2,866,471). For non-strand-specific samples (assuming all G-to-C events were sequencing or mapping errors), the FDR was 0.83% (72,500/8,717,781).

Collected RNA samples from macrophages and foam cells are near biological replicates (foam cells are derived from macrophages by incubating them with acetylated LDL for 48 h). Therefore, we assessed the reproducibility of RNA editing in 235 macrophage—foam cell pairs. We first assessed the percentage of editing sites detected in both samples in each pair; the median across all pairs was 24%. It may therefore be plausible that many RNA editing sites are not strictly reproducible even in samples that are close to biological replicates. This may reflect the stochastic nature of RNA editing or variability in sequencing depth. We also determined whether the editing ratios were stable for the ∼24% of sites found in both samples of each pair. Across each of the 235 macrophage-foam cell pairs, the median Spearman correlation was 0.83 (Fig. S14). Thus, despite known technical and biological confounders, the reproducibility was relatively high for sites that can be detected in both samples.

Gene expression

Gene expression was calculated as RPKM (Mortazavi et al., 2008) from reads counted with HTSeq (Anders, Pyl & Huber, 2014).

Association-testing of editing ratios and clinical parameters

Association-testing was conducted with the linear model in MatrixEQTL (Shabalin, 2012) v.2.1.1. Clinical parameters were first rank-normal transformed with the function rntransform in the R package GenABEL (Aulchenko et al., 2007) with ties broken randomly. Two biological covariates (age and sex) and two technical covariates (sequencing batch and laboratory) were included in the statistical model. Only sites with non-zero A-to-G(I) editing ratios in more than 50 samples per tissue were considered. Significant associations were defined as those satisfying the following criteria: (i) adjusted P < 0.05 (10,000 permutations were executed for each clinical parameter-tissue combination); (ii) R-sq > 0.4; and (iii) no correlation between gene expression (of the overlapping gene) and the tested clinical parameter (nominal P > 0.05).

RNA-editing quantitative trait loci (edQTLs)

Genotypic data were processed as described (Franzén et al., 2016); biallelic markers (MAF > 0.05) were converted to allele dosages (0, 1, and 2). edQTLs were generated separately for each tissue. Only A-to-G(I) sites found (i.e., editing ratio > 0) in more than 40 samples were included in the analysis (sites on sex chromosomes were excluded). Site-sample combinations with fewer than <5 reads supporting editing were removed. Editing ratios were then rank-normal transformed as described, and the null hypothesis of no association between allele dosages and editing ratios was tested with a linear regression model in MatrixEQTL (Shabalin, 2012). For each hypothesis, we required at least 20 subjects in at least 2 genotype-allele groups. Thus, the minimum number of individuals for any test was 40. All possible genotype-editing site combinations were tested regardless of the absolute genetic distance. We corrected for the following covariates: age, sex, laboratory, batch, and population stratification (i.e., cryptic relatedness of study subjects; represented as four genetic multi-dimensional scaling components). The FDR was controlled using a permutation procedure (Franzén et al., 2016): the RNA editing matrix was scrambled 1,000 times, i.e., breaking the link between sample identifiers and corresponding measurements. For each editing site, we generated 1,000 permutation p-values, which were used to adjust the real p-values. An adjusted p-value <0.20 was considered significant.

Overlap with eQTLs

Normalized gene expression data from (Franzén et al., 2016) were used in a linear regression model. P-values were adjusted with the Benjamini–Hochberg procedure and adjusted P < 0.05 was considered significant. For each editing site, all overlapping genes were tested.

GWAS SNPs

The NHGRI-EBI GWAS catalog (Welter et al., 2014) (r2018-01-01) and GWASdb (Li et al., 2016) v.2 were downloaded and merged into one database, keeping only SNPs with P < 1e − 08. GWAS SNPs that were also edSNPs were saved in a new list that was LD-pruned with PLINK (–indep-pairwise 400 kb 1 0.6) (Purcell et al., 2007).

General characteristics of the RNA editome

To characterize the human RNA “editome”, we analyzed RNA-seq data from the STARNET (Franzén et al., 2016) study, in which samples were obtained from up to 9 tissues (atherosclerotic aortic wall, n = 538; non-atherosclerotic arterial wall (internal mammary artery), n = 552; liver, n = 545; skeletal muscle, n = 533; visceral abdominal fat, n = 533; subcutaneous fat, n = 532; whole blood, n = 559; macrophages, n = 256; and foam cells, n = 235). In addition, 18 coronary artery samples were included. Sequencing was done by using 50-bp single-end reads (these libraries were strand-specific; n = 2, 267) and 100-bp single-end reads (these libraries were non-strand-specific; n = 2,034). Differences in read-length and strand specificity affected the sensitivity to detect RNA editing, and most analyses were therefore done separately on these two sets of samples. See Table S1 for a summary of analyzed samples. We applied a rigorous computational pipeline that includes internal and external filters to identify RNA editing in a total of 4,301 RNA-seq samples (Fig. 1A; Fig. S1). The number of samples per tissue ranged from 235 in foam cells to 559 in whole blood. For each putative editing site, we calculated the RNA editing ratio, defined as the count of edited bases divided by the total number of bases covering the site. The RNA editing ratio takes values within the semi-closed interval (0,1].

Most of the identified events were A-to-G mismatches (Figs. 1B and 1C; Table 1; Tables S2 and S3), consistent with A-to-I editing; henceforth we refer to these events as A-to-G(I) editing (including the reverse complement when the strand information was not preserved). We identified 1,484,357 and 450,002 A-to-G(I) events in non-strand-specific and strand-specific libraries, respectively; the number of non-redundant A-to-G(I) editing events in non-strand-specific and strand-specific libraries was 1,688,815. More A-to-G(I) events were detected in non-strand-specific libraries, owing to deeper sequencing and longer read length. Liver had the highest number of A-to-G(I) events, possibly reflecting higher transcriptional complexity. In total, 53.9% (910,943/1,688,815) of identified A-to-G(I) events had previously been reported in REDIportal (Picardi et al., 2017) or DARNED (Kiran & Baranov, 2010) (Fig. S2). Aortic wall had the highest number of novel A-to-G(I) sites (n = 117, 493). A-to-G(I) editing was detected in transcripts from 25,951 genes, of which 62.1% (16,131/25,951) were protein-coding. We intersected RNA editing events with genome annotations (Fig. S3), showing that 56.0% of A-to-G(I) events were in introns, 18.3% in intergenic regions, and 9.4% in 3′ UTRs; 67.7% of the A-to-G(I) editing sites were tissue-specific when comparing exact chromosome-positions without considering differences in gene expression (Fig. S4). To exclude bias from differences in gene expression we repeated the analysis and only examined genes being robustly expressed in 7 tissues (defined as median reads per kilobase of transcript per million mapped reads [RPKM] >10 across all samples; 2,639 genes), the result showed that most events were tissue-specific (Fig. S4).

In strand-specific samples, the second most common change was noncanonical T-to-C editing (16.2%, 116,318/714,372). Interestingly, in strand-specific libraries the number of A-to-G(I) events correlated with the number of T-to-C events (Spearman’s ρ = 0.67, P < 2.2e−16). T-to-C events have been reported in studies utilizing high-throughput sequencing data but it is unknown whether T-to-C events are authentic or artifacts (Bahn et al., 2012).

The relationship between sequencing depth and number of editing events was nearly linear (Fig. 1D and Fig. S5; liver samples: Spearman’s ρ = 0.86, P < 2.2e−16), suggesting that detection of RNA editing per sample was not saturated or that the mechanisms mediating RNA editing are nonspecific. If RNA editing is nonspecific, it can be assumed that deeper sequencing does not saturate detection. Analysis of the cumulative number of unique A-to-G(I) sites as a function of added samples did not reveal a clear tendency toward saturation (Fig. S6A). It is also possible that saturation cannot be detected because most editing takes place in introns, which are much less captured with RNA-seq. However, analysis of the number of genes in which RNA editing occurred did reveal a trend toward saturation (Fig. S6B).

Differential editing in Alu subfamilies

In the present version of the human genome assembly (GRCh38) there are 1,186,920 Alus covering 10.2% (315,412,481 bp/3,088,269,832 bp). Alu was the sole element significantly enriched in editing events (Fig. 1E; Figs. S7–S9). Overall, 70.3% (1,187,358/1,688,815) of all A-to-G(I) events were in Alus. Of the Alus in the human genome, 20.3% (241,226/1,186,920) had at least one A-to-G(I) editing event. Despite covering 20.7% (641,953,033 bp/3,088,269,832 bp) of the genome (Lander et al., 2001), LINE elements contained only 4.0% (67,810/1,688,815) of the A-to-G(I) events. The mean editing ratio per sample was slightly lower in non-Alu A-to-G(I) compared to editing sites in Alus (mean_non−Alu = 0.54, mean_Alu = 0.57, P = 3.5e−15), likely reflecting the strong preference of ADAR for double-stranded folds. When Alu elements were grouped by subfamilies (Y, J, and S) according to their established phylogeny (Bennett et al., 2008), it showed that 42.5% (718,653/1,688,815) of A-to-G(I) editing events were located within the AluS subfamily, followed by AluJ (23.2%) and AluY (4.5%), see Fig. 1F. The mean editing ratio per sample was highest in AluY (mean_AluS = 0.56, mean_AluJ = 0.56, mean_AluY = 0.64, P < 2.2e−16; Fig. 1G). Mean editing levels were also significantly higher within nongenic Alus compared to genic Alus (mean_genicAlus = 0.56, mean_nongenicAlus = 0.66, P < 2.2e−16), perhaps because purifying selection lowered editing levels in mRNAs.

ADAR1 is ubiquitously expressed and ADAR2 displays a vascular-specific signature

ADAR1 and ADAR2 catalyze A-to-I editing (Macbeth et al., 2005; Nishikura, 2016). ADAR3 is catalytically inactive and may have regulatory properties (Chen et al., 2000). The expression of ADAR1, ADAR2, and ADAR3 was compared across the tissues/cell types (Fig. 2; Fig. S10). ADAR1 was uniformly expressed in all tissues (median = 31 RPKM, median absolute deviation [MAD] = 14); the lowest levels observed were in skeletal muscle (median = 5 RPKM, MAD =1 ), see Fig. 2A. The latter was also observed in GTEx (Melé et al., 2015). Expression of ADAR1 correlated with A-to-G(I) editing events in every tissue except foam cells (q-value < 0.05), further corroborating the authenticity of called RNA editing events. Whole blood showed the strongest correlation between A-to-G(I) editing events and ADAR1 expression (strand-specific samples, Spearman’s ρ = 0.85, P < 2.2e−16; Fig. S11). Although sex-specific differences in ADAR1 expression have been reported in tumor samples (Paz-Yaacov et al., 2015), we found only modest sex-specific differences in foam cells (P = 0.015), liver (P = 0.049), and subcutaneous fat (P = 0.017; Fig. S12).

ADAR2 expression was low in all tissues (Fig. 2B; median = 3 RPKM, MAD = 2) except the three vascular tissues (atherosclerotic aortic wall, internal mammary artery, and coronary artery; median = 23 RPKM, MAD = 14); expression was highest in internal mammary artery (median = 42 RPKM, MAD = 12). ADAR2 expression correlated with the number of A-to-G(I) events only in atherosclerotic aortic wall (Spearman’s ρ = 0.18, P = 2.6e−05) and was stronger in females (females: ρ = 0.19, P = 0.012; males: ρ = 0.17, P = 5.2e−04). However, there were no sex-specific differences in editing in any tissues (P > 0.05). ADAR3 was essentially below the detection level in all tissues (median = 0.01 RPKM, MAD = 0.02).

Eleven novel recoding sites and one loss of stop codon

RNA editing provides a mechanism for cells to generate additional protein diversity (Nishikura, 2010). However, novel protein recoding events—those that cause amino acid changes—are rare. We systematically searched for A-to-G(I) editing events causing recoding. We found that 3.9% (67,429/1,688,815) of A-to-G(I) events were within known protein-coding exons (variants arising from antisense transcription were excluded), and only 3.0% (50,883/1,688,815) were non-synonymous substitutions. The latter appeared in a total of 9,739 protein-coding genes. The median number of non-synonymous substitutions per sample was 18 (MAD = 19).

Editing ratios had a long-tailed skewed distribution (median = 0.14), suggesting that these variants are present in a small number of transcripts. Most of the non-synonymous sites (90.8%, 46,245/50,883) were found in single samples. Thus, although variants are pervasive on a global scale, many reflect background noise.

To focus on likely functional sites, we looked for conserved non-synonymous substitutions, that were present in >20 tissue samples and had a median editing ratio >0.1 in any tissue. Using these criteria, we found 65 conserved editing sites in 44 genes (Table S4). Six sites were predicted to be deleterious according to PROVEAN (Choi et al., 2012) and SIFT (Ng & Henikoff, 2001) (chr4:57,110,146, chr6:44,152,612, chr9:33,271,197, chr16:5,044,722, chr16:57,683,958 and chr17:1,534,142). As expected, conserved editing sites had higher editing ratios (median_conserved = 0.30, median_notconserved = 0.05, P < 2.2e−16, Mann–Whitney U test). Moreover, 80% (52/65) of these sites had been reported in REDIportal (Picardi et al., 2017), and 35% (23/65) had been reported in the literature (Table S4). The classical example of recoding in GRIA2 (glutamate ionotropic receptor AMPA type subunit 2) was found at the known site (chr4:157,336,723; p.Q607R) (Wright & Vissel, 2012) in aortic wall and internal mammary artery. The median editing ratio in each tissue was 1, indicating that all transcripts carried the modified base.

Eleven recoding sites were novel (Table 2). Seven had low editing ratios (<0.5; median = 0.24), which in part may explain why they were not detected previously. This finding highlights the power of large RNA-seq datasets such as STARNET.

Two recoding sites (one previously reported Zhang et al., 2014a and one novel) resided in IGFBP7 (insulin-like growth factor-binding protein 7). The first site, chr4:57,110,062 changes a codon encoding lysine to arginine (amino acid 97); it was found specific to the internal mammary artery and had been reported in brain (Zhang et al., 2014a). The second, chr4:57,110,146, changes a codon encoding glutamic acid to glycine (amino acid 69); it was found in both atherosclerotic aortic wall and the internal mammary artery. Both recoding sites were in the IGF-binding domain of IGFBP7, possibly highlighting a novel regulatory mechanism.

Two recoding sites (chr2:219,483,602 in SPEG and chr3:179,375,226 in MFN1) were previously reported only in mouse (Danecek et al., 2012), suggesting these mRNA-recoding sites are of ancient origin since they were present in the common ancestor of extant human and mouse lineages (diverged ∼65 to 110 million years ago (Emes et al., 2003)).

Next, we investigated A-to-I recoding events resulting in the loss of stop codons (the opposite, the gaining of stop codons, is not biologically possible by A-to-I editing). Using the same criteria described above, we found one recoding site, at chr10:101,017,585 in PDZD7 (PDZ domain containing 7). It was specific to visceral abdominal fat and modifies the stop codon UAG to UGG (encodes tryptophan), extending PDZD7 by 18 amino acids (WSQIPGLKLSTRLSLPKC). This recoding site was recently reported in human brain tissue (Hwang et al., 2016). To our knowledge this is the first and perhaps only recoding site in humans that causes protein extension. We checked if the change was associated with measured traits (e.g., LDL-C) and we did not find any association.

pri-mir-27b and mir-605 editing in aortic wall and internal mammary artery are associated with increased plasma levels of high-density lipoprotein

In preparing the sequencing libraries, we used the Ribo-Zero protocol to eliminate ribosomal RNA, which enabled us to investigate small non-coding RNAs in 996 RNA-seq samples from atherosclerotic aortic wall and internal mammary artery. We identified A-to-G(I) editing sites in 102 unique microRNAs (miRNAs; annotations span the entire predicted stem-loop region) and 73 small nucleolar RNAs (snoRNAs). Stringent editing criteria (recurrence in > 50 samples) identified 18 editing sites in 10 miRNAs (mir-24-2, mir-27a, mir-27b, mir-605, mir-641, mir-1254, mir-1285-1, mir-1304, mir-3126, and mir-3138; Table S5) and 2 snoRNAs (U88 and SNORD17). RNA editing sites in mir-24-2, mir-605, mir-3126, and U88 had not previously been reported in REDIportal (Picardi et al., 2017) or DARNED (Kiran & Baranov, 2010). Six editing sites occurred in mature miRNAs: three in mir-605-3p and one each in mir-3138, mir-605-5p, and mir-24-2-3p. Two editing sites were associated with higher HDL cholesterol levels: chr9:95,085,457, which was located outside the mature miRNA (pri-mir-27b; P = 2.2e−05; FDR = 8.5%) and chr10:51,299,636 (mir-605-3p; P = 1.6e−06; FDR = 0.65%).

RNA editing in 3′ UTRs of several genes is associated with plasma lipid levels

Next, we looked for associations between 18,458 RNA editing sites and 208 clinical parameters (Franzén et al., 2016). Thirty-nine associations were found significant (FDR < 5%; 10,000 permutations) and involved 25 genes in five tissues (Table S6). 52% (13/25) of these were located in 3′ UTRs, implying they could affect transcript regulation. Of the remaining sites, eight were in introns, two were in known regulatory elements, one was exonic, and three were intergenic. A subset of all significant associations is listed in Table 3.

Although none of the 25 genes had been directly associated with coronary artery disease, 22 of 25 sites were associated with plasma lipoprotein levels. Two of the affected genes may have roles in lipid metabolism: LRP11 (low density lipoprotein receptor-related protein 11) and PLIN5 (lipid storage droplet protein 5). In subcutaneous fat, the LRP11 editing site (chr6:149,822,432), located in the most 3′ intron of the longer isoform, was associated with plasma low density lipoprotein (LDL) levels. RNA editing of two genes known for their roles in innate immune response—C3 (complement component 3) and C1RL (complement C1r subcomponent like)—occurred in liver and were both associated with plasma HDL levels.

QTL mapping of RNA editing identifies novel candidate disease genes

To uncover possible genetic regulation of RNA editing (Kurmangaliyev, Ali & Nuzhdin, 2015), we looked for RNA editing quantitative trait loci (edQTLs; Fig. 3A). Using data from the 1000 Genome Project (Durbin et al., 2010), imputation resulted in 6,246,842 common variants with minor allele frequency >5%. Next, we identified 17,718 edQTLs at FDR <20% (i.e., associations between allele dosages and RNA editing ratios; Fig. 3B) by performing ∼2.6 billion hypothesis tests (comprising 11,280 unique editing sites). We examined whether identified edQTLs were also expression QTLs (eQTLs). 19.5% (3,465/17,718) of RNA editing sites in edQTLs were eQTLs with respect to the overlapping gene. Across all tissues, we found 890 index edQTLs—the most significant association for each RNA editing site (Dataset S2). A total of 77% (690/890) of the index edQTLs were trans (the regulatory SNP and the editing site are located on different chromosomes).

Most of the editing sites with edQTLs were enriched in Alu elements (87%, 783/890); the Alu-subtype distribution was the same as shown in Fig. 1F. A total of 95% (853/890) of the edQTLs have previously reported editing sites (Kiran & Baranov, 2010; Picardi et al., 2017); 47% (425/890) were in 3′ UTRs, 25% (230/890) were in introns, and 12% (113/890) were in intergenic regions. RNA editing sites of edQTLs overlapped 417 genes (Dataset S2).

To assess a possible mechanistic role of edSNPs in disease, we intersected edSNPs with lead SNPs from GWAS. 90 edSNPs coincided with GWAS SNPs, encompassing 59 traits. LD-pruning collapsed the 90 SNPs to 30 non-linked SNPs. rs2573346 is associated with sarcoidosis (Hofmann et al., 2008) and had an edQTL in the 3′ UTR of FAM213A in visceral fat (Fig. 3C). rs2028299 is associated with type 2 diabetes (Kooner et al., 2011) and had two edQTLs in visceral fat (both located in the 3′ UTR of ARPIN). rs2028299 is associated with expression of ARPIN in fat tissue (P < 1e−10) (Kooner et al., 2011). A negative correlation was found between ARPIN expression in visceral fat and plasma levels of VLDL (r =  − 0.2; P < 1e−06). rs78579285 has been associated with joint mobility (Beighton score) (Pickrell et al., 2016) and we detected an edQTL in visceral fat within the intron of PIEZO1 (encodes a mechanosensitive ion-channel component). rs7546668 is associated with kidney function (glomerular filtration rate) (Gorski et al., 2017) and it formed an edQTL with an editing site detected in liver (Fig. 3D). The editing site overlapped DNAJC16 (encodes DnaJ Heat Shock Protein Family (Hsp40) Member C16) and AGMAT (encodes Agmatinase). rs4073054 is associated with HDL (Chasman et al., 2009) and had an edQTL with an RNA editing site in the 3′ UTR of F11R in liver. Expression of F11R was associated with plasma levels of VLDL (r =  − 0.2, P < 1e−08) and HDL (r =  − 0.2, P < 1e−07). rs1127311 is a pleiotropic locus and it is associated with C-reactive protein and triglyceride levels (Ligthart, Vaez & Hsu, 2016). Interestingly, rs1127311 is situated in the 3′ UTR of ADAR, and it forms an edQTL with an editing site in the second to the last intron of FGG. rs3947, situated in the 3′ UTR of the protease-encoding gene CTSB, is associated with blood protein levels (Suhre et al., 2017). rs3947 forms an edQTL with an editing site located in the 3′ UTR of CTSB. rs3947 is in strong LD with rs2740594 (R-sq = 0.86), which is associated with Parkinson’s disease (Chang et al., 2017). rs4957048 is associated with multiple sclerosis/ulcerative colitis (De Lange et al., 2017), and it was associated with an editing site located downstream of SLC22A25 (Fig. 3E). rs6756513 is associated with platelet count (Astle et al., 2016) and it had an edQTL in an intron of C3P1. rs10847434 is associated with coronary artery disease (Lee et al., 2013) and had an edQTL with an editing site in the most 3′ exon of APOC1P1 in liver. APOC1P1 encodes a long non-coding RNA without any known function or role. However, its locus (Jeemon et al., 2011) has been extensively linked to coronary artery disease (neighboring genes are APOE, APOC1, and TOMM40). In liver, there was a negative correlation between expression of APOC1P1 and APOE (r =  − 0.3; P = 7e−13) as well as APOC1P1 and APOC1 (r =  − 0.22; P = 1e−07). The negative correlation may indicate a regulatory function of APOC1P1. rs1883350 is associated with fatty liver disease (Feitosa et al., 2013) and we identified one liver edQTL slightly downstream of PNPLA3 (Patatin-like phospholipase domain-containing protein 3). PNPLA3 is associated with fatty liver disease (Speliotes et al., 2010) and encodes a protein known to be involved in liver fat storage. Finally, rs4739066 is associated with myocardial infarction (Wakil et al., 2016) and we found this SNP to have two edQTLs with editing sites in the 3′ UTR of the α-tocopherol transfer protein gene (TTPA). The editing sites are located 64 bp apart and are intriguing since they have opposite effect sizes. TTPA is associated with the severity of atherosclerotic lesions in the proximal aorta (Terasawa et al., 2000).