# Single Nucleotide Polymorphisms with Cis-Regulatory Effects on Long Non-Coding Transcripts in Human Primary Monocytes

Jonas Carlsson Almlöf, Per Lundmark, Anders Lundmark, Bing Ge, Tomi Pastinen, Alison H Goodall, François Cambien, Panos Deloukas, Willem H Ouwehand, Ann-Christine Syvänen

<figure><img src="image_1.png" /></figure>

To cite this version:

Jonas Carlsson Almlöf, Per Lundmark, Anders Lundmark, Bing Ge, Tomi Pastinen, et al. Single Nucleotide Polymorphisms with Cis-Regulatory Effects on Long Non-Coding Transcripts in Human Primary Monocytes. PLoS ONE, 2014, 9 (7), pp.e102612. 10.1371/journal.pone.0102612. hal-01344523

HAL Id: hal-01344523

https://hal.sorbonne-universite.fr/hal-01344523v1

Submitted on 12 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Single Nucleotide Polymorphisms with Cis-Regulatory Effects on Long Non-Coding Transcripts in Human Primary Monocytes

Jonas Carlsson Almlöf¹, Per Lundmark¹, Anders Lundmark¹, Bing Ge², Tomi Pastinen², Cardiogenics Consortium¶, Alison H. Goodall³,⁴, François Cambien⁵, Panos Deloukas⁶,⁷,⁸, Willem H. Ouwehand⁷,⁹, Ann-Christine Syvänen¹*

¹Department of Medical Sciences, Molecular Medicine and Science for Life Laboratory, Uppsala University, Uppsala, Sweden, ²Department of Human Genetics, McGill University, Montréal, Canada, ³Department of Cardiovascular Science, University of Leicester, Leicester, United Kingdom, ⁴Leicester NIHR Biomedical Research Unit in Cardiovascular Disease, Glenfield Hospital, Leicester, United Kingdom, ⁵INSERM UMRS 937, Pierre and Marie Curie University and Medical School, Paris, France, ⁶William Harvey Research Institute, Barts and The London School of Medicine and Dentistry, Queen Mary University of London, Heart Centre, Charterhouse Square London, United Kingdom, ⁷Wellcome Trust Sanger Institute, Hinxton, Cambridge, United Kingdom, ⁸Princess Al-Jawhara Al-Brahim Centre of Excellence in Research of Hereditary Disorders (PACER-HD), King Abdulaziz University, Jeddah, Saudi Arabia, ⁹National Health Service Blood and Transplant, Cambridge Centre, Long Road, Cambridge, United Kingdom

<figure><img src="image_3.png" /></figure>

## Abstract

We applied genome-wide allele-specific expression analysis of monocytes from 188 samples. Monocytes were purified from white blood cells of healthy blood donors to detect cis-acting genetic variation that regulates the expression of long non-coding RNAs. We analysed 8929 regions harboring genes for potential long non-coding RNA that were retrieved from data from the ENCODE project. Of these regions, 60% were annotated as intergenic, which implies that they do not overlap with protein-coding genes. Focusing on the intergenic regions, and using stringent analysis of the allele-specific expression data, we detected robust cis-regulatory SNPs in 258 out of 489 informative intergenic regions included in the analysis. The cis-regulatory SNPs that were significantly associated with allele-specific expression of long non-coding RNAs were enriched to enhancer regions marked for active or bivalent, poised chromatin by histone modifications. Out of the lncRNA regions regulated by cis-acting regulatory SNPs, 20% (n = 52) were co-regulated with the closest protein coding gene. We compared the identified cis-regulatory SNPs with those in the catalog of SNPs identified by genome-wide association studies of human diseases and traits. This comparison identified 32 SNPs in loci from genome-wide association studies that displayed a strong association signal with allele-specific expression of non-coding RNAs in monocytes, with p-values ranging from 6.7×10⁻⁷ to 9.5×10⁻⁸⁹. The identified cis-regulatory SNPs are associated with diseases of the immune system, like multiple sclerosis and rheumatoid arthritis.

**Citation:** Carlsson Almlöf J, Lundmark P, Lundmark A, Ge B, Pastinen T, et al. (2014) Single Nucleotide Polymorphisms with Cis-Regulatory Effects on Long Non-Coding Transcripts in Human Primary Monocytes. PLoS ONE 9(7): e102612. doi:10.1371/journal.pone.0102612

**Editor:** Tanja Zeller, Medical University Hamburg, University Heart Center, Germany

**Received** February 21, 2014; **Accepted** June 20, 2014; **Published** July 15, 2014

**Copyright:** © 2014 Carlsson Almlöf et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Funding:** The study was performed as part of WP5 of the Cardiogenics EU FP6 project (grant number LSHM-CT-2006-037593). The work leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreements no 262055, European Sequencing and Genotyping Infrastructure (ESGI) and 282510 BLUEPRINT. A full list of the investigators in the Blueprint consortium who contributed to the generation of the data is available from www blueprint-epigenome.eu. The study received funding from the Swedish Research Council (grants no A028001 and C0524801 to A-CS) and the Knut and Alice Wallenberg Foundation (KAW 2011.0073 to A-CS). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing Interests:** The authors have declared that no competing interests exist.

## Introduction

Determination of allele-specific gene expression (ASE) levels by quantitative genotyping of heterozygous SNPs on the RNA level [1], with genome-wide panels of SNPs [2,3] can be used as a guide for identifying cis-acting genetic variants that regulate gene expression. ASE analysis is more powerful for detecting cis-regulated gene expression than the total expression levels of genes determined by regular eQTL analysis [4]. The power and precision of ASE analysis to detect cis-regulatory SNPs (cis-rSNPs) stems from the fact that the differential expression of the two alleles of a SNP are measured in the same environment and have

been exposed to the same environmental conditions in the cells from which the RNA was extracted [4]. Thus trans-acting and environmental factors that may affect gene expression are controlled for in ASE analysis. To adjust for possible methodological differences in the efficiency of the genotyping assay to measure the two alleles of a SNP, such as the sequence context of a SNP or copy number variations, the relative amounts of the alleles measured in DNA, is used as a reference for quantification of the allelic expression. Genome-wide ASE analysis by SNP genotyping has been applied to map cis-regulation of protein coding genes associated with human diseases and traits in lymphoblastoid cell

<footnote>* Email: Ann-Christine.Syvänen@medsci.uu.se</footnote>

<footnote>¶ Members of the CARDIOGENICS Consortium are listed in the Acknowledgements.</footnote>
lines [2,5], osteoblasts [6], fibroblasts [5,7], T-cells [8], and
monocytes [5].

Long non-coding RNAs (lncRNAs) are involved in gene
regulation and other cellular processes. LncRNAs such as XIST
and TSIX, which are involved in X chromosome inactivation
[9,10], are well known. Other well known examples of lncRNAs
are AIR that suppresses gene expression via hypermethylation and
HOTAIR that interacts with the Polycomb Repressive Complex 2
to silence the HOXD locus [11]. New functions for lncRNAs are
continuously being discovered. LncRNAs can affect gene expres-
sion in many ways, as scaffolds or guides for chromatin
modifications, as decoys for reducing the amount of transcription
factors interacting with chromatin, as signaling molecules reflect-
ing active transcription factor complexes, or as reservoirs for
microRNAs [12]. The mechanism for regulation of gene
expression by lncRNAs usually involves formation of RNA-protein
complexes that influence the gene expression. Recent large studies
[13,14,15,16,17,18] have found that lncRNAs display positive
correlations with expression of protein-coding genes in *cis*, and
especially with genes that overlap with the antisense strand of a
lncRNA. By investigating the effect of regulatory SNPs on the
expression of lncRNAs using traditional microarray-based ex-
pressed quantitative trait locus (eQTL) analysis of peripheral blood
cells a recent study identified 112 *cis*-regulated lncRNAs [19].

The aim of our study was to apply ASE-analysis to identify cis-
acting SNPs that regulate the expression of lncRNAs and to
highlight previously unknown roles for lncRNAs in human
complex diseases. We reasoned that ASE-analysis that is highly
sensitive for detecting *cis*-rSNPs would be a powerful tool for
studying *cis*-regulation of lncRNAs that are expressed at lower
levels than protein-coding genes. To our knowledge ASE-analysis
has not been previously used for studies of lncRNAs.

We analysed RNA extracted from human monocytes purified
from white blood cells of 188 healthy blood donors using a
genome-wide panel of SNPs. Monocytes were selected for analysis
because they are a relevant cell type for multiple diseases. Our
analysis included 8929 genomic regions that harbor genes for
potential long non-coding RNA that were retrieved from the
GENCODE database [18]. Of these gene regions, 60% were
annotated as intergenic, and thus they do not overlap with protein-
coding genes. Using stringent criteria for the identification of SNPs
that regulate the expression of lncRNAs we identified 8267 cis-
acting regulatory SNPs out of which 3910 are located in intergenic
regions and 571 of these are located in enhancer regions marked
as active or poised by histone modifications in monocytes. We also
compared the indentified SNPs that regulate expression of
lncRNAs with risk SNPs for complex diseases and traits previously
identified in genome-wide association studies (GWAS). In this way
we were able to obtain new functional clues for these disease
associated GWAS loci. Additionally, we analysed co-expression
between lncRNA genes and nearby protein coding genes, and
found that 20% of the intergenic lncRNAs with a cis-rSNP were
co-regulated with a protein-coding gene.

**Materials and Methods**

Samples

Circulating monocytes were collected from healthy adult blood
donors of European origin (n = 188) recruited from the United
Kingdom National Blood Service Centre in Cambridge, UK as
part of the Cardiogenics Transcriptomic Study [20]. Volunteer
donors with a self-reported recent or acute illness were excluded.
Donors with measured full blood cell count and C-reactive plasma

protein levels outside the normal ranges were also excluded from
the study.

DNA was extracted from peripheral blood leukocytes using the guanidine hydrochloride - chloroform method. CD14+ magnetic microbeads (autoMACS Pro, Miltenyi Biotec, Bergisch Gladbach, Germany) were used to isolate monocytes from whole blood. RNA was extracted from cell pellets of freshly isolated monocytes by homogenization with Trizol-reagent (Invitrogen, Paisley, UK), chloroform-ethanol extraction and purification using Qiagen RNAeasy columns and reagents, followed by on-column DNase treatment. cDNA was synthesized using reagents from the Illumina TotalPrep RNA Amplification Kit, except that the poly-dT primers were substituted by random decamers (Applied Biosystems, Carlsbad, California, US).

Ethics statement

All participants in the study provided their written consent after
a personal meeting with a study nurse. Written consent from each
participant was recorded by name, date and signature on a
standardized consent form. Both the consent form and the study
has been reviewed and approved by the Cambridgeshire 1
Research Ethics Committee. The Research Ethics Committee is
independent from the research institute and part of the UK Health
Research Authority.

ENCODE lncRNA Regions

A total of 8929 genomic regions containing lncRNAs were
retrieved from the GENCODE version 7 database [18]. From
regions containing several transcripts, only the longest transcript
was included to avoid analyzing completely overlapping regions.
Of the lncRNA regions, 5346 were annotated as intergenic, which
implies that that there is no overlap between the lncRNA and any
known protein coding gene. The remaining 3583 regions overlap
with either exons, introns or encompass an entire gene on the
sense or anti- sense strand.

Allele-specific expression analysis

RNA (cDNA) and genomic DNA (gDNA) from the monocyte
samples were genotyped using the Infinium assay and Human
1.2 M Duo custom BeadChips (Illumina, San Diego, California,
USA) as described previously [4]. The genotype data from cDNA
is used as a quantitative measure for gene expression in the
calculation of the ASE-levels. The genotypes called in gDNA are
used for three different purposes i.e. (i) to phase the SNP alleles on
each chromosome to facilitate ASE-calling using multiple SNPs
per transcript, (ii) to correct the ASE-levels in cDNA for possible
allelic bias in gDNA due to sequence context of the SNPs, and (iii)
to test for associations between potential cis-rSNPs and the ASE-
levels of lncRNAs in a genomic region. Genotypes were called in
gDNA using Genome Studio version 2009.2 (Illumina) with a call
rate of 99% as the threshold for genotype calls for SNPs and 98%
for samples. SNPs were further filtered on deviations from Hardy-
Weinberg equilibrium with a p-value cutoff of 10<sup>-6</sup> (Chi-squared
test). One sample with higher than 40% call rate in cDNA were
removed due to possible DNA contamination. The raw two-colour
fluorescence signals from the Infinium assay were normalized to
remove dependency of the allele fractions on the signal intensities
of the fluorophores using a quadratic function with medians of
bins for the two fluorescence colors. From the regression, the
intensities were predicted based on the log10 values of the raw
signal, which were used to adjust the allele fraction, see Methods
S1 for further details and equations. ASE levels were calculated for
each heterozygote SNP as the difference in normalized allele
fractions between cDNA and gDNA: [Allele1<sub>cDNA</sub>/(Allele1<sub>cDNA</sub>+
Allele2<sub>cDNA</sub>] - [Allele1<sub>gDNA</sub>/(Allele1<sub>gDNA</sub>+Allele2<sub>gDNA</sub>)]. Here, gDNA, where the two SNP alleles are present in 1:1 ratio serves as a quantification standard for the relative expression levels of the two alleles in RNA.

A gene region without ASE will obtain an ASE-level close to 0, while a gene region with ASE, where the gDNA fraction is 0.5, will have an absolute ASE-level $0<X<=0.5$. The software IMPUTE 2 [21] (v2.1.0) was used to impute missing genotypes and to phase the genotypes across lncRNA gene regions in each individual sample. As a reference panel prefiltered haplotypes from HapMap3 release #2 [22] and 1000 Genomes pilot1 [23] available at the IMUPUTE 2 website were used. The ASE levels were assigned in each individual sample as the average ASE level for all heterozygous SNPs within a genomic region (ASE window), corresponding to the region of an annotated lncRNA. Windows of lncRNAs with less than three informative heterozygous SNPs were excluded. The association between SNPs and the ASE levels of lncRNAs was assessed by linear regression of the ASE levels in the groups of samples with heterozygous SNPs and homozygous SNPs [4]. In summary, homozygous SNPs (AA and BB) are in one group with an expected ASE-level close to zero, and heterozygous SNPs are either AB or BA with a stochastically determined direction during phasing. In the case of significant ASE the AB and BA group will on average obtain different signs of the ASE-level. SNPs with significant associations with allele specific expression are denoted as *cis*-rSNPs (Figure S1).

All analyses described above were performed and data is presented using the NCBI36/hg18 assembly as reference genome.

## Refseq protein coding regions

Protein coding regions extracted from the Refseq database were analysed to detect SNPs with association signals that overlap with those from lncRNAs. The total number of transcripts in Refseq (downloaded 12th of January 2013) was 42797. To retain only protein coding genes, all genes annotated as RNA genes in the GeneCards database [24] were removed, leaving 38387 transcripts. To avoid completely overlapping transcripts, 19654 transcripts that are a subsequence of another transcript were removed. As for lncRNAs, only regions that contained at least three informative SNPs were included in the analysis, leaving 10345 protein-coding regions for further analysis.

## Determination of total gene expression levels using the genotype data

The sum of the raw fluorescence signal intensities from both alleles of a SNP in cDNA was used as a measure of total gene expression in the position of the SNP. The average sum of the SNP signals across each lncRNA window and across all samples represents the total expression value for each lncRNA window. A stringent average signal threshold of 1000 fluorescence units for all 188 samples was set for the summed signal intensities, see Figure S2.

## Enhancer regions

Monocyte chromatin immunoprecipitation sequencing (ChIP-seq) data was retrieved from the Blueprint project [25]. We analysed the four samples available in the Blueprint data using two chromatin signatures for each sample; H3K27ac, corresponding to an active enhancer [26,27] and H3K4me1, corresponding to a poised enhancer [28,29]. The peaks supplied in the Blueprint data had been called by MACS2 using the standard parameters for each signature. For a ChIP-seq peak to be included in our analysis, we required a p-value of $10^{-5}$ in at least two out of four

overlapping peaks. To be counted as an overlapping peak the nucleotide overlap was required to be over 50%. This resulted in 22091 peaks for H3K27ac and 32692 peaks for H3K4me1.

## Co-expression of lncRNA and closest protein coding gene

We investigated to what extent expression of lncRNA genes with a *cis*-rSNP were correlated with gene expression of the closest protein coding gene. For this purpose we divided the samples into three groups based on three genotype combinations of the most significantly associated *cis*-rSNP (AA, AB/BA, BB). For each expressed gene region we calculated the average signal intensity for both alleles and used this value as a measure of gene expression levels. Next, we performed linear regression analysis with the three genotype groups as x-values and the average total expression levels as y-values in order to obtain a p-value for the co-expression.

## Results

### Allele-specific expression of lncRNA regions

We explored *cis*-regulation of 8929 genomic regions harboring lncRNA from the ENCODE project using allele-specific expression (ASE) analysis of 188 RNA samples from human primary monocytes. To determine ASE, we genotyped RNA (cDNA) and genomic DNA from the monocytes using a genome-wide panel of 1.2 million SNP markers [4]. To be informative for the detection of ASE, a SNP has to be heterozygous in DNA and expressed at a detectable level in RNA. Out of the 8929 lncRNA regions, 60% were annotated as intergenic. 1298 regions contain at least three informative SNPs and out of these 1122 regions (489 intergenic regions) showed a signal intensity above background and were considered to be expressed in the monocytes and were thus included in the ASE analysis. Of the individual SNPs in the lncRNA regions, two thirds have an intensity level above 1000 in at least 90% of the samples. We found that the mean expression levels of the lncRNAs were 1.5-fold lower than the expression levels of exons (Figure 1), while the lncRNAs were expressed at a 1.5-fold and 5.5-fold higher level than intronic and inter protein coding gene regions, respectively.

The definitions of exon and intron boundaries were taken from Refseq. As protein coding genes are on average expressed at higher levels than lncRNAs, the ASE signal for lncRNA that overlap with protein coding genes is hard to distinguish from that of a protein-coding gene. Because of this, we focus on the intergenic lncRNA regions where the ASE data is most reliable.

For the ASE analysis, the genotypes of 194530 SNPs located within 250 kb upstream and downstream of the lncRNA regions were tested for their association with the ASE levels of the corresponding expressed lncRNAs, using linear regression, see Manhattan plot in Figure 2. A p-value cut-off of $p<10^{-6}$, based on the number of tests performed and adjusted for SNPs in high LD ($>0.9$), was used to correct the association signals for multiple testing. In addition a threshold of 0.05 for the slope of the regression line was applied to exclude SNPs with low effects on ASE. This value corresponds to approximately a 20% expression difference between the two alleles and was set as a reasonable threshold for biological relevance. To avoid false positives with inflated p-values, we only retained associations that were based on at least four data points for each allele combination for the calculation of the p-value [4] (see example regression plots in Figure S3). Using these criteria, we detected 258 intergenic (53%) lncRNAs with at least one associated *cis*-rSNP (Table S1). The associated intergenic lncRNA regions varied in size from 0.9 kb to 547 kb. A total of 8267 of the 194530 tested SNPs (4.2%) were
**Intergenic lncRNA**

<figure><img src="image_4.png" /></figure>

**Intronic**

<figure><img src="image_5.png" /></figure>

**Exonic**

<figure><img src="image_6.png" /></figure>

**Intergenic**

<figure><img src="image_7.png" /><figcaption>**Figure 1. Comparison of total expression levels of regions annotated to intergenic lncRNAs, exons, introns, and intergenic regions.**</figcaption></figure>

The horizontal axes in the panels show bins of fluorescence signals from the genotyping data, summed for both alleles to give a measure for total expression. The average expression levels of annotated transcripts were 4900 fluorescence units in exons, 2100 in introns, 590 in intergenic regions, compared to 3300 in the intergenic lncRNA regions that were used in the ASE analysis. The vertical axes show the number of observations in each bin.

doi:10.1371/journal.pone.0102612.g001

associated with ASE of a lncRNA region (p-value<10⁻⁶) in the 188 analysed monocyte samples, and for the intergenic lncRNA regions 3910 out of 79628 tested SNPs (4.9%) were significantly associated *cis*-rSNPs.

To allow comparison between ASE of lncRNA genes and protein coding Refseq genes, we applied the same procedure for ASE analysis to the Refseq gene regions as for the lncRNAs. This analysis identified 8604 *cis*-rSNP out of 108881 tested SNPs (7.9%) in 471 out of 10345 regions (4.6%) that were associated with ASE of a Refseq gene.

**Associated SNPs in enhancer regions**

We identified 571 *cis*-rSNPs located in enhancer regions marked as active or poised by histone modifications in monocytes. To determine if the *cis*-rSNPs for lncRNAs are enriched in enhancer regions, we calculated the fold enrichment of significantly associated *cis*-rSNPs defined by ASE analysis of intergenic lncRNA regions to histone marks for active enhancers (H3K27ac) and bivalent, poised enhancers (H3K4me1) in monocytes. The fold enrichment was determined by comparing the fraction of significantly associated SNPs in enhancer regions with the fraction of non-significant (p-value>0.5) SNPs in enhancer regions. Of the
<figure><img src="image_8.png" /><figcaption>**Figure 2. Manhattan plot.** Manhattan plot with the p-values from ASE association tests between SNPs and lncRNAs on the vertical axis and the genomic lncRNA regions analysed in the study on the horizontal axis. The p-value cut-off of $10^{-6}$ is shown as a grey line.
doi:10.1371/journal.pone.0102612.g002</figcaption></figure>

significantly associated SNPs, 2.1% were located in H3K27ac peaks compared to 1.2% for non-significantly associated SNPs, which corresponds to a 1.72-fold enrichment (p-value $9.0 \times 10^{-6}$, Fisher's exact test). For H3K4me1 we observed a 1.54-fold enrichment (p-value $3.1 \times 10^{-16}$, Fisher's exact test) of significantly associated SNPs, with 12.5% of the significant SNPs located in the enhancer peaks, while this number was 8.1% for non-significant SNPs. For protein coding genes the percentages and fold enrichment are similar to those for lncRNA for H3K27ac with 2.2% and 1.5% for significant and non-significant SNPs, respectively, and 1.43 fold enrichment (p-value $7.0 \times 10^{-6}$, Fisher's exact test) and less pronounced for H3K4me1 with 9.8% and 8.5% for significant and non-significant SNPs, respectively, and 1.15 fold enrichment (p-value $3.5 \times 10^{-4}$, Fisher's exact test). This result shows that enhancers have a regulatory cis-effect on expression of not only protein coding genes, but also lncRNAs.

## Co-expression of lncRNA and closest protein coding gene

Next we investigated whether the lncRNA genes and nearby protein coding genes are co-expressed and co-regulated by a cis-rSNP. In our data 52 (20%) of the identified intergenic lncRNA regions to which a cis-acting regulatory SNP was associated were also correlated with the expression level of a protein-coding gene

at a false discovery rate (FDR) of 10%, suggesting co-regulation of the genes, see Table 1. We found that 21 of the cis-rSNPs are significantly associated both with ASE of a lncRNA and ASE of the closest protein coding gene, which implies that the expression of the lncRNA and the protein coding genes are regulated by the same cis-rSNP. For 19 cis-rSNPs the ASE signal is unique to a lncRNA, which could be an indication of involvement of the lncRNA in the regulation of the expression of an adjacent protein coding gene as has been suggested in previous studies [14,16]. The remaining 12 cis-rSNPs had a distance >250 kb to the nearest protein coding gene and the ASE signal could therefore not be evaluated.

## GWAS lead SNPs overlap with cis-rSNP for lncRNAs

Next we superimposed the cis-rSNP associated with ASE of intergenic lncRNAs in monocytes with SNPs identified in GWAS of human diseases and traits, listed in the National Human Genome Research Institute (NHGRI) GWAS Catalog [30]. The GWAS catalog (downloaded 12th of January 2013) includes 9617 entries, 641 traits and 7797 unique SNPs. We identified SNPs associated with 32 loci high-lighted by GWAS that are also cis-rSNP for lncRNAs. The cis-rSNP are the same as the GWAS associated SNP for 25 of the loci, while seven cis-SNPs are linkage disequilibrium (LD) proxies to the lead GWAS SNPs, with five
<table><caption>**Table 1.** Co-expressed lncRNA regions and protein-coding genes.</caption><thead><tr><th>SNP<sup>1</sup></th><th>LncRNA region</th><th>Refseq gene</th><th>p-value for co-expression<sup>2</sup></th><th>p-value for ASE of lincRNA<sup>3</sup></th><th>p-value for ASE of Refseq gene<sup>4</sup></th></tr></thead><tbody><tr><td colspan="6"><strong>Co-regulated lncRNA and protein-coding by the same cis-rSNP</strong></td></tr><tr><td>rs644234</td><td>9:135120874-135140438</td><td>ABO</td><td>5.38E-26</td><td>2.69E-66</td><td>2.69E-66</td></tr><tr><td>rs1056787</td><td>4:67965618-68016564</td><td>CENPC1</td><td>6.75E-06</td><td>9.89E-25</td><td>1.37E-57</td></tr><tr><td>rs7313235</td><td>12:9980445-9987361</td><td>CLEC12A</td><td>1.92E-18</td><td>2.23E-37</td><td>5.00E-56</td></tr><tr><td>rs2271101</td><td>4:178603127-178835274</td><td>AGA</td><td>1.81E-12</td><td>4.79E-23</td><td>4.68E-54</td></tr><tr><td>rs1112956</td><td>5:127304017-127446691</td><td>SLC12A2</td><td>1.62E-10</td><td>2.56E-70</td><td>7.74E-49</td></tr><tr><td>rs10240848</td><td>7:38347741-38384763</td><td>AMPH</td><td>3.45E-05</td><td>1.80E-59</td><td>2.24E-45</td></tr><tr><td>rs2071904</td><td>22:48328333-48336125</td><td>C22orf34</td><td>1.92E-08</td><td>3.14E-23</td><td>4.02E-43</td></tr><tr><td>rs2254177</td><td>17:45713937-45720215</td><td>TMEM92</td><td>1.02E-12</td><td>1.99E-16</td><td>2.02E-41</td></tr><tr><td>rs1044303</td><td>6:131190239-131199966</td><td>SMLR1</td><td>1.46E-05</td><td>2.26E-30</td><td>2.26E-30</td></tr><tr><td>rs2949192</td><td>7:1166536-1171429</td><td>ZFAND2A</td><td>1.24E-03</td><td>1.70E-26</td><td>2.60E-30</td></tr><tr><td>rs9355652</td><td>6:158623283-158653378</td><td>TULP4</td><td>1.38E-02</td><td>2.12E-55</td><td>2.31E-28</td></tr><tr><td>rs8112960</td><td>19:21562075-21568648</td><td>ZNF429</td><td>1.85E-24</td><td>4.41E-27</td><td>7.21E-24</td></tr><tr><td>rs7647643</td><td>3:158363935-158369616</td><td>CCNL1</td><td>2.53E-06</td><td>2.91E-31</td><td>1.77E-22</td></tr><tr><td>rs7336525</td><td>13:20775651-20820860</td><td>MIPEPP3</td><td>1.91E-11</td><td>3.63E-25</td><td>5.02E-21</td></tr><tr><td>rs12366</td><td>2:75020893-75023305</td><td>POLE4</td><td>1.83E-04</td><td>5.44E-91</td><td>3.18E-19</td></tr><tr><td>rs672527</td><td>1:180685904-180796316</td><td>RGSL1</td><td>1.40E-07</td><td>1.06E-18</td><td>1.06E-18</td></tr><tr><td>rs17802159</td><td>4:114420-147779</td><td>ZNF876P</td><td>6.86E-06</td><td>3.48E-59</td><td>6.90E-15</td></tr><tr><td>rs10463951</td><td>5:135493095-135498478</td><td>SMAD5</td><td>7.96E-03</td><td>6.44E-14</td><td>4.58E-14</td></tr><tr><td>rs178255</td><td>22:19641380-19648967</td><td>AIFM3</td><td>1.04E-04</td><td>2.34E-20</td><td>5.08E-14</td></tr><tr><td>rs12711793</td><td>2:114453616-114481349</td><td>ACTR3</td><td>2.34E-03</td><td>1.50E-25</td><td>5.58E-14</td></tr><tr><td>rs9605146</td><td>22:15462934-15514699</td><td>TPTEP1</td><td>3.68E-16</td><td>5.12E-08</td><td>5.12E-08</td></tr><tr><td colspan="6"><strong>Regulation of lncRNA by cis-rSNP, but not of protein-coding gene</strong></td></tr><tr><td>rs246105</td><td>5:108600720-108689969</td><td>PJA2</td><td>1.20E-03</td><td>5.12E-24</td><td>2.83E-06</td></tr><tr><td>rs3862666</td><td>11:60579785-60591583</td><td>CD5</td><td>7.26E-03</td><td>4.08E-08</td><td>6.24E-06</td></tr><tr><td>rs948421</td><td>8:61459701-61591893</td><td>RAB2A</td><td>8.47E-03</td><td>6.08E-30</td><td>9.84E-06</td></tr><tr><td>rs3737813</td><td>1:178162399-178177554</td><td>TOR1AIP2</td><td>2.40E-03</td><td>3.38E-10</td><td>4.38E-05</td></tr><tr><td>rs850942</td><td>12:13044652-13084880</td><td>KIAA1467</td><td>1.01E-02</td><td>9.33E-20</td><td>9.90E-05</td></tr><tr><td>rs269782</td><td>5:139517088-139528554</td><td>CYSTM1</td><td>3.41E-03</td><td>2.11E-32</td><td>1.20E-04</td></tr><tr><td>rs12192704</td><td>6:30874410-30906415</td><td>DDR1</td><td>6.61E-04</td><td>1.85E-37</td><td>1.55E-04</td></tr><tr><td>rs4978941</td><td>9:112401576-112407320</td><td>SVEP1</td><td>3.13E-08</td><td>1.61E-14</td><td>3.30E-03</td></tr><tr><td>rs8041057</td><td>15:38118804-38146783</td><td>SRP14</td><td>9.96E-04</td><td>1.60E-28</td><td>4.44E-03</td></tr><tr><td>rs6744457</td><td>2:216119166-216286863</td><td>FN1</td><td>3.71E-06</td><td>1.85E-10</td><td>9.65E-03</td></tr><tr><td>rs10167593</td><td>2:71083169-71145381</td><td>NAGK</td><td>1.53E-03</td><td>5.13E-43</td><td>3.61E-02</td></tr><tr><td>rs3809472</td><td>15:43458687-43481812</td><td>SPATA5L1</td><td>1.24E-07</td><td>2.28E-55</td><td>4.00E-02</td></tr><tr><td>rs10982360</td><td>9:116468537-116473850</td><td>C9orf91</td><td>9.63E-04</td><td>3.25E-24</td><td>1.74E-01</td></tr><tr><td>rs1887784</td><td>9:116458405-116464475</td><td>C9orf91</td><td>6.83E-05</td><td>2.50E-22</td><td>2.66E-01</td></tr><tr><td>rs628383</td><td>3:152071522-152094779</td><td>CLRN1</td><td>4.74E-34</td><td>2.20E-14</td><td>2.67E-01</td></tr><tr><td>rs2798686</td><td>1:113355832-113417250</td><td>LRIG2</td><td>1.07E-02</td><td>4.45E-30</td><td>3.36E-01</td></tr><tr><td>rs12343516</td><td>9:122645199-122654702</td><td>PHF19</td><td>2.96E-06</td><td>2.11E-94</td><td>6.18E-01</td></tr><tr><td>rs8112960</td><td>19:21561436-21697351</td><td>ZNF100</td><td>8.66E-06</td><td>7.45E-09</td><td>8.47E-01</td></tr><tr><td>rs4916908</td><td>5:87600489-87768258</td><td>TMEM161B</td><td>6.53E-03</td><td>1.30E-28</td><td>8.50E-01</td></tr><tr><td colspan="6"><strong>Unknown mode of regulation</strong></td></tr><tr><td>rs1870832</td><td>6:89065388-89187471</td><td>CNR1</td><td>7.91E-12</td><td>1.68E-28</td><td>NA</td></tr><tr><td>rs17739675</td><td>2:64418705-64422285</td><td>LGALSL</td><td>5.01E-08</td><td>5.71E-17</td><td>NA</td></tr><tr><td>rs841603</td><td>12:27149351-27205593</td><td>C12orf71</td><td>1.36E-03</td><td>8.20E-28</td><td>NA</td></tr><tr><td>rs1868841</td><td>8:58568012-58666672</td><td>FAM110B</td><td>3.55E-05</td><td>3.06E-19</td><td>NA</td></tr><tr><td>rs1348478</td><td>5:119609042-119697096</td><td>PRR16</td><td>1.69E-02</td><td>9.64E-14</td><td>NA</td></tr></tbody></table>
<table><caption>**Table 1.** Cont.</caption><thead><tr><th>SNP<sup>1</sup></th><th>LncRNA region</th><th>Refseq gene</th><th>p-value for co-expression<sup>2</sup></th><th>p-value for ASE of lincRNA<sup>3</sup></th><th>p-value for ASE of Refseq gene<sup>4</sup></th></tr></thead><tbody><tr><td>rs7046236</td><td>9:70345772-70437996</td><td>TMEM252</td><td>4.83E-06</td><td>6.12E-24</td><td>NA</td></tr><tr><td>rs10772397</td><td>12:11256393-11295570</td><td>PRB3</td><td>1.04E-28</td><td>1.04E-51</td><td>NA</td></tr><tr><td>rs9423393</td><td>10:5266321-5295165</td><td>AKR1C4</td><td>2.19E-03</td><td>8.29E-41</td><td>NA</td></tr><tr><td>rs928736</td><td>21:33352005-33359159</td><td>OLIG1</td><td>1.31E-10</td><td>1.00E-27</td><td>NA</td></tr><tr><td>rs7132674</td><td>12:11443582-11530879</td><td>PRB2</td><td>3.28E-33</td><td>2.82E-11</td><td>NA</td></tr><tr><td>rs17739675</td><td>2:64475858-64534436</td><td>LGALSL</td><td>5.01E-08</td><td>1.14E-10</td><td>NA</td></tr><tr><td>rs10956365</td><td>8:128420701-128474058</td><td>POU5F1B</td><td>3.22E-11</td><td>1.23E-07</td><td>NA</td></tr></tbody></table>

SNPs having LD = 1.0 and two having LD> = 0.88. Table 2 summarizes the data for overlapping GWAS SNPs and cis-rSNP associated with ASE of lncRNA regions as well as the phenotypes to which these SNPs were associated in GWAS. Figure 3 exemplifies a lncRNA region with association between a SNP reported in GWAS and allele-specific expression. For additional examples see Figure S4. The association signals highlight 26 unique traits, out of which 9 are related to the immune system, which is not unexpected given that monocytes function in the immune system. Notably, three of the cis-rSNP that were significantly associated with ASE of a lncRNA are located in an intron of a protein coding gene, illustrating a direct or indirect regulatory function of an intronic SNP. One locus is reported by GWAS to be associated with a non-coding gene (genomic region: 2:19067188-19369296), which we confirm here by a cis-rSNP for the same lncRNA gene. Furthermore, when searching for overlaps between cis-rSNP associated with lncRNA and GWAS SNPs that are located in active enhancer regions in monocytes, we detected three GWAS SNPs that are located in a ChIP-seq peak for H3K4me1, and seven additional SNPs that are in the 2.5 kb enhancer regions flanking the histone mark. For H3K27ac no overlap between the cis-rSNP and the actual ChIP-seq peak was found, but five SNPs were located in the 2.5 kb flanking region. The number of overlapping GWAS SNPs – cis-rSNPs that map to an enhancer region are enriched at the same level as cis-rSNPs in general. For these GWAS SNPs the effect on the enhancer could be the functional mechanism.

For Refseq protein coding regions we found 458 loci where a reported GWAS SNP overlaps with a cis-rSNP or a proxy cis-rSNP in high LD (31 SNPs with LD 1.0 and 29 with LD>0.81 (Table S2). As for the lncRNAs a high fraction (n = 112) of the traits are related to immune diseases.

## Discussion

The combination of a large number of samples from primary monocytes and the sensitive genotyping method for detecting cis-regulatory SNPs based on allele-specific gene expression [4] renders our study well powered for detecting cis-regulatory SNPs that affect expression of long non-coding RNAs, which are expressed at a lower level than protein-coding genes. In the ASE approach transcript are detected using a genome-wide panel of SNPs markers, and hence the use of predefined hybridization probes for each transcript is circumvented. This unbiased

detection of expressed transcripts is an advantage of ASE analysis over expression microarrays. This advantage is shared by ASE and “next” generation transcriptome sequencing (RNA-seq). An advantage of RNA-seq compared to ASE analysis using SNP genotyping is that alternatively spliced transcripts and strand specific gene expression can be detected, provided that the sequencing coverage is sufficient. However, ASE analysis using RNA-seq suffers from the drawback of being highly affected by the reference bias [31,32], although methods that can correct for most of the bias have been developed recently [33,34]. Moreover, lowly expressed genes obtain very low coverage in RNA-seq.

In this study we focused our ASE analysis of monocytes on 5346 lncRNAs that were annotated by the ENCODE project as intergenic without overlap with protein coding genes. Thus we ensured that our analysis targeted truly lncRNA genes and was not confounded by overlapping protein-coding genes. Using stringent criteria for calling ASE, we identified 258 lncRNAs that were regulated by at least one cis-rSNP. We observed correlated expression of 20% of the lncRNAs with their neighbouring protein-coding genes, and co-regulation of a lncRNA and a neighbouring protein coding gene by the same cis-rSNP for approximately half of these lncRNA genes. This finding is consistent with a recent study of peripheral blood cells, where 75% of 112 lncRNAs that were mapped using traditional eQTL analysis were found to be regulated independently of nearby protein-coding genes [19]. Because the ASE analysis used in our study only detects cis-regulation, ASE analysis allows the dissection of cis-effects by lncRNAs that influence the transcription of nearby genes. If the expression level of a lncRNA and a nearby protein coding gene are correlated, and they both show significant ASE association with the same cis-rSNP, it is likely that the cis-rSNP directly regulates the protein coding gene. If the expression levels between a protein-coding and lncRNA gene are correlated, but there is no shared association with a cis-SNP, the lncRNA may be a regulatory factor that affects the expression of the adjacent protein coding gene. In our ASE data set we found that these two putative mechanisms occur at similar frequencies. However, using the ASE approach, trans-acting regulatory mechanisms for co-expressed genes would remain undetected. A recent study used traditional eQTL analysis of monocytes found several trans-regulated modules of co-expressed protein coding genes. However, lncRNA expression was not addressed in this study [35].

<footnote n="1">¹All 52 significantly co-expressed lncRNA and protein coding genes with a lncRNA associated cis-rSNP are listed.</footnote>

<footnote n="2">²Multiple testing correction using FDR of 10%,</footnote>

<footnote n="3">³The p-value cut-off for significant ASE is $10^{-6}$,</footnote>

<footnote n="4">⁴ASE p-value for a Refseq gene is shown for all regions that are co-expressed with a lncRNA and have an overlapping ASE analysis window, NA otherwise.
doi:10.1371/journal.pone.0102612.t001</footnote>
<table><caption>**Table 2.** SNPs associated with allele-specific expression of lncRNA windows with published trait- or disease-associations from genome-wide association studies.</caption><thead><tr><th rowspan="2">Genomic region of lncRNA</th><th colspan="4">ASE</th><th colspan="2">GWAS</th></tr><tr><th>Distance to GWAS SNP (bp)<sup>1</sup></th><th>In or close to histone mark</th><th>ASE p-value</th><th>Slope<sup>2</sup></th><th>SNP<sup>3</sup></th><th>Trait<sup>4</sup></th></tr></thead><tbody><tr><td>1: 22224268-22229983</td><td>16228</td><td>H3K4m1<sup>5</sup></td><td>2.41E-7</td><td>0.13</td><td>rs2501276</td><td>Immune response to smallpox vaccine (IL-6)</td></tr><tr><td>1: 22224268-22229983</td><td>128475</td><td></td><td>3.03E-7</td><td>0.07</td><td>rs16826658</td><td>Endometriosis</td></tr><tr><td>1: 41480834-41523120</td><td>177376</td><td></td><td>2.66E-12</td><td>0.09</td><td>rs6686842<sup>6</sup></td><td>Height</td></tr><tr><td>1: 220153613-220218488</td><td>13083</td><td>H3K4m1</td><td>8.46E-17</td><td>0.22</td><td>rs6687758</td><td>Colorectal cancer</td></tr><tr><td>2: 19067188-19369296</td><td>68959</td><td></td><td>1.36E-7</td><td>0.06</td><td>rs1876040</td><td>Cognitive test performance</td></tr><tr><td>2: 101944950-101970061</td><td>59999</td><td></td><td>1.32E-7</td><td>0.08</td><td>rs2310173</td><td>Ankylosing spondylitis</td></tr><tr><td>2: 166359613-166410847</td><td>68123</td><td></td><td>5.60E-9</td><td>0.06</td><td>rs6710518</td><td>Bone mineral density</td></tr><tr><td>2: 166359613-166410847</td><td>50321</td><td></td><td>1.82E-9</td><td>0.06</td><td>rs1346004<sup>7</sup></td><td>Bone mineral density</td></tr><tr><td>3: 157947829-158017517</td><td>91268</td><td></td><td>3.72E-25</td><td>0.08</td><td>rs12638253</td><td>Multiple sclerosis (severity)</td></tr><tr><td>6: 28237538-28245351</td><td>184924</td><td>H3K4m1<sup>5</sup><br>H3K27ac<sup>5</sup></td><td>2.25E-8</td><td>0.06</td><td>rs6903823<sup>8</sup></td><td>Pulmonary function</td></tr><tr><td>6: 30050868-30054162</td><td>7639</td><td>H3K4m1</td><td>5.90E-11</td><td>0.08</td><td>rs3893464</td><td>Graves' disease</td></tr><tr><td>6: 29802357-29824805</td><td>257099</td><td>H3K27ac<sup>5</sup></td><td>6.18E-7</td><td>0.08</td><td>rs4313034<sup>9</sup></td><td>Graves' disease</td></tr><tr><td>6: 30874410-30906415</td><td>204180</td><td></td><td>4.05E-12</td><td>0.11</td><td>rs4248154</td><td>Graves' disease</td></tr><tr><td>6: 58380311-58395738</td><td>21176</td><td></td><td>3.60E-21</td><td>0.07</td><td>rs9500256<sup>10</sup></td><td>Eosinophilic esophagitis (pediatric)</td></tr><tr><td>7: 1166536-1171429</td><td>161007</td><td></td><td>3.40E-14</td><td>0.07</td><td>rs10256972</td><td>Longevity</td></tr><tr><td>7: 7261310-7283935</td><td>26354</td><td></td><td>4.86E-9</td><td>0.07</td><td>rs10259085</td><td>Multiple sclerosis (severity)</td></tr><tr><td>7: 7261310-7283935</td><td>0</td><td>H3K4m1<sup>5</sup></td><td>3.41E-7</td><td>0.06</td><td>rs1299548</td><td>Amount of visceral adipose tissue adjusted for body mass index (BMI)</td></tr><tr><td>8: 23138679-23144384</td><td>0</td><td>H3K4m1<sup>5</sup><br>H3K27ac<sup>5</sup></td><td>9.45E-89</td><td>0.27</td><td>rs13278062</td><td>Age-related macular degeneration</td></tr><tr><td>9: 122645199-122654702</td><td>25619</td><td>H3K4m1<br>H3K27ac<sup>5</sup></td><td>4.09E-16</td><td>0.13</td><td>rs1953126</td><td>Celiac disease and Rheumatoid arthritis</td></tr><tr><td>9: 122645199-122654702</td><td>38017</td><td></td><td>1.02E-16</td><td>0.13</td><td>rs881375</td><td>Rheumatoid arthritis</td></tr><tr><td>9: 122645199-122654702</td><td>75358</td><td>H3K4m1<sup>5</sup><br>H3K27ac<sup>5</sup></td><td>4.13E-17</td><td>0.13</td><td>rs3761847</td><td>Rheumatoid arthritis</td></tr><tr><td>9: 135120874-135140438</td><td>0</td><td></td><td>3.39E-50</td><td>0.43</td><td>rs687621</td><td>D-dimer levels</td></tr><tr><td>9: 135120874-135140438</td><td>0</td><td></td><td>2.69E-66</td><td>0.43</td><td>rs657152</td><td>Liver enzyme levels</td></tr><tr><td>9: 135120874-135140438</td><td>0</td><td></td><td>2.69E-66</td><td>0.43</td><td>rs643434</td><td>Inflammatory biomarkers</td></tr><tr><td>9: 135120874-135140438</td><td>0</td><td></td><td>5.66E-52</td><td>0.44</td><td>rs612169</td><td>Metabolic traits</td></tr><tr><td>9: 135120874-135140438</td><td>0</td><td></td><td>1.00E-50</td><td>0.44</td><td>rs505922</td><td>Protein quantitative trait loci</td></tr></tbody></table>
<table><caption>**Table 2.** Cont.</caption><thead><tr><th rowspan="2">Genomic region of lncRNA</th><th colspan="2">ASE</th><th colspan="4">GWAS</th></tr><tr><th>Distance to GWAS SNP (bp)<sup>1</sup></th><th>In or close to histone mark</th><th>ASE p-value</th><th>Slope<sup>2</sup></th><th>SNP<sup>3</sup></th><th>Trait<sup>4</sup></th></tr></thead><tbody><tr><td>9: 135120874-135140438</td><td>0</td><td></td><td>3.80E-22</td><td>0.43</td><td>rs507666</td><td>Soluble ICAM-1</td></tr><tr><td>9: 135120874-135140438</td><td>3551</td><td></td><td>8.83E-26</td><td>0.43</td><td>rs579459</td><td>Soluble E-selectin levels</td></tr><tr><td>9: 135120874-135140438</td><td>4250</td><td></td><td>7.93E-25</td><td>0.43</td><td>rs495828</td><td>Angiotensin-converting enzyme activity</td></tr><tr><td>15: 43458687-43481812</td><td>30170</td><td></td><td>1.24E-21</td><td>0.14</td><td>rs2453533<sup>11</sup></td><td>Chronic kidney disease</td></tr><tr><td>15: 82975686-82986699</td><td>245274</td><td></td><td>4.86E-13</td><td>0.05</td><td>rs3743162<sup>12</sup></td><td>Alzheimer's disease (age of onset)</td></tr><tr><td>17: 26060783-26121168</td><td>150673</td><td>H3K4m1</td><td>1.28E-8</td><td>0.06</td><td>rs3760318</td><td>Height</td></tr></tbody></table>

A drawback of our study is that this specific dataset was not replicated by an independent method. However, in a previous study we have shown that a similar method for genome-wide ASE-analysis using NS-12 BeadChips (Illumina) was highly reproducible between replicate samples, with correlation coefficients of 0.9969 for gDNA and 0.9956 for cDNA. ASE-levels determined using NS-12 BeadChips and those determined by quantitative Sanger sequencing were also strongly correlated (0.86) in nine

representative genes [3]. Furthermore, in the current study the consistency of the ASE-signal for individual SNPs across the lncRNA regions is high, with an average standard deviation of 0.05 over all regions and samples. The consistency of the ASE-signals can be observed in Figure S4, where almost all SNPs in the lncRNA regions show overexpression of the same allele.

In our study we found that 32 of the cis-rSNPs that regulate the expression of lncRNAs have been identified in GWAS as risk

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Illustration of a region with a SNP from genome wide association studies (GWAS) which is associated with ASE of lncRNAs.
The tracks are from top to bottom in each panel: Horizontal red bars represent lncRNA transcript windows (with genomic coordinates) used for determination of ASE levels; grey lines show p-values for the association of GWAS SNPs with ASE levels in the transcript window; a grey line overlayed with a red dotted line indicates that a cis-rSNP overlaps with the reported SNP in the GWAS catalog; red vertical lines are median ASE-levels for each SNP.
doi:10.1371/journal.pone.0102612.g003</figcaption></figure>

<footnote n="1">¹The distance to the lncRNA region is 0 if the SNP is located within the region and the smallest distance otherwise,</footnote>

<footnote n="2">²Slope is given in absolute numbers,</footnote>

<footnote n="3">³Listed are all cis-rSNPs that are also found in the GWAS catalog together with the associated lncRNA,</footnote>

<footnote n="4">⁴The trait is taken from the GWAS catalog,</footnote>

<footnote n="5">⁵Within 2.5 kb,</footnote>

<footnote n="6">⁶rs6663565 as proxy,</footnote>

<footnote n="7">⁷rs2303393 as proxy,</footnote>

<footnote n="8">⁸rs6922111 as proxy,</footnote>

<footnote n="9">⁹rs7739434 as proxy,</footnote>

<footnote n="10">¹⁰rs13214831 as proxy,</footnote>

<footnote n="11">¹¹rs1153862 as proxy,</footnote>

<footnote n="12">¹²rs12442557 as proxy.</footnote>

<footnote>doi:10.1371/journal.pone.0102612.t002</footnote>
variants for human diseases or traits. The diseases or phenotypic
traits associated with these cis-regulatory SNPs according to the
lncRNA ASE analysis, include diseases of the immune system;
multiple sclerosis, Graves' disease, eosinophilic esophagitis, celiac
disease, rheumatoid arthritis, immune response to small pox
vaccine, soluble E-selectin levels, soluble ICAM-1, and ankylosing
spondylitis, for which monocytes are a relevant cell type (Table 1).
Chronic kidney disease is also a relevant trait as it can be caused
both by the autoimmune disease type 1 diabetes and by
hypertension, which are diseases for which monocytes are a
relevant cell type [25]. By manually classifying the traits in the
GWAS catalog into two classes, *ie* traits related to the immune
system, and other traits, we tested for over-representation of cis-
rSNPs in the GWAS-catalog using Fisher's exact test. For traits
associated to the immune system, the p-value for overrepresenta-
tion is 0.00027. The SNP rs13278062 that confers risk of age-
related macular degeneration detected in our study of monocytes
using ASE analysis (Table 1) was also highlighted by traditional
eQTL analysis of peripheral blood samples [19], while none of the
other detected risk SNPs or traits overlapped between the two
studies.

We also found that the identified *cis*-rSNPs for lncRNA are enriched to active enhancer regions in monocytes, which suggests a mechanism for their *cis*-regulatory functions. The reported 32 risk SNPs from GWAS that are strongly associated with the expression of non-coding RNAs provides interesting leads for further characterization and functional clues into immune diseases. Some of the GWAS SNPs are located in enhancer regions, which could be the cause of the allele specific expression of the lncRNA. Thus our study suggests more complex functional mechanisms underlying findings from GWAS than regulatory variants or expression levels of nearby protein coding genes, and provides novel insights into the relationship between genetic variation and human diseases.

**Supporting Information**

**Figure S1 Schematic picture of the principles ASE analysis.** In ASE the allele-specific expression level is measured by the difference in fluorophore signal intensity between the two alleles in the same sample. The average ASE-level is calculated for all heterozygous SNPs in the region. This ASE-value is then used in an association test against the genotypes of *cis*-rSNPs (shown in yellow). Figure adapted from Almlöf *et al* [4].
(PDF)

**Figure S2 Signal intensity threshold.** The figure shows the distribution of signal intensities for the two alleles. The x-channel represents the A-allele and the y-channel represents the B-allele. In the xraw panel the fraction of SNPs with the BB genotype having intensity levels above the cutoff level of 1000 is very low and similarly for the yraw panel.
(PDF)

**Figure S3 Regression lines.** Regression lines for the 32 risk SNPs from GWAS that are significantly associated with a lncRNA region.
(PDF)

**Figure S4 lncRNA in the GWAS catalog.** Illustration of all lncRNA regions that have a significant association to a SNP that is also significantly associated in the GWAS catalog. The tracks are from top to bottom in each panel: Horizontal red bars represent transcript windows (with genomic coordinates) used for determination of ASE levels; grey lines show p-values for the association of GWAS SNPs with ASE levels in the transcript window; a grey line

overlaid with a red dotted line indicates that this is the SNP that
overlaps with the reported SNP in the GWAS catalog; red vertical
lines are median ASE-levels for each SNP; annotated transcripts
are shown in black below the tracks.
(PDF)

**Table S1 Significantly associated lncRNA regions.**
(DOCX)

**Table S2 SNPs associated with allele-specific expression of Refseq protein coding genes with published trait-or disease-associations from genome-wide association studies.**
(DOCX)

**Methods S1 Quadratic normalization.** Detailed explanation of the quadratic normalization performed on the intensity levels obtained from the genotyping.
(DOCX)

**Acknowledgments**

We thank all the volunteers from the Cambridge BioResource who donated blood for this study (http://www.cambridgebioresource.org.uk/). SNP-genotyping using the Human 1.2 M Duo custom BeadChips v1 was performed by the SNP&SEQ Technology Platform in Uppsala, Sweden (www.genotyping.se). We thank the UPPmax NExT-generation sequencing Cluster & Storage (SNIC-UPPNEX, Uppsala, Sweden, www.uppmax.uu.se) for data storage and CPU-time for analysis. **Members of the Cardiogenetics Consortium:** Lead author for the Cardiogenetics Consortium: Heribert Schunkert, Heribert.Schunkert@uk-sh.de.

Tony Attwood¹, Stephanie Belz², Peter Braund³, Jessy Brocheton⁴,
Jason Cooper⁵, Abi Crisp-Hihn¹, Patrick Diemert (formerly Linsel-
Nitschke)², Jeanette Erdmann², Nicola Foad¹, Tiphaine Godefroy⁴, Jay
Gracey³, Emma Gray⁶, Rhian Gwilliams⁶, Susanne Heimerl⁷, Christian
Hengstenberg⁷, Jennifer Jolley¹, Unni Krishnan³, Heather Lloyd-Jones¹,
Ulrika Liljedahl⁸, Ingrid Lugauer⁷, Seraya Maouche²,⁴, Jasbir S Moore³,
Gilles Montalescot⁴, David Muir¹, Elizabeth Murray¹, Chris P Nelson³,
Jessica Neudert⁹, David Niblett⁶, Karen O'Leary¹, Helen Pollard³, Carole
Proust⁴, Angela Rankin¹, Augusto Rendon¹⁰, Catherine M Rice⁶, Hendrik
Sager², Nilesh J Samani³,¹¹, Jennifer Sambrook¹, Heribert Schunkert²,
Gerd Schmitz¹², Michael Scholz⁹, Laura Schroeder², Jonathan Stephens¹,
Stefanie Tennstedt (formerly Gulde)², Chris Wallace⁵

¹Department of Haematology, University of Cambridge, Long Road, Cambridge, CB2 2PT, UK and National Health Service Blood and Transplant, Cambridge Centre, Long Road, Cambridge, CB2 2PT, UK;
²Medizinische Klinik II, Universität zu Lübeck Germany;
³Department of Cardiovascular Sciences, University of Leicester, Glenfield Hospital, Groby Road, Leicester, LE3 9QP, UK;
⁴INSERM UMRS 937, Pierre and Marie Curie University (UPMC, Paris 6) and Medical School, 91 Bd de l'Hôpital 75013, Paris, France;
⁵Juvenile Diabetes Research Foundation/Wellcome Trust Diabetes and Inflammation Laboratory, Department of Medical Genetics, Cambridge Institute for Medical Research, University of Cambridge, Wellcome Trust/MRC Building, Cambridge, CB2 0XY, UK;
⁶The Wellcome Trust Sanger Institute, Wellcome Trust Genome Campus, Hinxton, Cambridge CB10 1SA, UK;
⁷Klinik und Poliklinik für Innere Medizin II, Universität Regensburg, Germany;
⁸Molecular Medicine, Department of Medical Sciences, Uppsala University, Uppsala, Sweden;
⁹Trium, Analysis Online GmbH, Hohenlindenerstr. 1, 81677, München, Germany;
¹⁰European Bioinformatics Institute, Wellcome Trust Genome Campus, Hinxton, Cambridge, CB10 1SD, UK;
¹¹Leicester NIHR Biomedical Research Unit in Cardiovascular Disease, Glenfield Hospital, Leicester, LE3 9QP, UK;
¹²Institut für Klinische Chemie und Laboratoriums medizin, Universität, Regensburg, D-93053 Regensburg, Germany.

**Data availability.** Genotypes and raw intensity data can be accessed at the European Genome-phenome Archive (EGA, http://www.ebi.ac.uk/ega/) which is hosted at the EBI, under accession number EGAS00000000119 as part of a previous study by our group [4].
## Author Contributions

Conceived and designed the experiments: JCA PL TP PD AHG WHO FC ACS. Performed the experiments: AL. Analyzed the data: JCA PL AL BG

## References

1. Pastinen T, Hudson TJ (2004) Cis-acting regulatory variation in the human genome. Science 306: 647–650.

2. Ge B, Pokholok DK, Kwan T, Grundberg E, Morcos L, et al. (2009) Global patterns of cis variation in human cells revealed by high-density allelic expression analysis. Nat Genet 41: 1216–1222.

3. Milani L, Lundmark A, Nordlund J, Kiialainen A, Flaegstad T, et al. (2009) Allele-specific gene expression patterns in primary leukemic cells reveal regulation of gene expression by CpG site methylation. Genome Res 19: 1–11.

4. Almlof JC, Lundmark P, Lundmark A, Ge B, Maouche S, et al. (2012) Powerful identification of cis-regulatory SNPs in human primary monocytes using allele-specific gene expression. PLoS One 7: e52260.

5. Deloukas P, Kanoni S, Willenborg C, Farrall M, Assimes TL, et al. (2013) Large-scale association analysis identifies new risk loci for coronary artery disease. Nat Genet 45: 25–33.

6. Verlaan DJ, Ge B, Grundberg E, Hoberman R, Lam KC, et al. (2009) Targeted screening of cis-regulatory variation in human haplotypes. Genome Res 19: 118–127.

7. Lerner-Ellis JP, Anastasio N, Liu J, Coelho D, Suormala T, et al. (2009) Spectrum of mutations in MMACHC, allelic expression, and evidence for genotype-phenotype correlations. Hum Mutat 30: 1072–1081.

8. Murphy A, Chu JH, Xu M, Carey VJ, Lazarus R, et al. (2010) Mapping of numerous disease-associated expression polymorphisms in primary peripheral blood CD4+ lymphocytes. Hum Mol Genet 19: 4745–4757.

9. Brown CJ, Ballabio A, Rupert JL, Lafreniere RG, Grompe M, et al. (1991) A gene from the region of the human X inactivation centre is expressed exclusively from the inactive X chromosome. Nature 349: 38–44.

10. Lee JT, Davidow LS, Warshawsky D (1999) Tsix, a gene antisense to Xist at the X-inactivation centre. Nat Genet 21: 400–404.

11. Rinn JL, Kertesz M, Wang JK, Squazzo SL, Xu X, et al. (2007) Functional demarcation of active and silent chromatin domains in human HOX loci by noncoding RNAs. Cell 129: 1311–1323.

12. Wang KC, Chang HY (2011) Molecular mechanisms of long noncoding RNAs. Mol Cell 43: 904–914.

13. Guttman M, Amit I, Garber M, French C, Lin MF, et al. (2009) Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature 458: 223–227.

14. Khalil AM, Guttman M, Huarte M, Garber M, Raj A, et al. (2009) Many human large intergenic noncoding RNAs associate with chromatin-modifying complexes and affect gene expression. Proc Natl Acad Sci U S A 106: 11667–11672.

15. Kim TK, Hemberg M, Gray JM, Costa AM, Bear DM, et al. (2010) Widespread transcription at neuronal activity-regulated enhancers. Nature 465: 182–187.

16. Orom UA, Derrien T, Beringer M, Gumireddy K, Gardini A, et al. (2010) Long noncoding RNAs with enhancer-like function in human cells. Cell 143: 46–58.

17. Dunham I, Kundaje A, Aldred SF, Collins PJ, Davis CA, et al. (2012) An integrated encyclopedia of DNA elements in the human genome. Nature 489: 57–74.

FC. Contributed reagents/materials/analysis tools: PD. Wrote the paper: JCA PL FC ACS.

18. Derrien T, Johnson R, Bussotti G, Tanzer A, Djebali S, et al. (2012) The Gencode v7 catalog of human long noncoding RNAs: analysis of their gene structure, evolution, and expression. Genome Res 22: 1775–1789.

19. Kumar V, Westra HJ, Karjalainen J, Zhernakova DV, Esko T, et al. (2013) Human disease-associated genetic variation impacts large intergenic non-coding RNA expression. PLoS Genet 9: e1003201.

20. Garnier S, Truong V, Brocheton J, Zeller T, Rovital M, et al. (2013) Genome-wide haplotype analysis of cis expression quantitative trait loci in monocytes. PLoS Genet 9: e1003240.

21. Howie BN, Donnelly P, Marchini J (2009) A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet 5: e1000529.

22. Frazer KA, Ballinger DG, Cox DR, Hinds DA, Stuve LL, et al. (2007) A second generation human haplotype map of over 3.1 million SNPs. Nature 449: 851–861.

23. Abecasis GR, Altshuler D, Auton A, Brooks LD, Durbin RM, et al. (2011) A map of human genome variation from population-scale sequencing. Nature 467: 1061–1073.

24. Li X, Yu Z, Jin B (2002) [Treatment of the paranasal sinuses mycosis with endoscopic sinus surgery]. Lin Chuang Er Bi Yan Hou Ke Za Zhi 16: 412–413.

25. Abbott A (2011) Europe to map the human epigenome. Nature 477: 518.

26. Creyghton MP, Cheng AW, Welstead GG, Kooistra T, Carey BW, et al. (2010) Histone H3K27ac separates active from poised enhancers and predicts developmental state. Proc Natl Acad Sci U S A 107: 21931–21936.

27. Rada-Iglesias A, Bajpai R, Swigut T, Brugmann SA, Flynn RA, et al. (2011) A unique chromatin signature uncovers early developmental enhancers in humans. Nature 470: 279–283.

28. Heintzman ND, Ren B (2009) Finding distal regulatory elements in the human genome. Curr Opin Genet Dev 19: 541–549.

29. Heintzman ND, Stuart RK, Hon G, Fu Y, Ching CW, et al. (2007) Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. Nat Genet 39: 311–318.

30. Hindorff LA, Sethupathy P, Junkins HA, Ramos EM, Mehta JP, et al. (2009) Potential etiologic and functional implications of genome-wide association loci for human diseases and traits. Proc Natl Acad Sci U S A 106: 9362–9367.

31. Stevenson KR, Coolon JD, Wittkopp PJ (2013) Sources of bias in measures of allele-specific expression derived from RNA-seq data aligned to a single reference genome. BMC Genomics 14: 536.

32. Degner JF, Marioni JC, Pai AA, Pickrell JK, Nkadori E, et al. (2009) Effect of read-mapping biases on detecting allele-specific expression from RNA-sequencing data. Bioinformatics 25: 3207–3212.

33. Satya RV, Zavaljevski N, Reifman J (2012) A new strategy to reduce allelic bias in RNA-Seq readmapping. Nucleic Acids Res 40: e127.

34. Lappalainen T, Sammeth M, Friedlander MR, t Hoen PA, Monlong J, et al. (2013) Transcriptome and genome sequencing uncovers functional variation in humans. Nature 501: 506–511.

35. Rotival M, Zeller T, Wild PS, Maouche S, Szymczak S, et al. (2011) Integrating genome-wide genetic variations and monocyte expression data reveals trans-regulated gene modules in humans. PLoS Genet 7: e1002367.