Human umbilical vein endothelial cells (HUVECs) from 3 different donors were isolated from umbilical cords obtained from the maternity ward of the Kuopio University Hospital and used at passages 5–8. The 4th replicate represents a pool of HUVECs purchased from Gibco and used as passage 10 (potential of >16 population doublings guaranteed). The gene expression profiles of all donors correlated strongly. This work was carried out in accordance with the recommendations of the Research Ethics Committee of the Hospital District of Northern Savo, Kuopio, Finland. Informed written consent was received from all participants and the experiments were performed in accordance with the relevant guidelines and regulations.

HUVECs were maintained in endothelial cell growth medium (EGM; basal medium with SingleQuots supplements CC-4133; Lonza) on cell culture flasks coated with 10 g/ml fibronectin (Sigma, St Louis, MO, USA) and 0.05% gelatin and maintained at 37°C and 5% CO₂. Hypoxia was induced in Ruskinn Invivo2 400 hypoxia workstation (Baker Ruskinn) in the presence of 1% O₂ and 5% CO₂ for 8 h. The 8 h timepoint was chosen to provide representation of early and late responses (12). Moreover, short exposure to hypoxia tends to promote cell survival and growth, while prolonged exposure to hypoxia leads to cell death (13). For adenoviral overexpression, cells were treated with AdCMV (empty vector), AdHIF1α and AdHIF2α (constitutively active forms of HIFs) (14).

For GRO-Seq nuclei isolation, cells were treated with cycloheximide (0.1 mg/ml) for 10 min, PBS-washed and incubated in 10 ml of swelling buffer [10 mM Tris-HCl, 2 mM MgCl₂, 3 mM CaCl₂, and 2 U/ml SUPERase Inhibitor (Thermofisher, Waltham, MA, U.S.A.)] for 5 min on ice. Cells were scraped and pelleted for 10 min at 400 × g and resuspended in 500 μl of swelling buffer supplemented with 10% glycerol. Subsequently, 500 μl of swelling buffer supplemented with 10% glycerol and 1% Igepal was added drop by drop to the cells under gentle vortexing. Nuclei were washed twice with lysis buffer (10 ml of swelling buffer supplemented with 0.5% Igepal and 10% glycerol), and once with 1 ml of freezing buffer (50 mM Tris-HCl pH 8.3, 40% glycerol, 5 mM MgCl₂ and 0.1 mM EDTA). Nuclei were counted manually using a Bürker chamber after Trypan Blue staining, centrifuged at 900 × g for 6 min and suspended to a concentration of 5 million nuclei per 100 μl of freezing buffer, snap-frozen in liquid nitrogen and stored −80°C until run-on reactions. For the run-on reaction, the nuclear run-on reaction buffer [NRO-RB; 496 μM KCl, 16.5 μM Tris-HCl, 8.25 μM MgCl2 and 1.65was preheated to 30°C. Then each ml of the NRO-RB was supplemented with 1.5 mM DTT, 750 μM ATP, 750 μM GTP, 4.5 μM CTP, 750 μM Br-UTP (Santa Cruz Biotechnology, Inc., Dallas, Texas, U.S.A.) and 33 μl of SUPERase Inhibitor (Thermofisher Scientific, Waltham, MA, U.S.A). Fifty microliters of the supplemented NRO-RB was added to 100 μl of nuclei samples, thoroughly mixed and incubated for 5 min at 30°C. RNA was then harvested by phenol-chloroform extraction [TRIzol LS (Thermofisher Scientific, Waltham, MA, U.S.A)].

For total RNA isolation, cells were treated with cycloheximide (0.1 mg/ml) for 10 min, PBS-washed and scraped into lysis buffer [1x Mammalian Polysome Buffer (Epicenter, Madison, Wisconsin), 1% Triton X-100, 1 mM DTT, 250 U/ml SUPERase Inhibitor, 7.1 U/ml Turbo DNase (Thermo Fisher Scientific, Waltham, Massachusetts, U.S.A.) and 0.1 mg/ml Cycloheximide] on ice. To assure complete lysis, the lysates were drawn up and expelled 4 times through a sterile 22–25 gauge needle. The cleared lysate was then treated with 10% SDS, snap-frozen in liquid nitrogen and stored −80°C.

GRO-seq libraries were subsequently prepared as previously described (15). The run-on products were treated with DNAse I according to the manufacturer's instructions (TURBO DNA-free Kit, Thermofisher, Carlsbad, CA, U.S.A.), base hydrolysed (RNA fragmentation reagent, Thermofisher, Carlsbad, CA, U.S.A.), end-repaired, and immuno-purified using anti-Br-UTP agarose beads (Santa Cruz Biotechnology, Inc., Dallas, Texas, USA).

Total RNA with size >200 nt for sequencing libraries was purified using the Zymo RNA Clean and Conc kit (Zymo Research, Irvine, California, U.S.A.) and rRNAs were removed using the Ribo-Zero Gold rRNA Removal Kit (Illumina, San Diego, CA, U.S.A.). This was followed by fragmentation (RNA fragmentation reagent, Thermofisher) and dephosphorylation.

Subsequently, both RNA-seq and GRO-seq RNA fragments were poly-A tailed (PolyA polymerase, New England Biolabs, Ipswich, MA, USA) according to manufacturer's instructions, followed by circularization and re-linearization. The cDNA templates were PCR amplified (Illumina barcoding) for 11–16 cycles. The GRO-seq libraries were size selected to be 180–300 bp in length. The RNA-seq libraries were selected to be 190–350 bp in length. The final, strand-specific libraries were quantified (Qubit dsDNA HS Assay Kit on a Qubit fluorometer, Thermofisher, Carlsbad, CA, USA) and pooled for 50 bp single-end sequencing with Illumina Hi-Seq2000 (GeneCore, EMBL Heidelberg, Germany).

Sequencing results were trimmed to remove 3' A-stretches originating from the library preparation and poor quality reads were filtered out (minimum 97% of bp over quality cutoff 10). Reads were aligned to the hg19 genome using bowtie allowing up to two mismatches and reporting only one alignment for each read. RNA-Seq data was mapped using STAR v2.5.4b to GRCh37/hg19 reference genome (16) with ENCODE standard options for long RNA-seq pipeline (1). Each sequencing experiment was visualized using custom tracks for the UCSC Genome Browser. R version 3.4.0 was used to filter the data, make calculations and create plots.

De novo lncRNA detection from GRO-Seq was performed using using Homer V4.9 (17) “findPeaks.pl” - algorithm with “–groseq” option. To separate lncRNAs from overlapping Refseq genes, “mergePeaks”-command was used with “-strand” option. Refseq annotated ncRNAs were further separated from protein coding genes based on “NR_” accession prefix (and further exclusion of those with “protein coding” annotation). To further divide lncRNAs to known and novel lncRNAs, the coordinates were intersected with the non-codeV5 database (18). When the overlap with database was over 70%, the lncRNA was assigned the non-coding RNA ID. SnoRNAs, ScaRNAs, and mature miRNAs were removed from the data.

The differential gene expression analysis was performed for transcripts that were expressed (RPKM >0.5) in at least 3 samples among the 8 studied (n = 4 for normoxia and hypoxia) using EdgeR (19). Similarly, the correlation between the lncRNA expression and the closest coding gene was performed for the coding genes that were expressed (RPKM >0.5) in at least 3 samples among the 8 studied. This criteria has been defined in order to remove low counts in the libraries to improve the sensitivity and the precision of the differential genes expression (20). Moreover, this threshold was selected because, for GRO-Seq, reads are counted throughout the gene body, which represents more total reads per genes than RNA-seq (12, 15, 21, 22). Differentially expressed genes were defined as transcripts that exhibited over 2-fold change in expression compared to control and FDR < 5%. Ingenuity Pathway Analysis (IPA; QIAGEN) or DAVID 6.8 (23, 24) was used to analyse the pathways and gene ontologies enriched among the differentially regulated genes under hypoxia (25).

LncRNAs were divided to promoter- and enhancer associated transcripts by measuring the average signal of H3K4me3 and H3K4me1 calculated from the datasets GSE29611 and GSE39089 around 1 kb of transcriptional start site. A log₂-ratio of H3K4me3 to H3K4me1 was calculated and positive ratios were assigned as promoter-associated transcripts and negative ratio as enhancer-associated (enhancer RNA) transcripts. Normal enhancers and super-enhancers (SEs) were detected from the normoxia and hypoxia datasets of public H3K27ac ChiP-Seq data (26) using the homer algorithm “findPeaks.pl” with “-style histone” and “-style super” settings, respectively. The eRNA expression was quantified for the combination of normal enhancers and super-enhancers from GRO-Seq and differential expression was determined using edgeR. SEs exhibiting FDR < 0.05 were selected for further analysis (Table S1). The SEs were detected for hypoxia and normoxia and those not overlapping were defined as gained/lost SEs upon hypoxia. To identify the stable lncRNA, we correlated the GRO-Seq and RNA-Seq data from two matching HUVEC donors. Highly stable lncRNAs were defined as a transcript that exhibited a RPKM >2 in RNA-seq. This analysis was performed separately for the two conditions. The expression of stable lncRNAs was compared to previously published microarray data (27) for the lncRNA genes that were represented by the Affymetrix HGU133 Plus2 array.

Atherosclerotic segments of femoral arteries from 4 patients with primary atherosclerotic lesions (n = 4), 5 patients with restenotic lesions (n = 5) were compared with non-atherosclerotic arteries from 4 patients (n = 4). The samples were age matched and average age (+SD) of the patients for atherosclerotic plaques and normal mammary artery controls were 70.9(+7.3) and 67.2(+9.1) years, respectively (27). Atherosclerotic samples were collected at vascular atherectomy operations (28). The normal samples represent trimmed ends of mammary arteries isolated during cardiac bypass surgery. Total RNA was isolated, amplified, labeled and hybridized to Affymetrix HGU133 Plus2 microarrays comprising 54 675 probe sets) essentially as recommended by the manufacturer (Affymetrix, Santa Clara, CA, USA) exploiting 3′ IVT Express Kit (Affymetrix, Santa Clara, California, US) for probe preparation and GeneChip® Hybridization, Wash, and Stain Kit (Affymetrix, Santa Clara, California, US) for hybridization and detection. The analysis of gene expression followed standard procedures including Robust Multichip Average (RMA) data normalization. After RMA normalization (29) of the microarray raw data, a filtering step was applied to remove the weakest signals (intensities lower than 2 × above global background) and from initial 54,675 probe sets 17,589 were included in Significance Analysis of Microarrays (SAM); (30) to identify differentially expressed genes at False Discovery Rate <0.1. The studies were approved by Local Ethical Committee of Kuopio University Hospital under identification number 53/2011 and the subjects gave informed consent.

The GRO-Seq data is available under NCBI Gene Expression Omnibus under accession number GSE103945. The public ChIP-Seq datasets analyzed can be found under accession numbers GSE29611 (H3K4me1: GSM733690 and H3K4me3: GSM733673 and H3K27ac: GSM733696) and GSE39089 (H3K4me1: GSM955981 and H3K4me3: GSM955983). The RNA-Seq from AdHIF1α and AdHIF2α overexpression are available under accession GSE98060. The microarray data from human femoral atherosclerotic lesions collected during atherectomy can be found under accession GSE53274 and GSE7307. The RNA-Seq data of hypoxia treated HUVECs and GRO-Seq data of adenovirally overexpressed constitutively active HIF1α and HIF2α have been submitted to under accession number GSE118530.

To study the nascent transcriptomes of ncRNAs in primary human endothelial cells, we performed GRO-Seq on HUVECs subjected to 8 h of hypoxia (1% oxygen). Altogether, we identified 33,508 transcripts above the expression threshold of 0.5 RPKM in at least 3 samples. Among these, 468 protein coding transcripts were found upregulated and 57 downregulated at least 2-fold (with false discovery rate below 5%) in response to hypoxia (Table S2A). We further analyzed these genes for gene ontology and upstream regulators using Ingenuity Pathway Analysis. As expected, HIF1α signaling and glucose metabolism were among the top canonical pathways enriched (Table S3) and HIF1α was identified as the top upstream regulator of the differentially expressed genes (Figure 1A). This confirms that a significant fraction of hypoxia-regulated genes are regulated on a transcriptional level.

In addition to protein-coding genes, we identified 1,763 differentially regulated lncRNAs in our analysis (applying a false discovery rate threshold of 5%), with 544 upregulated and 350 downregulated more than 2-fold. (Table S2B, Figure 1B). Eight hundred and eight of these lncRNAs corresponded to known lncRNAs (RefSeq or NONCODE) (Figure 1C) and 955 were novel, previously uncharacterized lncRNAs. Altogether, this demonstrates that hypoxia regulates a markedly larger set of non-coding RNAs than protein coding transcripts.

To further classify the lncRNAs based on their epigenomic features, we divided the transcripts to promoter- and enhancer-associated lncRNAs according to higher enrichment for H3K4me3 and H3K4me1 histone marks, respectively (31). Altogether, we were able to define 707 promoter-associated lncRNAs (p-lncRNA) and 1056 enhancer RNAs (eRNAs) (Figure 2A). Surprisingly, majority (90%) of the differentially regulated p-lncRNAs were induced upon hypoxia stimulus, in contrast to eRNAs that showed an approximately equal distribution of expression change directions (Figure 2A). To identify possible mechanisms underlying the differences in expression profiles of different types of lncRNAs, we performed de novo motif analysis near their transcription start sites. Clear differences were seen among the top transcription factor motifs enriched, with JUN/AP-1, MAF and NKX3 transcription factor motifs being enriched at the start sites of eRNAs, whereas NRF1, E2F, and SP1 were enriched at the promoters of p-lncRNAs (Figure 2B).

The transcriptional response to hypoxia is known to be controlled by the two master regulators, HIF1α and HIF2α, which are able to collaborate with distinct transcription factors (14, 32). To see whether differences in the usage of different HIF α-subunits contributes to differential activation of lncRNA subtypes, we further assessed the regulation of p-lncRNAs and eRNAs by overexpressing the constitutively active HIF1α and HIF2α proteins in HUVECs for 48 h. Western blot analysis of the HUVEC lysates demonstrated high expression of HIF1α and HIF2α protein levels relative to the AdCMV-transduced control cells (Figure S1). The GRO-Seq analysis demonstrated that HIF2α overexpression led to more induced p-lncRNAs (54%−140/260) compared to HIF1α (40%−104/260) (Figure 2C). The higher induction of lncRNA expression by HIF2α was further confirmed using RNA-Seq (Figure 2D).

Recent studies have suggested that the expression of cell type-specific genes is controlled through clusters of enhancers called super-enhancers (33). There is also accumulating evidence that these enhancer-dense genomic regions play key roles in cellular response to stimulus (34, 35). To study how hypoxia affects SEs, we used H3K27ac profiles from HUVECs subjected to normoxia and hypoxia for 24 h. Altogether, 1058 and 799 SEs were identified under normoxia and hypoxia conditions, respectively, which were highly concordant with the number of published basal HUVEC SEs (33). As the ChIP-Seq data lacked sufficient replicates and GRO-Seq data can also serve as a reliable indicator of enhancer activity, we used our GRO-Seq data from 4 independent replicates to identify which of these SEs are hypoxia-regulated. The results show that hypoxia led to significant regulation of 573 SEs (Table S1, Figure 3A) by inducing the expression of 358 SEs and repressing the activity of 215. Interestingly, the induced SEs were found nearby differentially expressed coding genes related angiogenesis (e.g., PGF, MMP2, DLL4, EGFL7, and TGFB1), regulation of apoptosis (e.g., NOTCH1, SOX4, ANKRD1, and DUSP6), cell adhesion (e.g., CD34, CDH5, VWF, and CTNNB1), and transcriptional activation (e.g., MEF2A, SOX4, EGLN1, MAFK, ASH2L, and TCEA2), whereas the downregulated SEs were close to regulators of peptidyl-tyrosine dephosphorylation (e.g., DUSP6, PTP4A2, PTPN14, and PTPN1) and genes involved in signal transduction (e.g., IL15, IL15RA, and MAP2K6) and GTPase activity (e.g., ARHGAP7/18/24/29; Table S4, Figures 3B,C). Notably, for 20% of SEs the closest RefSeq gene was a lncRNA, as exemplified by MALAT1 (Figure 3D) and LUCAT1.

We next correlated the fold change of lncRNAs and their nearby coding genes in response to hypoxia. The results demonstrated a good correlation with the coding gene that was stronger for eRNAs than for p-lncRNAs (Figure 4A). This supports the concept that coding genes could be co-regulated through sharing of regulatory motifs or cis-regulation mediated by lncRNAs (36).

To infer functions of the identified lncRNAs, we applied the guilt-by-association principle: if lncRNA shows an expression profile that correlates with nearby protein-coding genes involved in a given function, the lncRNA may be involved in the same function. Gene ontology analysis indicates that genes expressed in close proximity to hypoxia-regulated lncRNAs are mostly involved in angiogenesis and cell migration, as exemplified by the induction of pathways such as ephrin signaling, HIF1α and the VEGF family ligand-receptor interactions pathway (Figure 4B). This analysis thus predicts that lncRNAs could be important in hypoxia-mediated changes in endothelial cell signaling and function.

LncRNAs could also mediate trans-effects, which likely require the lncRNA transcript to be abundant or long-lived (37–39). To identify this subgroup of stable lncRNAs, we compared the GRO-Seq expression levels to RNA-Seq data from matching HUVEC donors. Here, RNA-Seq is expected to reflect the stable lncRNA pool within cell, whereas GRO-Seq is able to display all nascent transcripts irrespective of their half-life. Altogether, 13% of the lncRNAs seen in GRO-Seq were detected using RNA-Seq, and majority of these were associated with promoter-signature (among the 236 detected lncRNAs, 180 were promoter-associated lncRNAs and 56 were enhancer-associated lncRNAs). This was in contrast to protein-coding genes, for which 84% of GRO-Seq-detected genes exhibited detectable expression in RNA-Seq as well (data not shown). The low fraction of lncRNAs detected using RNA-Seq highlights the advantage of GRO-Seq in capturing a large amount of unstable lncRNAs that would otherwise be missed. Interestingly, the hypoxia regulation of eRNAs (fold change in GRO-Seq) exhibited a higher correlation with RNA-Seq than that of p-lncRNAs, suggesting that larger fraction of stable p-lncRNAs could be further regulated at the post-transcriptional level (Figure 4C, Figure S2). Supporting this, similarly to mRNAs, the p-lncRNAs have been shown to be targeted by cytoplasmic miRNAs and even act as competing endogenous RNAs (40), whereas the eRNAs are less likely to leave the nucleus.

We subsequently investigated whether the expression of the stable lncRNAs is altered during atherosclerosis. To this end, the atherosclerotic segments of femoral arteries from 4 patients with primary atherosclerotic lesions and 5 patients with restenotic lesions were compared with non-atherosclerotic mammary arteries from 4 patients (27). Altogether 30 hypoxia-regulated lncRNAs from our HUVEC analysis were represented by the probe sets of the microarray and found to be differentially regulated in the atherosclerotic samples compared to the control tissue. This included for example the well-known lncRNA MALAT1 (Figure 5) and the lncRNA HYMAI, which has been studied in relation to diabetes (41, 42) but not previously reported in association with atherosclerosis or hypoxia. Several lncRNAs with the notable ability to distinguish between healthy and diseased vascular tissue, such as LOC730101, KIAA1656, and LOC339803, remain poorly studied. However, LOC730101 has been implicated in the activation of Wnt/β-catenin signaling (43), a pathway that is critical for endothelial cell proliferation and migration (44). Notably, restenosis samples following the intravascular surgery of the primary lesion formed a separate cluster from the primary lesions, demonstrating further altered lncRNA profiles during intimal hyperplasia and arterial remodeling in the treated arteries. Altogether, our results indicate significant regulation of lncRNA expression in human atherosclerotic lesions which could contribute to atherosclerotic lesion progression or serve as indicators of lesion hypoxia.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The reviewer UR and handling editor declared their shared affiliation at the time of the review.