1. Genes and Chromosomes
  2. Human Biology and Medicine
Download icon

Transcriptional networks specifying homeostatic and inflammatory programs of gene expression in human aortic endothelial cells

  1. Nicholas T Hogan
  2. Michael B Whalen
  3. Lindsey K Stolze
  4. Nizar K Hadeli
  5. Michael T Lam
  6. James R Springstead
  7. Christopher K Glass
  8. Casey E Romanoski Is a corresponding author
  1. University of California, San Diego, United States
  2. University of Arizona, United States
  3. University of Western Michigan, United States
Research Article
Cited
0
Views
979
Comments
0
Cite as: eLife 2017;6:e22536 doi: 10.7554/eLife.22536

Abstract

Endothelial cells (ECs) are critical determinants of vascular homeostasis and inflammation, but transcriptional mechanisms specifying their identities and functional states remain poorly understood. Here, we report a genome-wide assessment of regulatory landscapes of primary human aortic endothelial cells (HAECs) under basal and activated conditions, enabling inference of transcription factor networks that direct homeostatic and pro-inflammatory programs. We demonstrate that 43% of detected enhancers are EC-specific and contain SNPs associated to cardiovascular disease and hypertension. We provide evidence that AP1, ETS, and GATA transcription factors play key roles in HAEC transcription by co-binding enhancers associated with EC-specific genes. We further demonstrate that exposure of HAECs to oxidized phospholipids or pro-inflammatory cytokines results in signal-specific alterations in enhancer landscapes and associate with coordinated binding of CEBPD, IRF1, and NFκB. Collectively, these findings identify cis-regulatory elements and corresponding trans-acting factors that contribute to EC identity and their specific responses to pro-inflammatory stimuli.

https://doi.org/10.7554/eLife.22536.001

Introduction

Atherosclerosis is an inflammatory disease of large arteries mediated by the accumulation of plaque within the vessel wall. Through sequelae such as heart attack, stroke, and peripheral vascular disease, it is responsible for an immense burden of morbidity and mortality. The pathogenesis of atherosclerosis involves several cell types and environmental risk factors (Lusis, 2000; Glass and Witztum, 2001). One of the critical cell types is the arterial endothelial cell (EC). The onset of atherosclerosis involves the activation of ECs by pro-inflammatory micro-environmental exposures including hemodynamic turbulence, oxidized-specific epitopes, and inflammatory cytokines (Tabas et al., 2015). These inflammatory stimuli result in the expression of adhesion molecules on the luminal EC surface and rolling, attachment, and migration of leukocytes into the vessel wall. Sustained recruitment and accumulation of immune cells in the vessel wall leads to extracellular matrix remodeling, smooth muscle cell migration, and the development of necrotic debris. Acute plaque rupture can result in sudden vascular occlusion, leading to heart attack or stroke.

Genome-wide association studies have identified more than 50 loci that predispose humans to cardiovascular disease (CVD) (Nikpay et al., 2015), of which the major cause is atherosclerosis. The majority of CVD loci reside outside protein-coding regions of the genome, suggesting that the risk variants alter gene regulatory function (Hindorff et al., 2009; Manolio, 2010). Still, the target genes, pathways, and cell types of action are largely unknown due to challenges in linking regulatory variants to function. A major challenge is that mammalian genomes contain upwards of a million potential regulatory elements called enhancers, yet a given cell type only utilizes on the order of tens of thousands of active enhancers (ENCODE Project Consortium, 2012; Andersson et al., 2014). This makes it difficult to accurately predict the functional cell systems and units of regulation from sequence alone (Shlyueva et al., 2014).

An important insight into enhancer biology is the observation that unique combinations of a few transcription factors (TFs) together activate cell-type-specific enhancers. Enhancer priming by TFs is both collaborative, (such that one TF will not bind its DNA motif if the motif for a collaborating TF is mutated [Heinz et al., 2013]), and hierarchical (the majority of sites bound by newly abundant TFs occur at enhancers pre-bound by collaborating TFs [Heinz et al., 2015; Romanoski et al., 2015; Kaikkonen et al., 2013]). This model is perhaps best characterized in the hematopoietic system and with toll-like receptor 4 signaling (Heinz et al., 2013; Kaikkonen et al., 2013; Heinz et al., 2010). For example, myeloid-specific enhancer activation and cell differentiation requires the TF PU.1 in combination with C/EBPb, whereas B cells require PU.1 in combination with EBF and E2A (Heinz et al., 2010).

In the current study, we take a genome-wide approach using DNA variation, epigenetic, and transcriptomic data to identify the major TF families that coordinate human aortic endothelial cell (HAEC) gene expression in homeostasis and upon exposure to prototypic inflammatory stimuli characteristic of atherosclerosis. Using a combination of experimental and computational approaches, we find that members of the ETS and AP1 TF families bind EC enhancers and that removing ETS member ERG elicits an inflammatory profile. We demonstrate that many enhancers identified in ECs are cell type-specific and several enhancers overlap with SNPs that have been associated to coronary artery disease (CAD) and hypertension. In addition, we demonstrate that TFs NRF2, NFκB, CEBD, and IRF1 are signal-dependent TFs that mediate the EC response to inflammatory stimuli.

Results

Transcription factors in the AP1 and ETS families dominate the enhancer landscape in HAECs

A total of 16,929 high-confidence enhancer-like elements were mapped in HAECs (Figure 1a) using chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) to identify promoter-distal elements marked by significant levels of histone H3 di-methylation of lysine 4 (H3K4me2) and acetylation of lysine 27 (H3K27ac) that together mark active enhancers (Heintzman et al., 2007; Creyghton et al., 2010; Rada-Iglesias et al., 2011). Chromatin accessibility, measured by Assay for Transposase Accessible Chromatin with high-throughput sequencing (ATAC-seq [Buenrostro et al., 2013]), was used to center the enhancer-like regions. The position within the element with the maximum signal that reflects greatest accessibility was used for centering. Using public global nuclear run-on sequencing (GRO-seq) data in HAECs (Kaikkonen et al., 2014), we observed that our set of enhancer-like loci produced bi-directional nascent RNA transcripts, or enhancer RNAs (eRNAs), as evidenced by the red and blue strand-specific RNA signals in Figure 1a. The potential function of eRNAs is not understood; however, eRNA output is robustly correlated with enhancer activity (Kaikkonen et al., 2014; Lam et al., 2013; De Santa et al., 2010; Kim et al., 2010), further supporting our enhancer set as active EC enhancers (Figure 1a).

Figure 1 with 3 supplements see all
HAECs display a distinct repertoire of enhancers that nominate combinations of the AP1, ETS, SOX and GATA TF families as major orchestrators of HAEC gene expression.

(a) A heatmap of 16,929 enhancer-like regions were selected by: accessible chromatin (ATAC-seq), coincident with H3K4me2 and H3K27ac deposition (ChIP-seq) in gene-distal positions (≥3 kb from promoters). Rows are enhancer loci, repeated for each data type in columns. Bidirectional transcription of enhancer RNAs (eRNAs) is also evident (GRO-seq). (b) The top four enriched motifs that occur in enhancers from a are shown. The transcription factor (TF) family, de novo motif matrix, percentage motif occurrence at enhancer loci versus random loci, and enrichment -logPvalues are indicated. Enrichment was calculated from 200 bp sequence, centered on chromatin accessibility. (c) The positional enrichment of the enriched motifs are shown relative to the center of the enhancer-like elements from a, where 0 bp is the center of the accessibility signal. (d) Gene expression measured by RNA-seq and limited to TFs is ranked by FPKM values to nominate the most highly expressed TFs in the AP1, ETS, SOX, and GATA families. (e, f) RNA expression, histone modifications, and super enhancer (SE) definitions are shown for: (e) ERG and (f) JUN loci. SE regions are highlighted in yellow and were defined using H3K27ac. More related data in Figure 1—figure supplements 13.

https://doi.org/10.7554/eLife.22536.002

We hypothesized that the major TFs that select and maintain enhancers in HAECs would be evident via enrichment of binding motifs in enhancer DNA sequences. Thus, we performed de novo motif enrichment analysis and discovered that AP1, ETS, SOX and GATA motifs were significantly enriched (-logPvalues > 7.1e2) in HAEC enhancers compared to random GC-matched genomic background sequence (Figure 1b, comprehensive list in Figure 1—figure supplement 1). Based on previous evidence (Heinz et al., 2013), we expected for functional motifs to be enriched near the maximum signal for chromatin accessibility. Indeed, AP1 and ETS were most frequently observed near the signal maximum, whereas the relationships for SOX and GATA motifs were less pronounced (Figure 1c). These data nominate roles for AP1 and ETS TF members as important mediators of the HAEC enhancer landscape.

Tens of genes encode proteins of the AP1, ETS, SOX and GATA TF families. Within each family different members share nearly identical DNA-binding domains and thus bind the same motif. In addition, AP1 protein members bind AP1 motifs as homo- and hetero-dimers. These sources of redundancy make it challenging to identify the functional family member(s) without additional information. To narrow the search in HAECs, we hypothesized that the operational TFs would be highly expressed and that their genetic loci would contain super-enhancers (SEs), or unusually dense clusters of highly decorated enhancers (Hnisz et al., 2013). Characterization of enhancer marks across cell-types has found SEs to be frequently located at loci encoding lineage-defining TFs (Whyte et al., 2013; Dowen et al., 2014). We queried rank-ordered expression data for TF family members and found that multiple members of each group were highly expressed (Figure 1d). For example, SOX members SOX18, SOX4, and SOX17 were among the top 4% most expressed TFs in HAECs. AP1 family members JUND, JUN, and JUNB were also in the top 4%. RNA-seq from other HAEC donors and replicate samples confirmed these findings (Figure 1—figure supplement 1). Next, we defined SEs using H3K27ac ChIP-seq data and found that, among others, the genetic loci for ERG (an ETS member) as well as AP1 members JUN, JUND, and JUNB harbored SEs (Figure 1e,f, Figure 1—figure supplement 2). Taken together, these data suggest that while multiple TFs from each family probably bind HAEC enhancers, that JUN, JUNB, JUND, and ERG likely serve prominent roles.

Roughly half of HAEC enhancers are endothelial-specific

To investigate which enhancer-like elements discovered in HAECs were specific to ECs, we analyzed public H3K27ac ChIP-seq datasets from ENCODE (ENCODE Project Consortium, 2012) and Roadmap Epigenomics consortia (Kundaje et al., 2015). Considering ECs are present in nearly all tissues, we focused on data collected in single cell types with the exception of ‘aorta’, ‘right ventricle’, ‘left ventricle’, and ‘right atrium’ that were included to observe their relationship to aortic endothelium. A total of 61 datasets were analyzed (Supplementary file 1). Human umbilical vein ECs, or HUVECs, were the only other EC type in the analysis. H3K27ac ChIP-seq tags were counted and normalized in each experiment at the 16,929 HAEC-defined enhancer loci. Hierarchical clustering resolved three distinct clusters of enhancers: an endothelial-specific set (n = 7405), a set common across cell types (1575), and a mixed set where only some cell types exhibited H3K27ac modification (7949) (Figure 1—figure supplement 3). Motif analysis of these three sets revealed differential frequencies of AP1, ETS, SOX, and GATA motifs (Figure 1—figure supplement 3). AP1 and ETS motifs were least frequently observed in the common enhancer set, while the ETS, GATA, and SOX motifs were most frequently observed in the endothelial-specific enhancer set. These data are consistent with the model that different combinations of transcription factors maintain cell-specific gene expression programs.

Aortic endothelial enhancers overlap genome-wide association SNPs for CAD and hypertension

To investigate whether EC enhancers have utility to prioritize non-coding functional variants for the cardiovascular diseases CAD and hypertension, we overlapped physical coordinates of the 16,929 enhancers from Figure 1a with GWAS associated variants. SNPs meeting genome-wide significance for CAD or hypertension, which is a major risk factor for atherosclerosis and CAD, were downloaded from the NHGRI-EBI GWAS Catalog (Welter et al., 2014). To account for linkage disequilibrium (LD, the correlation of alleles) between closely spaced SNPs on the same chromosome, we used 1000 Genomes data (Auton et al., 2015) to retrieve SNPs in LD with the reported GWAS SNPs when r2 was greater than 0.8 based on European haplotype structure. We identified 16 SNPs that were within HAEC enhancers (Table 1) and represent 22 lead SNPs from GWAS studies. Fifty percent of overlapping SNPs were within EC-specific enhancers (as opposed to those common or mixed across cell types), whereas only 43% of enhancers in HAECs are EC-specific (Figure 1—figure supplement 3). These data provide a focused list of potential functional non-coding variants that affect predisposition to CAD and hypertension through EC gene regulation. Further studies will be required to establish the regulatory consequence and predisposing mechanisms of these variants. Nonetheless, our evidence that perturbed endothelial expression contributes to vascular disease underscores the importance of elucidating endothelial gene regulatory programs in homeostasis and inflammatory environments.

Table 1

Overlap of HAEC enhancers with GWAS loci reported for coronary artery disease (CAD) or hypertension (HT). Associated SNPs were downloaded from the NHGRI-EBI Catalog of published genome-wide association studies. SNPs in linkage disequilibrium (LD) to GWAS association traits were calculated when r2 >0.8 according to the European reference population of the 1000 Genomes Project. HAEC enhancers defined in Figure 1a were overlapped by physical position (hg19 genome build). The GWAS SNP, p-value, GWAS trait, gene reported, PMID, overlapping HAEC enhancer coordinates and enhancer type are shown.

https://doi.org/10.7554/eLife.22536.006
GWAS SNPHAEC enhancer
SNP in enhancerLD to lead SNP from studyp-Value of leadTraitReported gene of leadPubMed IDPosition
(chr, start bp, end bp)
Nearest geneType
 rs12091564Lead2.0E-07CADHFE2216261371, 145395579, 145395699LOC101928979common
 rs72701850LD, rs12091564, r2 = 0.953462.0E-07CADHFE2216261371, 145396840, 145397006LOC101928979common
 rs72701850LD, rs10218795, r2 = 0.953462.0E-07CADHFE2216261371, 145396840, 145397006LOC101928979common
 rs56348932LD, rs17114036, r2 = 0.9168234.0E-19CADPLPP321378990, 242623251, 56988477, 56988661PLPP3EC-specific
 rs56348932LD, rs9970807, r2 = 0.9428682.0E-09CADPLPP3263433871, 56988477, 56988661PLPP3EC-specific
 rs56348932LD, rs17114046, r2 = 0.9428683.0E-07CADPLPP321846871, 213789881, 56988477, 56988661PLPP3EC-specific
 rs10047079LD, rs2229238, r2 = 0.8668487.0E-07CADILR6223190201, 154468114, 154468189SHEEC-specific
 rs55916033LD, rs10496288, r2 = 12.0E-09HTintergenic216261372, 83278987, 83279062LOC1720EC-specific
 rs55916033LD, rs10496289, r2 = 12.0E-09HTintergenic216261372, 83278987, 83279062LOC1720EC-specific
 rs72836880LD, rs10496288, r2 = 12.0E-09HTintergenic216261372, 83308909, 83309314LOC1720EC-specific
 rs72836880LD, rs10496289, r2 = 12.0E-09HTintergenic216261372, 83308909, 83309314LOC1720EC-specific
 rs112798061LD, rs10496289, r2 = 12.0E-09HTintergenic216261372, 83308909, 83309314LOC1720EC-specific
 rs3748861LD, rs13420028, r2 = 0.9162661.0E-10HTGPR39216261372, 133196310, 133196505GPR39mix
 rs3748861LD, rs10188442, r2 = 0.9162661.0E-10HTGPR39216261372, 133196310, 133196505GPR39mix
 rs144505847LD, rs6725887, r2 = 11.0E-09CADWDR1221378990, 242623252, 203672243, 203672412ICA1LEC-specific
 rs144505847LD, rs7582720, r2 = 13.0E-08CADWDR12242623252, 203672243, 203672412ICA1LEC-specific
 rs56155140LD, rs17087335, r2 = 0.9791125.0E-08CADNOA1, REST263433874, 57824385, 57824541NOA1mix
 rs5869162LD, rs6452524, r2 = 0.9246982.0E-07HTXRCC4216261375, 82393827, 82393921XRCC4EC-specific
 rs5869162LD, rs6887846, r2 = 0.9246982.0E-07HTXRCC4216261375, 82393827, 82393921XRCC4EC-specific
 rs6475604LD, rs7865618, r2 = 0.9405972.0E-27CADMTAP216061359, 22052677, 22052823CDKN2BEC-specific
 rs17293632LD, rs72743461, r2 = 11.0E-07CADSMAD32634338715, 67442510, 67442670SMAD3common
 rs17293632LD, rs56062135, r2 = 0.9884895.0E-09CADSMAD32634338715, 67442510, 67442670SMAD3common
 rs17227883LD, rs17228212, r2 = 0.9814382.0E-07CADSMAD31763444915, 67442769, 67443128SMAD3common
 rs1563966LD, rs1231206, r2 = 0.8441519.0E-10CADintergenic2137899017, 2095878, 2096222LOC101927839mix
 rs1563966LD, rs216172, r2 = 0.9093151.0E-09CADSMG6, SRR21378990, 2634338717, 2095878, 2096222LOC101927839mix
 rs7408563LD, rs7246657, r2 = 0.9005127.0E-06CADZNF3832387019519, 37808501, 37809067HKR1common

TF expression dynamics across 97 HAEC donors nominates three major modules of TFs as coordinating gene expression

We next questioned how AP1, ETS, SOX, and GATA TFs were expressed in artery ECs across the human population. We postulated that the most prominent actors would be highly expressed with modest variation between people. By leveraging global transcript levels collected across 97 genetically distinct HAECs from healthy human donors (Romanoski et al., 2010), we found that JUN and JUND (AP1) and ERG (ETS) exhibited the greatest median expression values with relatively little variability across the EC donor population (Figure 2a).

Figure 2 with 2 supplements see all
Coordinate gene expression across 96 genetically distinct HAEC donors identifies three regulatory programs among ETS, AP1, SOX, and GATA family members.

(a) TF gene expression measured by AffyHU133A array is shown across a population of 97 unique HAEC donors. Array probe set IDs were manually confirmed to cover expressed transcript isoforms of the indicated TFs based on RNA-seq data. Boxplots midlines are medians, box edges are 1st and 3rd quartiles, and whiskers 95% confidence intervals. (b) A heatmap of clustered pairwise Spearman correlation coefficients across 97 HAEC donors and module designations (colored sidebar) is shown with modules and pairs of interest highlighted by colored outlines. (c–f) Pairwise Spearman correlations between indicated mRNAs from b where each dot is a genetically distinct HAEC donor. More related data in Figure 2—figure supplements 1 and 2.

https://doi.org/10.7554/eLife.22536.007

To gain insight into the behavior of the TF members with respect to each other, we measured co-variation in TF gene expression profiles across the human population. Co-variation, or co-expression, of TFs could result from one TF (in)directly regulating another, both (in)directly regulating each other, or from each being regulated by a common third mechanism. By clustering pair-wide correlation coefficients across all TFs of interest, we identified three main groups with similar co-expressed profiles: group 1 (in orange) with members FOSL1 and ETS1; group 2 (in green) with GATA2 and GATA6; and group 3 (in yellow) with the remaining factors (Figure 2b, detailed examples in Figure 2c-f). The degree of correlation between TFs is indicated by red intensity (anti-correlation with blue; no correlation with white). Notably, TF expression of groups 1 and 2 were mostly anti-correlated with group three members. A very similar grouping of these TFs was observed when their relationship to all expressed genes was used as the clustering parameter (Figure 2—figure supplement 1). The result that FOSL1 and ETS1 are anti-correlated in expression with the remaining family members, and to HAEC transcripts overall, suggests that they promote opposing gene expression profiles in HAECs.

Nominated factors, including ERG and JUN, bind HAEC enhancers at closely spaced motifs

To test whether the nominated TFs indeed bound HAEC enhancers, we performed the first chromatin immunoprecipitation sequencing (ChIP-seq) experiments for JUNB, JUN, and ERG in HAECs and analyzed GATA2(ENCODE Project Consortium, 2012) and ETS1 (Zhang et al., 2013) binding data from human umbilical vein endothelial cells (HUVECs). The JUND cistrome would also be informative in these studies; however, we proceeded with JUN and JUNB because the heterodimeric binding of AP1 factors makes it likely that JUN and JUNB profiles encompass a major portion of the overall AP1 landscape. JUN, JUNB, ETS1, and GATA2 were all confirmed to bind active HAEC enhancers in the open chromatin region (Figure 2—figure supplement 2). As determined by clustering of binding profiles, JUNB and JUN were similar and ERG and ETS1 were similar, supporting the role of canonical DNA motifs on factor recruitment. Next, we asked if there was enrichment of other motifs proximal to the bound motifs as was observed previously for TF pairs known to collaboratively activate cell-specific enhancers (Heinz et al., 2010; Gosselin et al., 2014). For this analysis, loci for the bound factor (e.g. ERG) were centered on the respective motif (e.g. ETS motif) and sites lacking the motif were omitted. Then, frequencies of other motifs were calculated as a function of distance to the reference motif. We found that AP1 motifs were most frequently oriented within 50 base pairs to ERG-bound ETS motifs and that AP1 motif presence decayed with distance from the ETS motif (Figure 3a). Reciprocally, ETS motifs were most frequently observed proximal (within 50 base pairs) to JUN/JUNB co-bound AP1 motifs (Figure 3b). Both AP1 and ETS motifs were frequently observed near GATA2-bound GATA motifs; however, neither GATA nor SOX motifs were prominent in the vicinity of ETS or AP1-bound motifs (Figure 3c). These data support that AP1 and ETS factors collaborate to determine the active chromatin landscape in HAECs with GATA and SOX serving less active roles genome-wide. This observation is consistent with GATA and SOX motifs only having enrichment in EC-specific enhancers (Figure 1—figure supplement 3).

Figure 3 with 2 supplements see all
ERG and JUN co-bind EC enhancers and are enriched at EC-specific genes.

(ac) Promoter-distal regions bound by ERG, JUN, and JUNB, or GATA2 are shown in a, b, and c respectively in a one kilobase window. Each set was centered on the corresponding binding motif, and the frequency of other enriched motifs are shown on the y-axis. GATA2 binding was measured in HUVECs. (d) Allele-specific JUN binding (y-axis) as a function of allele-specific ETS motif mutations (colored lines). Each vertical line represents a single JUN peak identified via ChIP-seq. For a more complete explaination see the Figure 3—figure supplement 1 legend and methods section. (e) TF binding for JUN and ERG, the histone modifications H3K4me2 and H3K27ac and RNA abundance are shown at the genetic loci for VWF and NOS3. Promoters are highlighted in pink and JUN/ERG co-bound enhancers are highlighted in yellow. (f) The HAEC enhancer set from Figure 1a is annotated for JUN and/or ERG binding. Ingenuity Pathway Analysis results are shown for the genes nearest the 4248 JUN/ERG co-bound enhancers in the right panel. More related data in Figure 3—figure supplements 1 and 2.

https://doi.org/10.7554/eLife.22536.010

Allele-specific binding to chromosomes lacking motif mutations supports collaborative binding between AP1 and ETS factors

One approach to study collaborative binding between TFs is to knock-down/out a TF of interest and observe a shift in binding or activity at the regulatory element. To avoid complications in interpretability caused by potential redundancy of TF members, we took an alternative approach. As applied in inbred mouse strains previously (Heinz et al., 2013; Gosselin et al., 2014), we utilized naturally occurring genetic variation as a genome-wide source of motif mutations. The hypothesis is if the motif for a collaborative transcription factor is mutated then it should affect binding of the collaborating transcription factor whose motif remains in tact. To test this, whole-genome sequencing (WGS) of one HAEC donor was performed at an average of 40X coverage and the identified SNPs were phased with the appropriate 1000 Genomes (Auton et al., 2015) reference population (see Materials and methods). To quantify JUN binding to distinct homologous chromosomes at heterozygous loci, JUN ChIP-seq reads were iteratively mapped to human genome builds containing the appropriate allele. Sequence tags with discrepant mappings were omitted to avoid bias. For all loci with a JUN peak containing at least one heterozygous SNP, ChIP-seq tags were counted that could be uniquely assigned to one homologous chromosome. These data were then analyzed with respect to loci where only one SNP allele mutated either the AP1 or ETS motif.

Results showed that JUN binding was significantly affected by mutations in the AP1 motif (p = 5.0e-10) such that binding was predominant on the chromosome lacking AP1 motif mutations and diminished on chromosomes containing the mutation (Figure 3—figure supplement 1). Interestingly, JUN binding was also significantly affected by mutations in the ETS motif that occurred within 100 base pairs of the JUN peak center (Figure 3d p-value = 2.6e-3). We would expect to observe the reciprocal relationship, in which AP1 motif mutations alter ERG binding, but the ERG ChIP-seq experiment in the sequenced HAEC donor yielded less than ten thousand peaks and more information is necessary for this analysis. Taken together, these data support a collaborative relationship between AP1 and ETS factors at endothelial enhancers.

JUN and ERG co-occupy multiple elements near endothelial-specific genes

To interrogate the gene targets of JUN and ERG, we began with loci for genes expressed specifically or predominantly in ECs. All the genes queried, including vascular endothelial cadherin, (CDH5 or VE-cadherin), epidermal growth factor-like protein 7 (EGFL7), von Willebrand Factor (VWF), endothelial nitric oxide synthase (NOS3), and TEK receptor tyrosine kinase (TEK, or TIE2), exhibited between three and seven ERG/JUN co-bound enhancers across their genetic loci (Figure 3e, Figure 3—figure supplement 2). By cross-referencing cistromes of ERG and JUN that individually bound 35,559 and 63,312 genomic loci respectively, we found that 10,919 of the 16,292 (65%) high confidence enhancers from Figure 1a. were bound by one or both of ERG and JUN (Figure 3f). Each enhancer was assigned a target gene(s) based on the following criteria. Since nearest genes are not necessarily the target of enhancer activity, we incorporated expression quantitative trait loci, or eQTL, that were identified in HAECs (Romanoski et al., 2010). eQTL are SNP-gene pairs that describe a genetic locus whose alleles are associated with quantitative levels of gene expression values of the target gene, and thus provide a functional link between DNA sequence and gene regulation. Only 7% of enhancers harbored an eQTL SNP, in which case the associated gene was considered the target. In the remaining cases, the nearest gene was used. Pathway analysis for the resulting 4396 target genes for the 4248 ERG and JUN co-bound loci revealed significant enrichment in ‘Cardiovascular system development and function’ (p = 6.1e-37, Figure 3f). Together, these data support that ERG and JUN are major TFs at EC enhancers, and that their collaborative binding regulates expression of endothelial-specific genes important in vascular development and function.

ERG knockdown elicits a pro-inflammatory gene expression profile in HAECs

To test the functional importance of ERG on target gene expression, we knocked-down ERG using siRNA in HAECs and measured gene expression changes with RNA-seq and RT-qPCR (Figure 4, Figure 4—figure supplements 1 and 2). ERG RNA was reduced to less than 40% of normal levels in three independent experiments and resulted in differential expression of up to 1000 transcripts (>4-fold, FDR < 5%) by RNA-seq. Functional enrichment analysis demonstrated that ERG target genes are significantly annotated for ‘cell movement’, ‘breast or ovarian cancer’, ‘angiogenesis’, ‘development of vasculature’, ‘leukocyte migration’, and other pro-inflammatory functions (p-values from 1e-4 to 1e-17, Figure 4a,b).

Figure 4 with 2 supplements see all
ERG knockdown elicits a pro-inflammatory gene profile in HAECs.

(a) Genes up-regulated by ERG knockdown (≥4 fold, 5% FDR) were enriched in the indicated functional pathways using Ingenuity Pathway Analysis (IPA). Results are from three experiments using cells from different HAEC donors. ERG levels are shown in c. (b) A network of molecules regulated by ERG knockdown (red = up-regulated; blue = down-regulated) with membership (lines) to five functional processes. The predicted state of the functional processes is shown by color (activation in red; and inhibited in blue). Network analysis was performed using IPA and p-values correspond to pathway enrichments. (c) The log2 fold change caused by siERG knock-down is shown for endothelial-specific genes (blue) and pro-inflammatory genes (pink). Plotted is mean ± SD of fold changes measured in three HAEC donors. Per donor/gene effects are included in Figure 5—figure supplement 1. Asterisk indicates p-value<0.05 from a paired t-test of siERG to siSCR values across three donors. More related data in Figure 4—figure supplements 1 and 2.

https://doi.org/10.7554/eLife.22536.013

Among the most up-regulated genes caused by ERG knockdown were cytokines interleukin one alpha (IL1α; 16-fold), interleukin one beta (IL1β; 68-fold), leukemia inhibitory factor (LIF; 13-fold), interleukin 6 (IL6; fourfold), granulocyte colony stimulating factor (CSF3; 41-fold), transforming growth factor beta 2 (TGFβ2; fivefold) and other pro-inflammatory molecules including tissue factor (F3, 8-fold) (Figure 4b,c, Figure 4—figure supplements 1 and 2). In addition, EC-enriched genes that had multiple elements bound by ERG, such as CDH5, VWF, PECAM1, EGFL7, NOS3, and TEK were down-regulated upon ERG knockdown. To ensure that the inflammatory gene profile elicited by ERG knock-down was not a consequence of transfection itself or off-target effects, the profile resulting from six individual siERG oligos was measured along with two non-targeting scrambled siRNA controls and non-transfection controls (Figure 4—figure supplement 2). These data were reproducible and consistent with ablated ERG expression as the cause of pro-inflammatory expression profiles. Together, these data suggest that ERG normally functions to maintain EC-specific gene functions such as development and proliferation while at the same time suppressing inflammatory gene expression.

Oxidized phospholipids and inflammatory cytokines alter HAEC gene expression through signal-dependent changes to HAEC enhancer landscapes

To identify regulatory factors that instruct gene expression responses to inflammatory environments of atherosclerosis, we exposed HAECs to three pro-inflammatory treatments: (i) the oxidized products of 1-palmitoyl-2-arachidonoyl-sn-glycero-3-phosphocholine (oxPAPC) that are components of oxidized low-density lipoproteins (oxLDL) (Lee et al., 2012; Romanoski et al., 2011), (ii) tumor necrosis factor alpha (TNFα) that is a cytokine secreted largely by macrophages, and (iii) interleukin one beta (IL1β) that is released by many cell types including macrophages. Between 322 and 1174 genes were regulated by these exposures (>2-fold; 5% false discovery rate) (Table 2), and these genes were enriched in known response pathways (Figure 5—figure supplement 1). For example, 4 hr exposure to 40 µg/ml oxPAPC treatment resulted in up-regulation of genes belonging to the ‘NRF2-mediated oxidative stress response’ (p-value=2.e-13), and the ‘unfolded protein response’ (p-value=7.0e-7). Genes regulated by 4 hr exposures to TNFα and IL1β were highly enriched in the same pathways including, ‘TNF receptor signaling’ (p-values from 1.7e-14 to 1.8e-12), ‘granulocyte adhesion and diapedesis’ (p-values from 2.2e-9 to 7.1e-13), and ‘macrophage, fibroblast and EC roles in rheumatoid arthritis’ (p-values from 2.1e-12 to 1.7e-13).

Table 2

Molecular trait changes observed upon HAEC exposure to pro-inflammatory stimuli. Differential expression was determined in DESeq with duplicate RNA-seq experiments.

https://doi.org/10.7554/eLife.22536.016
StimulusRegulated genes
(>2 fold, 5%FDR)
de novo enhancers increased H3K27ac and accessibility upon stimulation
 oxPAPC 4 hr (40 µg/ml) versus control322 (242 up, 80 down)839
 TNFα 4 hr (2 ng/ml) versus control840 (611 up, 229 down)266
 IL1β 4 hr (10 ng/ml) versus control1174 (807 up, 367 down)3199

To better understand the program that coordinates HAEC response to oxPAPC, TNFα, and IL1β, we measured enhancer elements genome-wide. We next defined de novo or latent, enhancer-like elements that are genomic regions that became accessible and gained H3K27ac modification upon treatment. De novo enhancers result when signal-dependent transcription factors (SDTFs) play a critical role in the enhancer activation process and, in this study, were used in this study to identify SDTFs. We identified between 266 and 3199 de novo enhancers across treatments (Figure 5a, Table 2). To identity the SDTFs, we performed motif enrichment that revealed differential enrichment of TF motifs across treatments (Figure 5b; comprehensive list in Figure 5—source data 1). For example, the C/EBP, NFκB, and IRF motifs were preferentially enriched in TNFα and IL1β enhancers, whereas the anti-oxidant response element, or ARE, was enriched in the oxPAPC enhancer set. In all sets, we found AP1 and ETS motifs were highly enriched (p<1e-6), consistent with the model that the predominant AP1 and ETS endothelial factors collaborate with newly activated SDTFs to activate responsive enhancers and direct dynamic gene expression (Figure 5a,b). Enrichment of κB motifs at TNFα and IL1β enhancers is consistent with previous work demonstrating that NFκB is a master transcription factor of inflammatory gene programs in s as well as other cell types (Brown et al., 2014). Likewise, oxPAPC-induced enrichment of the ARE motif, to which the TF nuclear factor, erythroid 2-like 2, or NFE2L2/NRF2 binds, is consistent with reports of single gene targets that the transcription factor NRF2 regulates the response to oxidative stress (Ma, 2013; Jyrkkänen et al., 2008).

Figure 5 with 4 supplements see all
Inflammatory signals activate enhancer-like elements with distinct motifs and suggest that CEBP and IRF members mediate responsiveness to TNFα and IL1β.

(a) A schematic for de novo enhancer activation by combinations of signal-dependent and collaborative TFs. (b) Enrichments of the NRF2, NFκB, IRF, C/EBP, AP1, and ETS motifs (y-axis) are shown for de novo enhancer sets activated after 4 hr by either 40 µg/ml oxPAPC (ox), 10 ng/ml TNFα, or 2 ng/ml IL1β (x-axis). Blue bars show the percent of de novo enhancers containing the motif and red bars indicate the percentage of GC-matched random genome sequence containing the motif. Motifs were considered within 200 basepair windows of enhancers. (c) ChIP-seq analysis of p65-binding sites shows regions of the genome that are co-bound by p65 and JUN, p65 and ERG, all three, or solely p65 post treatment by 10 ng/ml IL1β and 2 ng/ml TNFα treatment. (d) Gene expression, measured by RNA-seq is shown by heatmap for TFs of the C/EBP and IRF families. Expression values are shown from two HAEC donors and replicate samples. (e) RT-qPCR analysis shows CEBPB, CEBPD, and IRF1 expression from 0 to 24 hr post treatment by 10 ng/ml IL1β and 2 ng/ml TNFα treatment. Plotted are mean ± S.D., experiment performed with biological triplicates. * represents p<0.05 by t-test. More related data in Figure 5—figure supplements 14.

https://doi.org/10.7554/eLife.22536.017

To directly test the hypothesis that NRF2 and NFκB were in fact SDTFs responsible for inflammatory gene responses, we measured the NRF2 cistrome upon oxPAPC treatment and the NFκB cistrome (using ChIP-seq with κB-component p65) in TNFα and IL1β 4-hr-treated ECs. NRF2 and p65 binding were associated with increases in H3K27ac on adjacent nucleosomes consistent with a role for NRF2 and NFκB as activators of target gene expression (Figure 5—figure supplements 2 and 3). Genome-wide binding of NRF2 with oxPAPC treatment identified 6923 NRF2-bound peaks that overlapped with 21% of oxPAPC-elicited de novo enhancers. These included gene loci for known targets including heme oxygenase 1 (HMOX1), thioredoxin reductase 1 (TXNRD1), NAD(P)H quinone dehydrogenase 1 (NQO1), and glutamate-cysteine ligase modifier subunit (GCLM) (Figure 5—figure supplement 2). The top four enriched motifs in the NRF2 cistrome were the AP1 motif, an ARE, a nuclear receptor element (NRE) and the ETS motif. These data suggest that, whereas NRF2 is a significant component of the oxPAPC gene response, several additional TFs likely coordinate gene responses at the chromatin level.

For NFκB, we identified tens of thousands of bound loci after TNFα and IL1β treatment (75,937 and 51,040, respectively). Upon cytokine treatment, 52–63% of the elements that gained p65 were pre-bound in untreated HAECs by ERG or JUN (Figure 5c), consistent with the working model that enhancers are selected by lineage-restricted combinations of factors that direct signal-dependent transcription factor binding profiles (Romanoski et al., 2015). As for de novo enhancers, NFκB had a major presence and bound to 86% and 55% of those elicited by TNFα and IL1β treatments, respectively. Motif enrichment of NFκB cistromes in HAECs identified a prominent role of the AP1 motif as well as roles for ETS, IRF, and CEBP factors (Figure 5—figure supplement 3, discussed below).

In addition, we analyzed allele-specific NFκB binding as a function of motif mutations in the κB motif itself as well as mutations in the AP1 and ETS motif. Interestingly, we observed that NFκB binding was diminished at loci where κB, ETS, or AP1 motifs were mutated (Figure 5—figure supplement 4). These data are consistent with collaborative interactions between ETS, AP1, and NFκB in establishing inflammatory expression profiles in human ECs, and offer a mechanism whereby ECs may generate cell-specific transcriptional responses to environmental stimuli. Furthermore, these data demonstrate that allele-specific binding in heterogeneous human cells is a useful means to reveal collaborative interactions between transcription factors.

Non-random spacing of ETS and κb motifs at co-bound elements suggest interplay between ERG and NFκB in inflammation

The transcriptional signature resulting from ERG knockdown in aortic ECs (Figure 4) is evidence that ERG performs an anti-inflammatory role. This has been demonstrated previously, where ERG suppressed IL-8 and NFκB-mediated inflammation in HUVECs (Dryden et al., 2012; Sperone et al., 2011; Yuan et al., 2009). One model to explain the relationship between anti-inflammatory effects of ERG and the pro-inflammatory effects of NFκB at co-bound loci involves the observation that ERG levels are decreased upon inflammatory stimuli ([Yuan et al., 2009] Figure 6—figure supplement 1a). Depletion of ERG, simultaneous with increased NFκB concentrations, could result in a functional switch caused by stoichiometric competition between factors. At the pro-inflammatory target gene ICAM-1, Sperone et. al. demonstrated putative ETS-binding sites within the NFκB motif at the gene promoter(Sperone et al., 2011). We also observe co-occupancy of ERG and NFκB (p65) at the ICAM-1 promoter (Figure 6—figure supplement 1b).

To examine the genome-wide relationship between ERG and NFκB occupancy upon inflammatory signaling, we compared binding profiles. Twenty-nine percent of ERG binding sites in untreated ECs gained p65 binding upon 4-hr treatment with TNFα, and conversely, 25% of TNFα-elicited p65 binding occurred at loci pre-bound by ERG (Figure 6—figure supplement 1c). All co-bound loci were centered on the presumed ETS motif to which ERG binds and the distribution of NFκB motifs within 100 base pairs was calculated. This analysis revealed a distance relationship between ETS and NFκB motifs that is consistent with ERG and NFκB affecting each other’s binding and activity to promote pro-inflammatory gene expression (Figure 6—figure supplement 1d). Importantly, however, we do not observe reduction in ERG binding after TNFα treatment, as would be expected if the factors were competing for the same element. This is exemplified at the ICAM-1 promoter and in the variability of TNFα-induced changes to ERG and NFκB binding genome-wide (Figure 6—figure supplement 1e). Together, these data support coordinated regulation of inflammatory pathways by ERG and NFκB that likely involves multiple mechanisms.

A role for CEBPD and IRF1 in the HAEC response to inflammatory cytokines

The result that CEBP and IRF motifs were preferentially enriched in TNFα and IL1β-induced enhancers suggested that TFs in these families direct the EC response to cytokines (Figure 5b). To identify the likely members, we examined relative expression levels with and without cytokine exposure and found expression of CEBP and IRF factors to be relatively low in untreated HAECs. Upon TNFα and IL1β exposure, however, CEBPD and IRF1 transcription were highly induced (Figure 5d). In order for CEBPD and IRF1 to participate in enhancer activation, we reasoned they would need to be expressed prior to 4 hr when de novo motifs were measured. Indeed, a time-course experiment confirmed that both CEBPD and IRF1 RNAs were induced after cytokine treatment with a peak expression at 2 hr (Figure 5e).

Next, we measured binding of IRF1 and CEBPD genome-wide after cytokine treatment (Figure 6a). For IRF1, we observed binding to the minority (<10%) of de novo enhancers, and when we did observe binding it was most frequently co-bound with either NFκB or CEBPD. On the other hand, CEBPD bound to a significant proportion of IL1β and TNFα de novo enhancers (34% and 28%, respectively), mostly in conjunction with NFκB. Of all de novo enhancers, 24–26% were co-bound by NFκB and CEBPD within 200 basepairs of each other, strongly suggesting that these factors together regulate target gene expression.

Figure 6 with 1 supplement see all
CEBPD and IRF1 knockdown effect on gene expression response to TNFα.

(a) Binding of p65, CEBPD, and IRF1 was measured by ChIP-seq and is shown at de novo enhancers elicited by IL1β (e) and TNFα (f). Factor binding was measured by CHIP-seq upon 4-hr cytokine treatment, with binding of CEBPD measured after IL1β only. (b) Factor binding, H3K27ac, and mRNA expression is shown at the CEBPD locus in control-treated and TNFα-treated HAECs. CEBPD mRNA is shown in the bottom three tracks as a function of control knockdown (siSCR), CEBPD knockdown, and ERG knockdown. (c) A global view of mRNA responses to TNFα as a function of CEBPD knockdown (y-axis) compared to control (x-axis); mean ± S.D. (d) The response to TNFα is shown for molecules of interest in control (siSCR) and CEBPD knockdown HAECs. (e) The IFIT locus, highlighting elements co-bound by IRF1 and NFκB (p65). More related data in Figure 6—figure supplement 1.

https://doi.org/10.7554/eLife.22536.023

Notably, ERG and NFκB bind to elements at the CEBPD locus and ERG knockdown induces CEBPD expression (Figure 6b). These data are consistent with a model whereby cytokine treatment causes down-regulation of ERG as well as nuclear entry and binding of NFκB. These two events induce CEBPD expression and enable CEBPD and NFκB to co-bind enhancers of target genes.

Lastly, we tested whether reducing CEBPD and IRF1 levels with siRNA would dampen the endothelial response to cytokines. As hypothesized, CEBPD knockdown to less than 20% control RNA levels resulted in dampened fold-change ratios for genes up-regulated by TNFα (250 genes up-regulated in siCEBPD compared to 497 in scrambled control; Figure 6c, bottom triangle). Genes in this set contained inflammatory molecules including vascular cell adhesion molecule 1 (VCAM1), E-selectin (SELE), IL1A, IL1β, and F3 (Figure 6d). However, lesser induction was partly due to an increase in these molecules in the untreated state (Figure 6—source data 1) suggesting that CEBPD may play an anti-inflammatory role at baseline. In the IRF1 knockdown, the most prominent locus affected was a stretch of related antiviral proteins called ‘interferon-induced proteins with tetratricopeptide repeats’ (IFITs) on chromosome 10q23 (Figure 6e). We found that IRF1 bound six elements along this locus including at the promoters of IFIT2, IFIT3, IFIT1, and IFIT5. IRF1 knockdown reduced the basal levels of RNA from these genes and likewise prevented induction upon TNFα treatment, suggesting that IRF1 plays a critical role in their regulation. Together, these data provide evidence that CEBPD and IRF1 tune cytokine-induced regulatory function that is dominated by NFκB. However, further studies will be required to understand the precise role that IRF1 and CEBPD have in aortic EC gene regulation.

Discussion

In this study, we provide the first detailed characterization of human aortic endothelial enhancers at baseline and under inflammatory conditions using a combination of genome-wide approaches. We observe that ETS, AP1, GATA, and SOX motifs are enriched in active EC enhancers, and we provide evidence that ETS and AP1 factors bind the majority of active enhancers in aortic ECs. Allele-specific binding paired with allele-specific motif mutations provided further evidence of collaborative binding between AP1, ETS, and κB factors. Our work demonstrates that knockdown of the ETS factor ERG results in a pro-inflammatory expression profile and corresponding down-regulation of EC-enriched genes. Further, we identify several hundred de novo enhancers formed in response to pro-inflammatory molecules abundant in the atherosclerotic plaque: namely, oxidized phospholipids and the cytokines TNFα and IL1β. Motif enrichments paired with expression changes prioritized NFκB, CEBPD, and IRF1 as the responsible coordinators of cytokine response and NRF2 as a coordinator of response to oxidized phospholipids. Our work provides valuable context to the role of these transcription factors within the regulatory networks controlling the endothelium in physiologic and disease states.

Our approach was to learn the regulatory lexicon of ECs beginning with H3K4me2-positive, H3K27ac-positive, accessible DNA sequence. In doing so we identified a set of TF families, each of which possesses many highly expressed members and inter-correlated expression patterns across aortic ECs from 97 people (Figure 2). Of note, whereas SOX members (SOX4/17/18) were among the most abundant of all TFs in aortic ECs, the SOX motif was not as highly enriched at EC enhancers as AP1 and ETS motifs (Figures 1b, c and 3a). An explanation for this could be that SOX factors are critical for the selection of key endothelial enhancers in lineage development, and that only some of these enhancers remain open and available for other factor binding. This would explain the moderate enrichment of SOX motifs and is consistent with the role of SOX17 and SOX 18 in angiogenesis and vascular integrity (reviewed by [De Val and Black, 2009]). However, we acknowledge that this theory does not explain why ECs persistently express such high levels of SOX mRNAs. Another explanation, consistent with the selective enrichment of SOX and GATA motifs only in EC-specific enhancers (Figure 1—figure supplement 3) is that these factors play localized but important roles on target genes.

The finding that the ETS and AP1 factors ERG and JUN together bind the majority of EC enhancers, and co-bind key EC loci, supports a prominent and collaborative role for these factors in selecting endothelial specific enhancers. Our analysis of motif mutations further exemplifies a co-dependence of these factors in enhancer binding. While we focus on ERG and JUN in particular, other combinatorial interactions between other ETS and AP1 members almost certainly play prominent roles in endothelial regulation. Still, our findings are consistent with previous reports that ERG is among several ETS transcription factors critical for endothelial lineage development (De Val and Black, 2009; Nikolova-Krstevski et al., 2009; Shah et al., 2016) and regulation of candidate endothelial genes such as CDH5, VWF, and eNOS (Yuan et al., 2009; Birdsey et al., 2008; Laumonnier et al., 2000). Knockdown of ERG in aortic ECs revealed two notable effects, namely a reduction in EC-specific gene expression, and production of a pro-inflammatory transcriptional profile. The first of these effects is in keeping with our results demonstrating that ERG and AP1 co-bind promoters and enhancer elements at key endothelial gene loci. This supports our model whereby ERG binds collaboratively with AP1 factors to drive a basal lineage-defining transcriptional network in ECs.

Functional interactions between ETS and AP1 family members have been previously described (Bassuk and Leiden, 1995), and cooperative binding of ETS factors and FOS/JUN occurs at adjacent ETS and AP1 motifs with variable spatial orientation (Kim et al., 2006). In the current study, we also observe a spike in AP1 motifs near ETS motifs consistent with precise spatial orientation at some loci (Figure 3a,b), although this motif relationship does not describe the majority of genome-wide ERG and JUN co-binding we observe.

The transcriptional signature resulting from ERG knockdown in aortic EC also further supports the anti-inflammatory role for ERG, which has been demonstrated in HUVECs (Dryden et al., 2012; Sperone et al., 2011; Yuan et al., 2009). Our observation of regularly spaced ETS and κB motifs at co-bound loci suggests interplay between these two TF families. However, we do not observe coordinated change in genome-wide binding for ERG and NFκB upon TNFα treatment, indicating that these factors act in ways other than competitors for binding (Figure 6—figure supplement 1). In addition, we demonstrate ERG to regulate expression of numerous transcription factors (Supplementary file 2) that likely regulate secondary transcriptional targets. Recent work in murine EC has also demonstrated a role for ERG in promoting vascular integrity through promotion of canonical Wnt-signaling, via stabilization of β-catenin (Birdsey et al., 2015). This highlights the potential that ERG influences endothelial function through multiple mechanisms. Our present work underscores a role for ERG in promoting endothelial homeostasis and in regulating the inflammatory response. It also broadens our perspective of its regulatory function and cooperation with other factors genome-wide.

Our finding that the CEBP motif was enriched at TNFα and IL1β-induced enhancers and that CEBPD transcript levels were highly induced after treatment supported its role in binding inflammatory enhancers (Figure 5). CEBPD knockdown in aortic ECs in the absence of TNFα or IL1β caused modest up-regulation of pro-inflammatory molecules including ICAM1, SELE, IL1α, and IL1β (Figure 6, Figure 6—source data 1), suggesting that CEBPD maintains an anti-inflammatory expression profile in resting ECs. This is consistent with an anti-inflammatory role of CEBPD in pancreatic beta cells in the setting of cytokine stimulation (Moore et al., 2012). However, upon TNFα and IL1β stimulation in aortic ECs, CEBPD knockdown also dampened this inflammatory response as compared to untreated ECs suggesting that in fact CEBPD is required for full responsiveness to cytokines (Figure 6). This is consistent with previous reports showing CEBPD to be up-regulated in many inflammatory settings including atherosclerosis (Takata et al., 2002; Ko et al., 2015). Going forward, a key part of elucidating the role of SDTF such as CEBPD will be to quantify their loss-of-function effect on de novo enhancers that arise upon inflammatory stimuli. Given our hypothesis of hierarchical transcription factor binding, we would expect to observe loss of some of these enhancers following SDTF inhibition. Complicating this approach, however, is that inflammatory stimuli often induce expression of SDTFs. For example, NFκB binds the promoter of CEBPD (Figure 6b). Therefore, targeted mutagenesis strategies, such as the CRISPR-Cas9 system, will be required to fully eliminate SDTF function and interpret the data. Overall, these findings motivate further experiments to fully understand the regulatory function of CEBPD.

Cell type and context-specific enhancer mapping is a critical approach toward fully understanding regulatory disease mechanisms, as many disease loci reside in non-coding DNA sequence. Nonetheless, enhancers are a challenge to measure in pure human cell types relevant to many diseases. We provide a list of candidate SNPs with potential EC-specific non-coding regulatory function (Table 1). This is an important step toward fully understanding the etiology and molecular pathogenesis of CAD in humans. In conclusion, our study of dynamic endothelial enhancer elements is an advancement toward a systems-levels understanding of vascular inflammatory diseases.

Materials and methods

Cell culture

HAEC were isolated as described (Navab et al., 1988) from aortic trimmings of donor hearts at the University of California, Los Angeles (UCLA). All HAECs were de-identified and exempt from consideration as human subjects research by institutional regulatory boards at UC San Diego and The University of Arizona. Cells were grown in culture in M-199 (ThermoFisher Scientific, Waltham, MA, MT-10–060-CV) supplemented with 1.2% sodium pyruvate (ThermoFisher Scientific, Catalog# 11360070), 1% 100X Pen Strep Glutamine (ThermoFisher Scientific Cat# 10378016), 20% fetal bovine serum (FBS, GE Healthcare, Hyclone, Pittsburgh, PA), 1.6% Endothelial Cell Growth Serum (Corning, Corning, NY, Product #356006), 1.6% heparin, and 10 μL/50 mL Amphotericin B (ThermoFisher Scientific #15290018). Cells were grown to 90% confluence in either 10-cm or 15-cm plates, and used primarily at passages 6 to 10. Cells were then treated with M-199 containing 1% FBS (control) or additionally containing either 40 μg/mL Ox-PAPC, 2 ng/mL human recombinant TNFα, or 10 ng/mL human recombinant IL1β (cytokines from R&D Systems, Minneapolis, MN).

Small-interfering RNAs and qPCR Primers

Knockdown of ERG, CEBPD, and IRF1 were performed using 1 nM siRNA oligonucleotides in Opti-MEM (ThermoFisher Scientific) with Lipofectamine 2000 (ThermoFisher Scientific). Transfections were performed in serum-free media for 4 hr, then cells were grown in full growth media for 48 hr. All siRNAs and qPCR primers used in this study are listed in Supplementary file 3.

RNA-seq

HAECs were resuspended in RNA Lysis Buffer and RNA was extracted from cells using the Quick-RNA Micro Prep kit from ZymoResearch (Irvine, CA, #R1051), including optional DNase I treatment. mRNA was selected through poly-A isolation using Oligo d(T)25 beads (New England BioLabs, Ipswich, MA, #S1419S). Selected RNA was fragmented, followed by single strand cDNA synthesis using a SuperScript III First-Strand Synthesis System (ThermoFisher Scientific # 18080051), followed by second strand synthesis using DNA Polymerase I (Qiagen/Enzymatics, Beverly, MA, #P7050L). dsDNA ends were repaired with T4 DNA Polymerase (Enzymatics #P7080L). Barcode adapters (BIOO Scientific NEXTflex, Austin, TX, #514104) were ligated onto the ends of sequences using T4 DNA Ligase (Enzymatics #L-6030-HC-L) and samples were treated with Uracil DNA Glycosylase (UDG) (Enzymatics #G5010L). Libraries were then amplified by PCR (Phusion Hot Start II, ThermoFisher Scientific, #F549L) and purified (Zymo #D5205) for high-throughput sequencing.

Chromatin immunoprecipitation sequencing (ChIP-seq)

ChIP-seq was performed as previously described (Gosselin et al., 2014). Briefly, HAECs were fixed at room temperature with 1% paraformaldehyde in PBS for 10 min, and then quenched with glycine. ChIPs for p65, JUN, and JUNB were performed from chromatin cross-linked by 2 nM Disuccinimidyl Glutarate Crosslinker (DSG) (ProteoChem, Hurricane, UT, #c1104) in PBS for 30 min followed by 1% paraformaldehyde in PBS for 15 min, and then quenched with glycine. Between 2 and 10 million cells were used for each ChIP-seq. Cell lysates were sonicated using a BioRuptor Standard or BioRuptor Pico (Diagenode, Belgium), and then immunoprecipitated using antibodies bound to a 2:1 mixture of Protein A Dynabeads (Invitrogen #10002D) and Protein G Dynabeads (Invitrogen #10004D). Antibodies used included H3K4me2 (EMD Millipore, Billerica, MA, #07–030), H3K27ac (Active Motif, Carlsbad, CA , #39135), CEBPD (Santa Cruz Biotechnology, Dallas, TX, #sc-636X), IRF1 (Santa Cruz #sc-497x), p65 (Santa Cruz #sc-372X), NRF2 (Santa Cruz #sc-1694X), JUN (Santa Cruz #sc-13032X), ERG (Santa Cruz #sc-354X), and JUNB (Santa Cruz #sc-73). Following immunopreciptation, crosslinking was reversed and libraries were prepared using the same method described for RNA-seq beginning with dsDNA end repair and excluding UDG. For each sample condition, an input library was also created using an aliquot of sonicated cell lysate that had not undergone immunoprecipitation. These samples were sequenced as below and used to normalize ChIP-seq results.

Transposase-accessible chromatin using sequencing (ATAC-seq)

ATAC-seq was performed on 50,000 HAEC nuclei according to the original published protocol (Buenrostro et al., 2013) with the exception of size selection (125–175 base pairs on TBE gel) prior to sequencing to enrich for enhancer elements.

Sequencing data samples, mapping, and normalization

Libraries were sequenced on an Illumina HiSeq 4000 according to manufacturer’s specifications at the University California San Diego and at the University of Chicago. Public data was downloaded from public repositories and processed exactly as new data in this study (see below). Reads from ChIP-seq and ATAC-seq were mapped to the hg19 build of the human genome with Bowtie2 (Langmead and Salzberg, 2012) and RNA-seq reads were mapped with STAR (Dobin et al., 2013). For ATAC-seq, reads mapping to mitochondrial DNA were discarded. Mapped reads were organized into HOMER’s preferred data structure called Tag Directories using the ‘makeTagDirectory’ command.

ATAC-seq and RNA-seq experiments that measured accessibility and expression in different treatment stimuli (control, oxPAPC, IL1β, and TNFα) were all conducted at least twice. ERG knockdown was performed in three different HAEC donors. CEBPD and IRF knockdowns were performed in one HAEC donor. The Benjamini-Hochberg false discovery rate (FDR) method was used to correct for all multiple testing in this study. No explicit power analysis was used to compute sample size. Instead, genome-wide features such as enhancer elements served as the replicates for motif analysis and duplicate ChIP-seq experiments confirmed that binding peaks were consistent. With the exception of CEBPD and JUNB that were performed one time, all ChIP-seq experiments were performed in biological duplicate, meaning separate cell expansion and collection.

RNA-seq analysis

For quantification, RNA-seq was normalized using Reads Per Kilobase of transcript per Million mapped reads (RPKM) procedure in HOMER (Heinz et al., 2010). RNA sequencing tags were only considered when they mapped to the same DNA strand as indicated by RefSeq annotation. Further, only tags in exons of genes were incorporated as to remove bias created by variable intron sizes. Together, RNA quantification was achieved using the HOMER command ‘analyzeRepeats rna –strand + -count exons –rpkm’. In this study, RPKM is synonymous with FPKM (Fragments per kilobase mapped). Statistically significant differential expression from RNA-seq experiments was determined first by unnormalized counts (with analyzeRepeats -noadj) followed by statistical testing in DESeq (Anders and Huber, 2010) and a restriction to a 5% False Discovery Rate. Pathway enrichment analysis was performed using Ingenuity Pathway Analysis software (Qiagen).

Peak calling

ChIP-seq and ATAC peaks were identified using un-immunoprecipitated chromain, or ‘input’, as a negative control. Inputs from the corresponding crosslinking condition were used for each ChIP. No input was used to call ATAC peaks. Peaks were identified in HOMER with the findPeaks program according to the data type. Transcription factor peaks were called using the ‘findPeaks -style factor –size 200’ option and histone peaks called using the ‘findPeaks -style histone’ option. ATAC-seq peaks were called using findPeaks with ‘-L 8 F 8 -style histone -size 75 -minDist 75 -minTagThreshold 6’ options. Differential peaks between experiments were determined using the ‘getDifferentialPeaks program with default parameters’ (normalized tag count difference >4 fold and poisson enrichment p-value<0.0001).

Peak merging was performed in HOMER using the ‘mergePeaks’ program. For enhancers defined in Figure 1a, the ATAC-seq peak file was merged with peak regions defined in H3K4me2 and H3K27ac ChIP-seq experiments using the option ‘-d given’ that requires overlap of genomic coordinates between the three peak files. Because the ATAC set was listed first, enhancers were centered on the center of the accessible region. Distal peaks throughput the study were defined as being at least three kilobases from the transcriptional start site of a gene (RefSeq hg19 definitions) using HOMER’s ‘getDistalPeaks.pl’ command.

De novo enhancers and super-enhancers

De novo enhancers were defined as loci that gained ATAC-seq and H3K27ac ChIP-seq in the same cell simulation. Individual gained peak sets was determined by ‘getDifferentialPeaks’, explained above, and significant sets were intersected using ‘mergePeaks’ where the ATAC-seq gained peaks were listed first to maintain centering on accessibility.

Super-enhancers were defined in HOMER from H3K27ac ChIP-seq experiments and corresponding inputs using ‘findPeaks -style super -L 0’. This procedure follows the same logic as the original definition proposed by the Young laboratory (Whyte et al., 2013). Briefly, the implementation in HOMER identifies ‘typical’ ChIP-seq peaks, stitches proximal enhancers together, ranks the resulting enhancers by normalized tag counts over input, and thresholds enhancers above a flex-point (slope >1) as super-enhancers.

Motif enrichment and distance analysis

Motif enrichment analysis was performed on peak sets using HOMER’s ‘findMotifsGenome.pl’ program. As a background control, this program selects a set of sequences from the same genome build that are matched in size and GC content to the peak set of interest. In Figure 5b, to enable enrichment comparison for motifs across multiple peak sets, we used HOMER’s ‘findMotifsGenome.pl –mknown <motifs>’ option iteratively across each de novo enhancer set. For each motif analysis, the amount of sequence analyzed depended on the data type and is indicated in the text.

Analysis of motif distances in Figure 3a–c was performed in peak regions identified by ChIP-seq (e.g. ERG-bound in 3a; AP1 bound in 3b; GATA2-bound in 3c). Each peak set was centered on the highest-scoring TF motif that matched the factor immuno-precipitated, so that 0 bp on the x-axis was the beginning of the likely bound motif (e.g. the ETS motif for ERG-bound in 3a). Peaks lacking consensus motifs for the respective factors were excluded from this analysis. Next, the frequency of the other motifs queried (e.g. AP1, GATA and SOX motifs in 3a) were calculated in HOMER by ‘annotatePeaks.pl –hist –m’ options and the frequency of the additional motifs tested in the vicinity were plotted to show positional relationships among motifs and surrounding genomic sequence.

Whole-genome sequencing, motif mutation analysis and allele-specific factor binding

For the HAEC line sequenced for analysis as in Figure 3d, genomic DNA was prepared for paired-end WGS on Illumina HiSeqX (Illumina, CA) by Novogene (Sacramento, CA) according to manufacturer specifications. More than 80 million raw sequencing reads and a raw depth of 41X were obtained. Data were mapped to the hg19 genome using Burrows-Wheeler Aligner (Li and Durbin, 2009). Single nucleotide variants (SNVs) were identified using GATK (DePristo et al., 2011). SNVs were then compared to 1000 Genomes reference populations (Auton et al., 2015) and the HAEC donor clustered with the Mexican American reference (MXL) population by multi-dimensional scaling analysis in PLINK (Purcell et al., 2007). SNVs were then phased according to the 1000 Genomes MXL reference population in BEAGLE (Browning and Browning, 2016; Browning and Browning, 2007).

Allele-specific binding of JUN and NFκB/p65 was quantified using the WASP pipeline (van de Geijn et al., 2015). To avoid mapping bias caused by alternate alleles, reads that mapped discordantly to reference and alternate alleles at heterozygous SNPs were discarded. Next, sequencing reads from ChIP experiments were aligned to each haplotype at polymorphic loci and summed within ChIP-seq peaks for the respective factor.

Allele-specific motif mutations were identified with the same method as reported previously (Heinz et al., 2013; Gosselin et al., 2014) with slight modifications. In summary, reference genome sequence (hg19) was extracted at peaks of interest and intersected with SNVs in the HAEC donor of interest. Phased variant data was pulled for each homologous chromosome and alleles were inserted in turn to the genome sequence. Motifs were located by alignment to position weight matrices of ETS, AP1, and κB motifs in HOMER (Heinz et al., 2010). Motif mutations were defined as instances where motifs were only identified on one of the two homologous chromosomes.

For visualization, plots as in Figure 3d show the relationship between allele-specific TF binding and mutations in motifs of interest within the TF peaks containing at least one variant. TF binding was transformed to the log2 scale and peaks containing 0 reads on one of the two alleles were scaled to ±7.5 (or else the ratio would be ±Infinity). Boxplots showing the distribution of read ratios in each category (no mutations, mutation on allele 1, mutation on allele 2) exclude peaks with 0 counts on either chromosomal pair. This is a conservative approach, as many peaks have zero reads on the mutated chromosome.

Data visualization

Heatmap-syle histograms of sequencing tags (e.g. Figure 1a), were generated using HOMER’s ‘annotatePeaks.pl -ghist’ option and plotted in R using the ‘heatmap.2()’ function of the gplots library. Cumulative histograms of tag frequencies by position to peak center (e.g. Figure 1c), were generated using HOMER’s ‘annotatePeaks.pl –hist’ option. Standard heatmaps (e.g. Figure 2b) were plotted in R using ‘heatmap.2()’ with default clustering parameters (Hierarchical clustering and Euclidian distance).

Enhancer overlap with GWAS data

Coronary artery disease, hypertension, and related traits were downloaded from the NHGRI-EBI Catalog of published genome-wide association studies (Welter et al., 2014). SNPs in linkage disequilibrium (LD) to GWAS association traits were calculated when r2 >0.8 in PLINK (Purcell et al., 2007) according to the European reference population of the 1000 Genomes Project (Auton et al., 2015). HAEC enhancers defined in Figure 1a were overlapped by physical position (hg19 genome build). The studies reporting associations that overlapped EC enhancers are as follows: (Nikpay et al., 2015; Samani et al., 2007; Coronary Artery Disease (C4D) Genetics Consortium, 2011; Schunkert et al., 2011; Wild et al., 2011; Slavin et al., 2011; Aouizerat et al., 2011; Mehta, 2011; Davies et al., 2012; Wojczynski et al., 2013; Dichgans et al., 2014; Kertai et al., 2015).

Public datasets

With the exception of data used to generate Figure 1—figure supplement 3 (see Supplementary file 1 for list), raw sequencing data was downloaded from Gene Expression Omnibus (GEO) as short read archive files and converted to fastq files using ‘fastq-dump’. Fastq files were mapped to the hg19 genome build in Bowtie2 and processed according to methods outlined for corresponding data types above. The following publicly available datasets were analyzed: GSE52642 (HAEC GRO-seq), GSE20060 (HAEC microarrays), GSE41166 (ETS-1 HUVEC ChIP-seq), and GSE31477 (GATA2 HUVEC ChIP-seq).

Data generated in this study is available under GEO accession GSE89970 and WGS via NCBI SRA Archive, BioProject PRJNA381088.

Cell lines

The HAECs used in this study were primary isolates of ECs from aortic trimmings collected from trimmings of donor hearts transplanted through the UCLA heart transplant program as previously described (Navab et al., 1988). No transformations were performed and they were used at low passage (p6-p10).

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
    A comprehensive 1,000 Genomes-based genome-wide association meta-analysis of coronary artery disease
    1. M Nikpay
    2. A Goel
    3. HH Won
    4. LM Hall
    5. C Willenborg
    6. S Kanoni
    7. D Saleheen
    8. T Kyriakou
    9. CP Nelson
    10. JC Hopewell
    11. TR Webb
    12. L Zeng
    13. A Dehghan
    14. M Alver
    15. SM Armasu
    16. K Auro
    17. A Bjonnes
    18. DI Chasman
    19. S Chen
    20. I Ford
    21. N Franceschini
    22. C Gieger
    23. C Grace
    24. S Gustafsson
    25. J Huang
    26. SJ Hwang
    27. YK Kim
    28. ME Kleber
    29. KW Lau
    30. X Lu
    31. Y Lu
    32. LP Lyytikäinen
    33. E Mihailov
    34. AC Morrison
    35. N Pervjakova
    36. L Qu
    37. LM Rose
    38. E Salfati
    39. R Saxena
    40. M Scholz
    41. AV Smith
    42. E Tikkanen
    43. A Uitterlinden
    44. X Yang
    45. W Zhang
    46. W Zhao
    47. M de Andrade
    48. PS de Vries
    49. NR van Zuydam
    50. SS Anand
    51. L Bertram
    52. F Beutner
    53. G Dedoussis
    54. P Frossard
    55. D Gauguier
    56. AH Goodall
    57. O Gottesman
    58. M Haber
    59. BG Han
    60. J Huang
    61. S Jalilzadeh
    62. T Kessler
    63. IR König
    64. L Lannfelt
    65. W Lieb
    66. L Lind
    67. CM Lindgren
    68. ML Lokki
    69. PK Magnusson
    70. NH Mallick
    71. N Mehra
    72. T Meitinger
    73. FU Memon
    74. AP Morris
    75. MS Nieminen
    76. NL Pedersen
    77. A Peters
    78. LS Rallidis
    79. A Rasheed
    80. M Samuel
    81. SH Shah
    82. J Sinisalo
    83. KE Stirrups
    84. S Trompet
    85. L Wang
    86. KS Zaman
    87. D Ardissino
    88. E Boerwinkle
    89. IB Borecki
    90. EP Bottinger
    91. JE Buring
    92. JC Chambers
    93. R Collins
    94. LA Cupples
    95. J Danesh
    96. I Demuth
    97. R Elosua
    98. SE Epstein
    99. T Esko
    100. MF Feitosa
    101. OH Franco
    102. MG Franzosi
    103. CB Granger
    104. D Gu
    105. V Gudnason
    106. AS Hall
    107. A Hamsten
    108. TB Harris
    109. SL Hazen
    110. C Hengstenberg
    111. A Hofman
    112. E Ingelsson
    113. C Iribarren
    114. JW Jukema
    115. PJ Karhunen
    116. BJ Kim
    117. JS Kooner
    118. IJ Kullo
    119. T Lehtimäki
    120. RJ Loos
    121. O Melander
    122. A Metspalu
    123. W März
    124. CN Palmer
    125. M Perola
    126. T Quertermous
    127. DJ Rader
    128. PM Ridker
    129. S Ripatti
    130. R Roberts
    131. V Salomaa
    132. DK Sanghera
    133. SM Schwartz
    134. U Seedorf
    135. AF Stewart
    136. DJ Stott
    137. J Thiery
    138. PA Zalloua
    139. CJ O'Donnell
    140. MP Reilly
    141. TL Assimes
    142. JR Thompson
    143. J Erdmann
    144. R Clarke
    145. H Watkins
    146. S Kathiresan
    147. R McPherson
    148. P Deloukas
    149. H Schunkert
    150. NJ Samani
    151. M Farrall
    152. CARDIoGRAMplusC4D Consortium
    (2015)
    Nature Genetics 47:1121–1130.
    https://doi.org/10.1038/ng.3396
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
    Large-scale association analysis identifies 13 new susceptibility loci for coronary artery disease
    1. H Schunkert
    2. IR König
    3. S Kathiresan
    4. MP Reilly
    5. TL Assimes
    6. H Holm
    7. M Preuss
    8. AF Stewart
    9. M Barbalic
    10. C Gieger
    11. D Absher
    12. Z Aherrahrou
    13. H Allayee
    14. D Altshuler
    15. SS Anand
    16. K Andersen
    17. JL Anderson
    18. D Ardissino
    19. SG Ball
    20. AJ Balmforth
    21. TA Barnes
    22. DM Becker
    23. LC Becker
    24. K Berger
    25. JC Bis
    26. SM Boekholdt
    27. E Boerwinkle
    28. PS Braund
    29. MJ Brown
    30. MS Burnett
    31. I Buysschaert
    32. JF Carlquist
    33. L Chen
    34. S Cichon
    35. V Codd
    36. RW Davies
    37. G Dedoussis
    38. A Dehghan
    39. S Demissie
    40. JM Devaney
    41. P Diemert
    42. R Do
    43. A Doering
    44. S Eifert
    45. NE Mokhtari
    46. SG Ellis
    47. R Elosua
    48. JC Engert
    49. SE Epstein
    50. U de Faire
    51. M Fischer
    52. AR Folsom
    53. J Freyer
    54. B Gigante
    55. D Girelli
    56. S Gretarsdottir
    57. V Gudnason
    58. JR Gulcher
    59. E Halperin
    60. N Hammond
    61. SL Hazen
    62. A Hofman
    63. BD Horne
    64. T Illig
    65. C Iribarren
    66. GT Jones
    67. JW Jukema
    68. MA Kaiser
    69. LM Kaplan
    70. JJ Kastelein
    71. KT Khaw
    72. JW Knowles
    73. G Kolovou
    74. A Kong
    75. R Laaksonen
    76. D Lambrechts
    77. K Leander
    78. G Lettre
    79. M Li
    80. W Lieb
    81. C Loley
    82. AJ Lotery
    83. PM Mannucci
    84. S Maouche
    85. N Martinelli
    86. PP McKeown
    87. C Meisinger
    88. T Meitinger
    89. O Melander
    90. PA Merlini
    91. V Mooser
    92. T Morgan
    93. TW Mühleisen
    94. JB Muhlestein
    95. T Münzel
    96. K Musunuru
    97. J Nahrstaedt
    98. CP Nelson
    99. MM Nöthen
    100. O Olivieri
    101. RS Patel
    102. CC Patterson
    103. A Peters
    104. F Peyvandi
    105. L Qu
    106. AA Quyyumi
    107. DJ Rader
    108. LS Rallidis
    109. C Rice
    110. FR Rosendaal
    111. D Rubin
    112. V Salomaa
    113. ML Sampietro
    114. MS Sandhu
    115. E Schadt
    116. A Schäfer
    117. A Schillert
    118. S Schreiber
    119. J Schrezenmeir
    120. SM Schwartz
    121. DS Siscovick
    122. M Sivananthan
    123. S Sivapalaratnam
    124. A Smith
    125. TB Smith
    126. JD Snoep
    127. N Soranzo
    128. JA Spertus
    129. K Stark
    130. K Stirrups
    131. M Stoll
    132. WH Tang
    133. S Tennstedt
    134. G Thorgeirsson
    135. G Thorleifsson
    136. M Tomaszewski
    137. AG Uitterlinden
    138. AM van Rij
    139. BF Voight
    140. NJ Wareham
    141. GA Wells
    142. HE Wichmann
    143. PS Wild
    144. C Willenborg
    145. JC Witteman
    146. BJ Wright
    147. S Ye
    148. T Zeller
    149. A Ziegler
    150. F Cambien
    151. AH Goodall
    152. LA Cupples
    153. T Quertermous
    154. W März
    155. C Hengstenberg
    156. S Blankenberg
    157. WH Ouwehand
    158. AS Hall
    159. P Deloukas
    160. JR Thompson
    161. K Stefansson
    162. R Roberts
    163. U Thorsteinsdottir
    164. CJ O'Donnell
    165. R McPherson
    166. J Erdmann
    167. NJ Samani
    168. Cardiogenics
    169. CARDIoGRAM Consortium
    (2011)
    Nature Genetics 43:333–338.
    https://doi.org/10.1038/ng.784
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71

Decision letter

  1. Manuel Garber
    Reviewing Editor; University of Massachusetts Medical School

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Transcriptional networks specifying homeostatic and inflammatory expression programs in human aortic endothelial cells" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Harry Dietz as the Senior Editor. The following individual involved in review of your submission has agreed to reveal her identity: Lynn Butler (Reviewer #2).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Hogan et al. performed an extensive analyses of transcriptional regulation in human aortic endothelial cells (HAECs). They authors provide a comprehensive map of the gene regulatory landscape of HAEC cells under basal (part1) and inflammatory conditions (part2). Using a combination of enhancer annotation, in-silico motif analysis and transcription factor ChIP the authors find that a Jun (an AP1 factor) and ERG (an ETS factor) cooperate to establish much of the cell enhancer landscape. In the second part, the authors study the transcriptional network in response to inflammatory stimulus.

Essential revisions:

1) It is unclear why the focus is only on JUNB. The reviewers felt there the authors should explain why JUND is not pursued further.

2) Two reviewers, one of them from prior experience, are worried that some of the observed effects are due to the transfection itself, or other off-target effects. The authors should carefully establish this before they conclude that ERG, or any other siRNA knock-down is the cause of the (observed) inflammatory response. The reviewers strongly suggest a combination of appropriate controls and replication with different siRNAs to ensure that the result is not due to secondary effects not related to the result of the ERG and all other knock-down experiments.

3) The collaborative nature of JUNB and ERG function would be strengthen by showing the effect of loss of function of JUNB.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Transcriptional networks specifying homeostatic and inflammatory expression programs in human aortic endothelial cells" for further consideration at eLife. Your revised article has been favorably evaluated by Harry Dietz (Senior editor) and a Reviewing editor.

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

There is one remaining point you may want to discuss. While the AP1 factors are highly expressed, the ETS factors are not. Figure 4—figure supplement 1 shows ERG expression to be less than 2.5 FPKM. Your data shows ERG and JUN co-binding on most regions. Is the protein of ERG that stable? Are the stoichiometry relationships reasonable? Clearly, at their mRNA levels ERG and JUN are at very different, is there evidence that ERG is a more stable protein? A related point: In Figure 3—figure supplement 2, Two of the three examples (ERK and TEK) have very small fold changes in siERG. What about other genes that show stronger effects?

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled "Transcriptional networks specifying homeostatic and inflammatory expression programs in human aortic endothelial cells" for further consideration at eLife. Your revised article has been favorably evaluated by Harry Dietz (Senior editor), a Reviewing editor, and two reviewers.

The manuscript has been improved but there is one remaining issue that need to be addressed before acceptance, as outlined below:

While the methods for quantification of RNA-Seq data have been clearly improved since the original submission, there was consensus that these methods are not adequately described in the revised manuscript. Please expand and refine this description in the Materials and methods section.

https://doi.org/10.7554/eLife.22536.041

Author response

Essential revisions:

1) It is unclear why the focus is only on JUNB. The reviewers felt there the authors should explain why JUND is not pursued further.

We agree that JUND is likely to be one of the important AP-1 family members involved in regulating endothelial transcription. This is evident from motif analysis of the 17 thousand HAEC enhancers that revealed a prominent role in AP1 transcription factors (Figure 1A, B, Figure 1—figure supplement 1A), high levels of JUND mRNA in HAECs (Figure 1D, Figure 1—figure supplement 1B) and the presence of a super-enhancer at the JUND locus (Figure 1—figure supplement 2, middle panel). We proceeded without JUND binding in the original manuscript because we had information about the binding of the other two AP1 factors (JUN and JUNB) and reasoned that with the redundancy in binding between factors and the heterodimeric binding property of AP1 factors that these two factors would be sufficient to provide considerable insight into functional roles of AP-1 as a whole. Consistent with this, we observe that the cistromes for JUN and JUNB are nearly identical (Figure 2—figure supplement 1). We did, however, fail to explicitly state our reasoning – an oversight that has now been corrected and is evident in the following sections of the revised manuscript:

Results section: “Tens of genes encode proteins of the AP1, ETS, SOX and GATA TF families. Within each family different members share nearly identical DNA binding domains and thus bind the same motif. In addition, AP1 protein members bind AP1 motifs as homo- and hetero-dimers. These sources of redundancy make it challenging to identify the functional family member(s) without additional information.”

Results section: “AP1 family members JUND, JUN, and JUNB were also in the top 4%. RNA-seq from other HAEC donors and replicate samples confirmed these findings (Figure 1—figure supplement 1). Next, we defined SEs using H3K27ac ChIP-seq data and found that, among others, the genetic loci for ERG (an ETS member) as well as AP1 members JUN, JUND, and JUNB harbored SEs (Figure 1E, F, Figure 1—figure supplement 2). Taken together, these data suggest that while multiple TFs from each family probably bind HAEC enhancers, that JUN, JUNB, JUND and ERG likely serve prominent roles.”

Results section: “The JUND cistrome would also be informative in these studies; however, we proceeded with JUN and JUNB because the heterodimeric binding of AP1 factors makes it likely that JUN and JUNB profiles encompass a major portion of the overall AP1 landscape. Consistent with this, the JUN and JUNB cistromes are highly concordant (Figure 2—figure supplement 2).”

2) Two reviewers, one of them from prior experience, are worried that some of the observed effects are due to the transfection itself, or other off-target effects. The authors should carefully establish this before they conclude that ERG, or any other siRNA knock-down is the cause of the (observed) inflammatory response. The reviewers strongly suggest a combination of appropriate controls and replication with different siRNAs to ensure that the result is not due to secondary effects not related to the result of the ERG and all other knock-down experiments.

This point is well taken. Additional experiments were performed (refer to Figure 4—figure supplement 2) by each of the 4 individual siRNAs that comprised the pooled siERG mix used in the original submission (siERGs #1-4) as well as two additional single ERG-targeting siRNAs (siERG#5-6). In addition, two non-targeting scrambled oligos were used for comparison (siSCR#1-2) along with negative controls including Lipofectamine2000 only, oligo only, and optiMEM only, and full growth media changes only. qPCR for ERG was used to confirm knock-down and 1 of the 6 siERG oligos failed to reduce ERG message (siERG #6). In addition, inflammatory genes LIF, IL8, IL6, F3, IL1B, IL1A, CSF3, and CCL2 were measured across experimental conditions in three biological replicates. Data show that complete transfection complexes, e.g., siSCRs, slightly increased inflammatory messages for some genes. For example, IL8 levels were ~2.5 fold induced by siSCR compared to the non-transfected controls. This was in contrast to the robust up-regulation of inflammatory gene expression that was reproducibly observed across individual oligos targeting ERG. siERG6, which failed to reduce ERG levels, served as a control in that inflammatory expression was not induced for this condition. In addition, we measured 4 endothelial-specific genes in this experiment: VWF, PECAM1, EGFL7 AND NOS3. VWF was reproducibly significantly reduced upon ERG knock-down and levels of the other genes was also reduced – however, less dramatically than VWF. Our conclusion from these experiments is that ERG knock-down was the cause of the pro-inflammatory gene expression profile we (and others) have observed in endothelial cells. This information has now been added to the revised manuscript as follows:

Results section: “To ensure that the inflammatory gene profile elicited by ERG knock-down was not a consequence of transfection itself or off-target effects, the profile resulting from six individual siERG oligos was measured along with two non-targeting scrambled siRNA controls and non-transfection controls (Figure 4—figure supplement 2). These data were reproducible and consistent with ablated ERG expression as the cause of pro-inflammatory expression profiles. Together, these data suggest that ERG normally functions to maintain EC-specific gene functions such as development and proliferation while at the same time suppressing inflammatory gene expression.”

3) The collaborative nature of JUNB and ERG function would be strengthen by showing the effect of loss of function of JUNB.

We agree that an implication of collaborative TF binding is that loss-of-function of a participating TF would decrease the ability of the other factor to bind target elements. The expectation to observe this effect genome-wide is complicated, however, by other TFs in the same families that can bind the same DNA-binding motif and provide a surrogate collaborative role. Considering that we observe multiple AP1 and ETS members with high expression in HAECs, we decided to undertake an alternative approach to test for collaborative interactions.

This approach stems from our previous work in genetically diverse inbred mouse strains where millions of naturally occurring SNPs between strains were considered nature’s mutagenesis experiment (Heinz, Romanoski et al. Nature 2013; Gosselin et al., 2014). The current study is the first to implement this approach in human cells. The approach is as follows: SNPs are identified across the genome within a single individual. We achieved this using whole-genome sequencing (40X average depth) in one aortic endothelial cell donor. Next, we identified loci where the alleles of the SNPs mutated a canonical DNA binding motif for a TF of interest (AP1, ETS, or NFkB) only on one homologous chromosome. These are the potential functional SNPs that are used to test for effects on TF binding. Next, we developed an analytical pipeline using a combination of public tools (WASP, HOMER, R) to quantify allele-specific TF binding for JUN and NFkB/p65 that takes into account personal genome sequence and avoids caveats of mapping bias (described in the Materials and methods section).

Results show that allele-specific mutations in AP-1 motifs genome-wide (n = 237) are coincident with reduced JUN binding in the expected allele-specific manner (Figure 3D, Figure 3—figure supplement 1, p = 5e-10). This demonstrates that motif mutations have an effect on binding of the respective factor. Perhaps more interestingly, we found that ETS motif mutations were coincident with reduced JUN binding within 100 base pairs genome-wide (n = 150, p=3e-3). This provides further evidence that ETS and AP1 factors collaborate in binding to endothelial enhancers. We extend this approach to the later portion of our manuscript where we find that altered binding of NFkB/p65 associated with motif mutations in the kB motif itself as well as AP1 motifs and ET motifs (Figure 5—figure supplement 4). These data support a hierarchical role where AP1 and ETS factors prime the chromatin landscape for binding of signal-dependent factors like NFkB. These data and analytical details have been added to the manuscript as follows:

Results:

“Allele-specific binding to chromosomes lacking motif mutations supports collaborative binding between AP1 and ETS factors […] Taken together, these data support a collaborative relationship between AP1 and ETS factors at endothelial enhancers.”

Results: “In addition, we analyzed allele-specific NFkB binding as a function of motif mutations in the kB motif itself as well as mutations in the AP1 and ETS motif. […]Furthermore, these data demonstrate that allele-specific binding in heterogeneous human cells is a useful means to reveal collaborative interactions between transcription factors.”

Materials and methods:

“Whole-Genome Sequencing (WGS), Motif Mutation Analysis and Allele-Specific Factor Binding For the HAEC line sequenced for analysis as in Figure 3D, genomic DNA was prepared for paired-end WGS on Illumina HiSeqX (Illumina, CA) by Novogene (Sacramento, CA) according to manufacturer specifications[…]This is a conservative approach, as many peaks have zero reads on the mutated chromosome.”

[Editors' note: further revisions were requested prior to acceptance, as described below.]

The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

There is one remaining point you may want to discuss. While the AP1 factors are highly expressed, the ETS factors are not. Figure 4—figure supplement 1 shows ERG expression to be less than 2.5 FPKM. Your data shows ERG and JUN co-binding on most regions. Is the protein of ERG that stable? Are the stoichiometry relationships reasonable?

When quantifying gene expression from RNA-seq data as all tags between the transcription start site to the termination site, the expression levels between ETS and AP-1 factors are different. This is largely because ETS factors tend to have multiple large introns (e.g., ERG in Figure 1E) and AP-1 factors lack introns (e.g., JUN in Figure 1F). In both cases, however, the tag counts over the exons is comparable. This is further shown by comparison of microarray data measuring ETS and AP-1 factors (Figure 2A), which shows similar levels of expression between ERG and JUN.

To enable more accurate comparisons between genes with different numbers of introns, we have now updated the definition by which RNA-seq data is counted (now restricted only to exons). This revision is reflected in Figure 1D, Figure 1—figure supplement 1B, Figure 4C, and Figure 4—figure supplement 1. The results remain largely the same with the main exception that ERG is rank ordered higher among transcription factors in general (Figure 1C). The results and conclusions of ERG knock-down experiments remain the same (Figures 4C, and Figure 4—figure supplement 1).

While our study highlights a collaborative role for ERG and JUN specifically, we acknowledge that this combination is not likely the sole modulating combination. In fact, other collaborative interactions involving other AP-1 and ETS factors almost certainly regulate EC gene expression, and the cistromes of these factors and their functional effects remain a target of future investigation.

We have now revised our manuscript to reflect these changes in the following ways:

Discussion:

“While we focus on ERG and JUN in particular, other combinatorial interactions between other ETS and AP-1 members almost certainly play prominent roles in endothelial regulation.”

Clearly, at their mRNA levels ERG and JUN are at very different, is there evidence that ERG is a more stable protein? A related point: In Figure 3—figure supplement 2, Two of the three examples (ERK and TEK) have very small fold changes in siERG. What about other genes that show stronger effects?

We believe this comment refers to loci shown for TEK and CDH5. PECAM1, what was more significantly affected by siERG, is now included (Figure 3—figure supplement 2) that displays many ERG/JUN co-bound enhancers.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

The manuscript has been improved but there is one remaining issue that need to be addressed before acceptance, as outlined below:

While the methods for quantification of RNA-Seq data have been clearly improved since the original submission, there was consensus that these methods are not adequately described in the revised manuscript. Please expand and refine this description in the methods section.

We now present a more detailed description of the RNA-seq analysis that reflects the precise options used to quantify gene expression and determine differential expression. This change can be found in subsection “RNA-seq Analysis” as follows:

“RNA-seq Analysis For quantification, RNA-seq was normalized using Reads Per Kilobase of transcript per Million mapped reads (RPKM) procedure in HOMER[…] Pathway enrichment analysis was performed using Ingenuity Pathway Analysis software (Qiagen).”

https://doi.org/10.7554/eLife.22536.042

Article and author information

Author details

  1. Nicholas T Hogan

    Department of Cellular and Molecular Medicine, University of California, San Diego, San Diego, United States
    Contribution
    NTH, Conceptualization, Investigation, Writing—original draft, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared.
  2. Michael B Whalen

    Department of Cellular and Molecular Medicine, University of Arizona, Tucson, United States
    Contribution
    MBW, Investigation, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared.
  3. Lindsey K Stolze

    Department of Cellular and Molecular Medicine, University of Arizona, Tucson, United States
    Contribution
    LKS, Formal analysis, Investigation, Writing—original draft, Writing—review and editing
    Competing interests
    No competing interests declared.
  4. Nizar K Hadeli

    Department of Cellular and Molecular Medicine, University of Arizona, Tucson, United States
    Contribution
    NKH, Formal analysis, Validation, Investigation
    Competing interests
    No competing interests declared.
  5. Michael T Lam

    Department of Cellular and Molecular Medicine, University of California, San Diego, San Diego, United States
    Contribution
    MTL, Conceptualization, Investigation
    Competing interests
    No competing interests declared.
  6. James R Springstead

    Department of Chemical and Paper Engineering, University of Western Michigan, Kalamazoo, United States
    Contribution
    JRS, Resources, Methodology
    Competing interests
    No competing interests declared.
  7. Christopher K Glass

    Department of Cellular and Molecular Medicine, University of California, San Diego, San Diego, United States
    Contribution
    CKG, Conceptualization, Resources, Writing—review and editing
    Competing interests
    CKG: Reviewing editor, eLife
    ORCID icon 0000-0003-4344-3592
  8. Casey E Romanoski

    Department of Cellular and Molecular Medicine, University of Arizona, Tucson, United States
    Contribution
    CER, Conceptualization, Resources, Data curation, Formal analysis, Supervision, Funding acquisition, Investigation, Visualization, Methodology, Writing—original draft, Project administration, Writing—review and editing
    For correspondence
    cromanoski@email.arizona.edu
    Competing interests
    No competing interests declared.
    ORCID icon 0000-0002-0149-225X

Funding

National Institutes of Health (R00123485)

  • Casey E Romanoski

Sarnoff Cardiovascular Research Foundation (Research Fellowship)

  • Nicholas T Hogan

National Institutes of Health (1R15HL121770-01A1)

  • James R Springstead

National Institutes of Health (DK091183)

  • Christopher K Glass

National Institutes of Health (DK063491)

  • Christopher K Glass

National Institutes of Health (GM085764)

  • Christopher K Glass

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

CER was supported by NIH grant R00123485 and CKG by NIH grants DK091183, DK063491, and GM085764. NH was supported as a Research Fellow by the Sarnoff Cardiovascular Research Foundation. JRS was supported by NIH grant 1R15HL121770-01A1.

Reviewing Editor

  1. Manuel Garber, Reviewing Editor, University of Massachusetts Medical School

Publication history

  1. Received: October 20, 2016
  2. Accepted: May 22, 2017
  3. Version of Record published: June 6, 2017 (version 1)

Copyright

© 2017, Hogan et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 979
    Page views
  • 187
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Scopus, Crossref.

Comments

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Biochemistry
    Martin Steger et al.
    Research Advance Updated