A dynamic and integrated epigenetic program at distal regions orchestrates transcriptional responses to VEGFA

  1. Bing Zhang1
  1. 1Key Laboratory of Systems Biomedicine, Shanghai Center for Systems Biomedicine, Xin Hua Hospital, Shanghai Jiao Tong University, Shanghai 200240, China;
  2. 2Department of Biostatistics and Computational Biology, Dana-Farber Cancer Institute and Harvard T.H. Chan School of Public Health, Boston, Massachusetts 02215, USA;
  3. 3School of Electronic Information and Electrical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China;
  4. 4Department for Biomedical Informatics, Harvard Medical School, Boston, Massachusetts 02115, USA;
  5. 5Department of Cardiology, Boston Children's Hospital, Boston, Massachusetts 02115, USA;
  6. 6Renji-Med Clinical Stem Cell Research Center, Renji Hospital, School of Biomedical Engineering, Shanghai Jiao Tong University, Shanghai 200127, China;
  7. 7Department of Molecular Medicine, University of Texas Health Science Center, San Antonio, Texas 78229, USA;
  8. 8Harvard Stem Cell Institute, Cambridge, Massachusetts 02138, USA

  1. ↵9 These authors are co-first authors and contributed equally to this work.

  • Corresponding authors: gcyuan{at}jimmy.harvard.edu, wpu{at}pulab.org, bingzhang{at}sjtu.edu.cn
    •  
    Next Section

    Abstract

    Cell behaviors are dictated by epigenetic and transcriptional programs. Little is known about how extracellular stimuli modulate these programs to reshape gene expression and control cell behavioral responses. Here, we interrogated the epigenetic and transcriptional response of endothelial cells to VEGFA treatment and found rapid chromatin changes that mediate broad transcriptomic alterations. VEGFA-responsive genes were associated with active promoters, but changes in promoter histone marks were not tightly linked to gene expression changes. VEGFA altered transcription factor occupancy and the distal epigenetic landscape, which profoundly contributed to VEGFA-dependent changes in gene expression. Integration of gene expression, dynamic enhancer, and transcription factor occupancy changes induced by VEGFA yielded a VEGFA-regulated transcriptional regulatory network, which revealed that the small MAF transcription factors are master regulators of the VEGFA transcriptional program and angiogenesis. Collectively these results revealed that extracellular stimuli rapidly reconfigure the chromatin landscape to coordinately regulate biological responses.

    Divergent gene programs control distinct cell identities and biological functions. Environmental signals guide cell behavior by modulating gene expression, but the transcriptional and epigenetic mechanisms that underlie rapid, signal-induced gene expression changes are incompletely understood. As an extracellular growth factor that controls almost every step of angiogenesis, vascular endothelial growth factor A (VEGFA) exemplifies the powerful effect of environmental cues on cellular gene expression and function (Leung et al. 1989). Although VEGFA-induced angiogenesis is essential for vertebrate organ development and tissue repair, and abnormalities of angiogenesis and VEGFA signaling are linked to diseases with high morbidity and mortality like myocardial infarction, stroke, and macular degeneration, the gene program temporally controlled by VEGFA and its transcriptional regulatory mechanisms are incompletely understood (Carmeliet 2005).

    Diverse combinations of histone modifications generate an epigenetic code that governs gene activation and repression (Strahl and Allis 2000; Hake et al. 2004). This code is established by epigenetic enzymes that read and write histone modifications, and by sequence-specific transcription factors (TFs), which recruit epigenetic enzymes to specific genomic loci. Targeted studies over the past decade have demonstrated essential roles of histone modifications, epigenetic enzymes, and TFs in regulating angiogenesis in development and disease. For example, EP300 and CBP, acetyl-transferases that deposit activating acetyl-marks on histone residues, including lysine residues 4, 9, and 27 of histone H3 (H3K4ac, H3K9ac, and H3K27ac), are essential to vascular development and VEGFA responses (Yao et al. 1998). Their action is counter-balanced by histone deacetylases, including HDAC6, -7, and -9, which likewise are essential for normal angiogenesis (Zhang et al. 2002; Chang et al. 2006; Birdsey et al. 2012). EZH2, the catalytic subunit of polycomb repressive complex 2 (PRC2), represses genes by trimethylating lysine 27 of histone H3 (H3K27me3) and is required for promoting angiogenesis in tumors (Lu et al. 2010). EZH2 is dispensable for developmental angiogenesis (Yu et al. 2017b), pointing out important differences in the epigenetic regulation of these distinct angiogenic programs. A number of TFs, including members of the ETS, GATA, FOX, and SOX TF families, have been shown similarly to have essential roles for angiogenesis in development and disease (De Val and Black 2009). In particular, members of the ETS TF family are key regulators of angiogenesis, often through combinatorial interactions with other TFs, most notably Forkhead family members (De Val and Black 2009). Our recent study showed that one ETS family member, ETS1, broadly regulates endothelial gene expression to promote angiogenesis (Chen et al. 2017).

    Despite these advances in identifying essential roles of histone modifications and TFs in the regulation of angiogenesis, there is a paucity of information about how these factors control the responses of endothelial cells to extracellular signals, which underlies the intricate process of angiogenesis. A major barrier has been the lack of a global map of the transcriptional and epigenetic landscape of endothelial cells responding to key angiogenic factors, such as VEGFA. In this study, we used multiple genome-wide approaches to unveil the time-dependent effect of VEGFA on the epigenetic and transcriptional landscape of endothelial cells.

    Previous SectionNext Section

    Results

    VEGFA induces a temporal change in transcription

    To identify the genes regulated by VEGFA in endothelial cells, we measured mRNA and lncRNA expression by RNA-seq in human umbilical vein endothelial cells (HUVECs) at 0 (unstimulated), 1, 4, and 12 h after addition of VEGFA. Eight hundred seventy-four mRNAs and 61 lncRNAs were differentially expressed (absolute fold change ≥2 and FDR ≤ 0.1) at 1, 4, or 12 h compared with 0 h (Fig. 1A; Supplemental Tables S1, S2). We validated eight differentially expressed genes (DEGs) by RT-qPCR and found similar dynamic changes to RNA-seq (Supplemental Fig. S1A). Many of these DEGs, such as FOXC1, MEF2C, DLL4, HDAC9, and other genes highlighted in Figure 1A, are known to participate in angiogenesis (Prasitsak et al. 2015; Sacilotto et al. 2016; Wilhelm et al. 2016; Hasan et al. 2017).

    Figure 1.
    View larger version:
    Figure 1.

    Dynamic transcriptome changes induced by VEGFA in HUVECs. (A) Heatmap of differentially expressed genes [DEGs; fold change compared with baseline (H0) ≥2, two biological replicates per time point]. HUVECs were treated with VEGFA for 0 (untreated), 1, 4, and 12 h, and gene expression was profiled by mRNA-seq. Genes were separated into seven clusters according to their temporal expression pattern. Gene number, representative genes, and selected enriched Gene Ontology terms for each cluster are listed to the right. (B) Heatmap of differentially expressed TFs. Representative TFs are listed on the right. (C) Volcano plots of expressed lncRNAs. Dots are colored by differential expression compared to baseline: (green) FC ≥ 2 or FC ≤ 0.5 and P-value <0.05; (beige) FC ≥ 2 or FC ≤ 0.5 and P-value ≥0.05; (red) 0.5 < FC < 2 and P-value <0.05. (D) Correlation of expression of lncDEGs to their potential target genes. The curves indicate the distribution of Pearson correlation coefficients (PCCs) between expression of lncDEGs and their potential target genes. The Kolmogorov–Smirnov test P-value comparing the distributions for neighboring target genes and randomly selected genes (excluding neighboring genes) is shown. (E) Expression of lncRNA SNHG15 and its adjacent gene CCM2 during the VEGFA stimulation time course, illustrating a positive correlation between a lncRNA-adjacent gene pair over time.

    The DEGs were grouped into seven clusters (G1-7) according to their temporal expression pattern (Fig. 1A). Genes in distinct functional classes tended to segregate into distinct clusters. More than half (55%) of the genes in clusters G4-5 were TFs and rapidly up-regulated at 1 h (Fig. 1A,B), suggesting that extracellular signals rapidly induce genes that impact cellular behavior by modifying their transcriptional output. Among these TFs were FOXC1, NR4A1, FOSL1, KLF4, and HLX, which are known to regulate angiogenesis (Prasitsak et al. 2015). DUSP1 and DUSP5, which encode phosphatases that antagonize activated MAPK signaling, were also up-regulated in these two clusters. Because VEGFA signal transduction requires RAS-MAPK activation, up-regulation of DUSP establishes a negative feedback circuit that antagonizes VEGFA signal transduction. Genes implicated in vascular remodeling, such as CXCR4, SEMA6D, UNC5B, and ADAMTS9, were predominantly expressed in late-responsive clusters G3, G6, and G7. Moreover, genes implicated in inflammatory responses, including those in the NFκB pathway, were found in clusters G1, G5, and G7, highlighting an under-recognized proinflammatory role of VEGFA.

    Long noncoding RNAs (lncRNAs) have been demonstrated recently to function in many biological processes including angiogenesis (Kumar and Goyal 2017). In our RNA-seq of polyadenylated RNAs, we found 2814 well-annotated lncRNAs expressed in HUVEC cells (FPKM ≥ 1 in at least one time point) (Fig. 1B; Supplemental Table S2). The lncRNAs MALAT1, XIST, and MEG3, which are functionally implicated in tumor angiogenesis (Michalik et al. 2014; Kumar and Goyal 2017; Yu et al. 2017a), were robustly detected in our RNA-seq data. Expression of 61 lncRNAs significantly changed with VEGFA treatment (lncDEG: absolute fold change ≥2 and P-value ≤0.05) (Fig. 1C; Supplemental Fig. S1B). Among them, SNHG8, MIR22HG, MIR503HG, and SNHG15 were previously reported to become differentially expressed in endothelial cells under ischemic conditions, which activate VEGFA signaling (Voellenkle et al. 2016).

    A subset of lncRNAs regulates expression of neighboring protein-coding genes in cis (Sigova et al. 2013). Therefore, we tested the hypothesis that differentially expressed lncRNAs regulate VEGFA responsiveness by altering transcription of neighboring genes. We calculated the expression correlation between lncDEGs and their adjacent genes and found that it was significantly stronger compared to lncDEGs and randomly selected, nonneighboring genes (P-value = 0.00096, Kolmogorov–Smirnov test) (Fig. 1D). Consistent with this result, multiple lncDEGs had dynamic expression patterns that were correlated or anticorrelated to adjacent coding genes with critical roles in angiogenesis. Examples of such lncDEG–DEG pairs are SNHG15–CCM2, LINC01013–CTGF, RP11–1100L3.8-NR4A1, and MIR22HG–WDR81 (Fig. 1E; Supplemental Fig. S1C). Together, these results suggest that VEGFA regulates the transcription of some coding genes by modulating the expression of their adjacent lncRNAs.

    Epigenetic dynamics at promoters do not correlate with VEGFA responsiveness

    Epigenetic modifications of histones play a pivotal role in regulating gene transcription, yet information on the epigenetic chromatin landscape and its relationship to gene expression in angiogenesis is lacking. To gain insights, we used chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) to profile the genome-wide chromatin occupancy of multiple histone modifications in HUVEC cells stimulated by VEGFA. The marks that we profiled are H3K4me1, a marker of active and poised transcriptional enhancers; H3K4me2, a marker of both enhancers and promoters; H3K4me3, a marker of active promoters; H3K27ac, a marker of active enhancers and promoters; H3K27me3, a marker of repressed promoters; and H3K36me3, a marker of actively transcribed genes (Goldberg et al. 2007; Heintzman et al. 2007). These histone marks at hour 0 were well correlated with ENCODE data sets profiled in HUVEC cells, indicating a high quality of our ChIP-seq data (R = 0.55–0.89) (Supplemental Fig. S2). In addition, we profiled the occupancy of RNA polymerase II as a marker of transcriptional activity. The characteristics of our data sets for these seven markers at the four points in the VEGFA stimulation time course (0, 1, 4, and 12 h) are listed in Supplemental Tables S1 and S5.

    We first examined the temporal patterns of these seven chromatin features at promoters, defined as transcription start sites (TSS) ± 1 kb. Using k-means clustering, we identified nine temporal clusters (T1–T9) (Fig. 2A). Clusters T1–T3, T5, and T7 had signatures of active promoters (occupancy by active histone modifications H3K4me2, H3K4me3, and H3K27ac, and by RNAP II), cluster T6 was occupied predominantly by the repressive mark H3K27me3, and cluster T4 was occupied by both activating (H3K4me3) and repressive (H3K27me3) marks, typical of “bivalent” promoters (Bernstein et al. 2006). The remaining two clusters contained either low signal for all markers (T8) or signal limited to H3K36me3 and RNAP II (T9). We first interrogated the relationship between these temporal epigenetic promoter patterns and gene expression by calculating the enrichment of DEG clusters (G1–G7) within the epigenetic promoter clusters (Fig. 2B). DEGs were generally overrepresented in active clusters T1–T5, which shared high H3K4me2 and H3K4me3. Notably, cluster (T4) with both repressive H3K27me3 and active H3K4me3 marks was enriched with small portion of DEGs from all expression clusters.

    Figure 2.
    View larger version:
    Figure 2.

    The epigenetic changes at promoters do not configure VEGFA transcriptional responses. (A) Heatmap showing promoter epigenetic signatures in HUVEC cells stimulated by VEGFA. Promoters were grouped by k-means clustering based on their temporal epigenetic signature for six histone modifications and RNAP II, as determined by ChIP-seq. Heatmap shows mean input subtracted values for promoters in each cluster. (B) Enrichment of DEG clusters (G1–G7) among genes associated with specific temporal patterns of epigenetic marks at the promoter. (T1–T9) Promoters clustered by temporal epigenetic signature. (G1–G7) DEGs clustered by temporal expression pattern (see Fig. 1A). Enrichment P-value was calculated using Fisher's exact test. Blue rectangle highlights the lack of enrichment of T3 genes in DEG clusters G4–G6. (C) Relationship between promoter histone mark signal strength and DEG frequency. Promoters were ranked by H3K27ac, H3K27me3, H3K4me2, or H3K4me3 signal strength and divided into 300 equal sized bins. The number of DEGs within each bin was plotted in gray. Red, orange, purple, and yellow color bars indicate ChIP-seq signal strength within each bin. The blue spline is the trend line of DEG counts. (D) Relationship between promoter histone mark signal variation and DEG frequency. Plot is comparable to C, except that promoters were ranked by the coefficient of variation of H3K27ac or H3K4me3 signal and divided into 300 equal sized bins. The number of DEGs within each bin was plotted in gray. Red, orange, purple, and yellow color bars indicate ChIP-seq signal variation within each bin. The blue spline is the trend line of DEG counts.

    Given the combinations of histone modifications did not show particular association with DEG, we further examined the contribution of both the strength and variation of epigenetic marks at the promoter to VEGFA responsiveness. For H3K27ac, H3K4me3, H3K4me2, and H3K27me3, we plotted the relationship between DEG frequency and signal strength of modified histones averaged over the four time points (Fig. 2C). This revealed a tendency for genes with higher mean promoter occupancy by the three activating histone modifications and with lower mean promoter occupancy by repressive H3K27me3 to have greater DEG frequency. In contrast, we did not observe a relationship between the variation of modified histone signal strength across time points and DEG frequency (Fig. 2D). These findings were further supported by logistic regression, in which we used the mean strength and variation of these four modified histones to develop a predictor of whether or not a gene was differentially expressed. In 10-fold cross-validation, the classifier based on these parameters had moderate predictive accuracy (AUC = 0.7723) (Supplemental Fig. S3A). The logistic regression coefficients indicated that the most significant predictors were the mean signal strength of the four histone marks, whereas their coefficients of variation did not carry significant predictive power (Fig. 2C,D). These results suggest that the epigenetic strength at the proximal promoter, and not its dynamic variation, contributes to genes’ transcriptional response to an extracellular signal.

    The dynamics of distal epigenetic landscape associated with VEGFA responsiveness

    To gain greater insight into the epigenetic regulation of the VEGFA transcriptional response, we next interrogated epigenetic marks on distal chromatin regions (±1–100 kb from TSS). Although distal regions have been illustrated to have a substantial function in regulating lineage-specific gene expression, their role in the regulation of signal-driven gene expression remains equivocal. Distal transcriptional regulatory regions were categorized into three functional groups based on their chromatin features: active enhancers (AEs; all H3K27ac peaks, with or without H3K4me1 co-occupancy), poised enhancers (PEs; H3K4me1 peaks without H3K27ac co-occupancy), and repressive elements (REs; H3K27me3 peaks) (Fig. 3A). In total, 57,189 AEs, 28,152 PEs, and 48,072 REs were defined. Each main cluster was further stratified based on the temporal pattern of their respective histone modifications (see Methods; Supplemental Table S3). For example, H3K27ac peaks present at all four time points were denoted AE: H0-1-4-12, whereas H3K27ac peaks present only at baseline (0 h) were named AE: H0. We found 51.98%, 85.88%, and 95.74% of AEs, PEs, and REs, respectively, changed their chromatin occupancy during the 12-h VEGFA treatment time course (Fig. 3A). To validate these findings, we used ChIP-qPCR to measure H3K4me1, H3K27ac, and H3K27me3 occupancy of randomly selected dynamic sites. The ChIP-qPCR results were largely consistent with the ChIP-seq results (Supplemental Fig. S3B). The high percentage of regions that changed over time illustrates the broad dynamics of epigenetic modifications present on distal regions after VEGFA treatment.

    Figure 3.
    View larger version:
    Figure 3.

    VEGFA altered the chromatin landscape at distal regulatory regions. (A) Heatmap of chromatin signatures at distal regulatory regions. Distal chromatin regions were classified into active enhancers (AEs), poised enhancers (PEs), and repressive elements (REs) based on H3K27ac, H3K4me1, and H3K27me3 occupancy (see text and Methods). Within each class of regulatory element, regions were grouped by their temporal pattern of occupancy by H3K27ac (AEs), H3K4me1 (PEs), or H3K27me3 (REs). These clusters were named by the time points in which the histone mark was present, for example, AE: H1-4 indicates an AE with H3K27ac present at 1 and 4 h. For each cluster, the average signals of H3K27ac, H3K4me1, H3K4me3, RNAP II, and H3K27me3 at each time point are displayed as a heatmap. The DEG heatmap shows the enrichment of DEGs within 10, 50, or 100 kb of the regions in the cluster. Enrichment P-values were calculated using Fisher's exact test. The super-enhancer (SE) heatmap shows the number of SEs that overlap the regions in each cluster. The GWAS heatmap shows the enrichment (odds ratio [OR]) of previously defined blood vessel–related GWAS loci in each distal regulatory element cluster. The number of regions within each cluster is displayed as a histogram at the right of the plot. (B) Percentage distribution of DEG targeted by dynamic and stable AE or PE enhancers. (C) Enrichment of members of DEG clusters among genes linked to clusters in which AEs were found at a single time point (AE clusters H0, H1, H4, and H12). Enrichment P-value was calculated using Fisher's exact test. (D) Impact of dynamic REs on gene transcription. For genes associated with indicated set of REs, the ratios of expressions were shown. The plot showed that the reduction of dynamic REs (H0, H1, H1-12, H0-1-12, H12) up-regulated their associated target genes, compared to stable REs present at all time points (H0-1-4-12). Fold change represented the ratio of RE-targeted gene FPKM at the time of RE disappearance relative to that ahead of or after RE disappearance. Wilcoxon rank-sum test was used to calculate the significance. (**) P-value <0.001.

    DEGs were associated with more enhancers than non-DEGs (Supplemental Fig. S3C). We further analyzed the relationship between the histone landscape at distal regulatory regions and VEGFA-induced gene expression changes. The enrichment of DEGs within 10, 50, or 100 kb of distal regulatory regions within each temporal cluster was calculated using Fisher's exact test (Fig. 3A, columns labeled “DEG”). A subset of regulatory element clusters was significantly enriched with DEGs; among these were dynamic AEs and PEs, especially clusters with peaks at a single time point (outlined by black rectangles in Fig. 3A). Indeed, the large majority of DEGs were potentially regulated by AEs (88%) or PEs (81%) that belonged to “dynamic” clusters (Fig. 3B), that is, clusters other than H0-1-4-12 that were occupied by AE or PE epigenetic marks at some but not all time points. Moreover, the temporal pattern of H3K27ac on AEs was concordant to the temporal pattern of DEG expression (Fig. 3C), whereas this relationship was not observed for PEs or REs (Supplemental Fig. S3D). Though REs were not significantly associated with DEGs, VEGFA did stimulate dramatic H3K27me3 alterations across the genome, especially at distal and DNA repeat regions (Fig. 3A; Supplemental Fig. S3E). The expression changes of genes associated with dynamic REs were significantly higher than those with static REs, evidence of the functional impact of these REs on target gene regulation (Fig. 3D). Moreover, the dynamic AEs and PEs along with static regions were also associated with sequence polymorphisms associated with diseases relevant to blood vessels in genome-wide association studies (GWAS), suggesting that these regions are functionally relevant to these diseases (see Supplemental Table S3).

    Taken together, these results indicated that VEGFA induced widespread alterations of the epigenetic landscape at distal regulatory regions, which was associated with dynamic gene transcription and critical for regulating VEGFA responses.

    Super-enhancer dynamics mediates angiogenic responses

    Dense collections of enhancers, termed super-enhancers (SEs), activate the expression of master regulators (MRs) of cell lineage specification and function in development and disease (Whyte et al. 2013). Little is known about the identity and function of SEs in angiogenesis and VEGFA signaling. Based on H3K27ac ChIP-seq signal, we identified 1738 SEs over the four time points (Fig. 4A; Supplemental Fig. S4A; Supplemental Table S3). In line with the canonical role of SEs in other systems, endothelial cell SE-associated genes were highly expressed and function in pathways related to angiogenesis and blood vessel development, and many SE-associated genes, such as VWF, NOTCH1, CDH5, ERG, and DUSP1, encode known key endothelial cell regulators (Fig. 4A,B; Supplemental Fig. S4A–C). Forty-eight percent of SEs were present at all four time points, while the remaining 52% were dynamic under VEGFA stimulation, changing at least at one time point (Fig. 4C), and the average H3K27ac signal at all SEs (Fig. 4D and SE targeted genes; Supplemental Fig. S4D) changed as well across the time points. These dynamic SEs were associated with dynamic AEs (Fig. 3A) and notably were associated with dynamic PEs, particularly those clusters in which PEs were found at a single time point (Fig. 3A). All these together suggested that VEGFA-induced global changes in SEs are coordinated.

    Figure 4.
    View larger version:
    Figure 4.

    SEs preferentially regulate DEG. (A) HUVEC SEs at baseline. AEs were ranked by H3K27ac ChIP-seq signal, and nearby enhancers were stitched together using the ROSE algorithm (see Methods). Representative SE genes with critical vascular functions are highlighted. (B) Gene Ontology terms enriched among SE genes at baseline. Angiogenesis-related terms were significantly overrepresented. (C) Venn diagram showing the overlap between SEs at the four time points following VEGFA treatment. (D) Average plot of H3K27ac on SE regions at hours 0 (H0), 1, 4, and 12. VEGFA changed H3K27ac occupancy of SEs. (E) DEGs were enriched among SE genes. P-value was calculated using the χ2 test. (F) Chromatin landscape of the DUSP1 SE during VEGFA stimulation. S1–S9 represent sites of component enhancers of the SE (E1–E9). (G) 3C experiments showing the connection between SE component enhancers (E1–E9) and DUSP1 promoter. E1, E2, and E4–E9 interacted with the promoter, and E2 and E4–E8 increased interaction frequency at 1 h after VEGFA stimulation, when DUSP1 expression peaked. n = 3. Student's t-test: (*) P-value <0.05 compared to H0. (H) Luciferase assay testing the transcriptional activity of DUSP1 component enhancers. n = 3. Student's t-test: (*) P-value <0.05 compared to H0.

    A disproportionate fraction of DEGs was associated with SEs (22%; χ2 test P-value < 0.0001) (Fig. 4E). To gain additional insights into the participation of SEs in signal-induced gene regulation, we studied the regulation of DUSP1, a SE-associated gene that is rapidly up-regulated at 1 h after VEGFA treatment and then rapidly down-regulated (expression cluster G5) (Figs. 1A, 4F). The DUSP1 SE has nine component enhancers (E1–E9) (Fig. 4F). Although the DUSP1 SE was present at all time points, H3K27ac occupancy of some component enhancers (E2, E5, E7, E8) changed between time points (Fig. 4F). We used chromatin conformation capture to measure the contact frequency between the component enhancers and the DUSP1 promoter (E7 and E8 were not separately resolved in this assay and are collectively referred to as E7/8) (Fig. 4G). Eight of the nine component enhancers (all but E3) looped to the DUSP1 promoter, suggesting that these component enhancers contributed to regulating DUSP1 expression. Among these eight, four (E2, E5, E6, and E7/8) increased their promoter contact frequency after VEGFA treatment, and the chromatin conformation and gene expression kinetics matched, with each peaking at 1 h. Moreover, the component enhancers that displayed dynamic H3K27ac occupancy all exhibited dynamic promoter contacts. To further analyze the contribution of these component enhancers to DUSP1 transcriptional regulation following VEGFA treatment, we measured their transcriptional activity using the luciferase assay. We cloned the nine component enhancers upstream of the DUSP1 promoter and luciferase. Six component enhancers (E2, E4, E5, E6, E7, E8) had transcriptional activity that increased after VEGFA treatment, and the kinetics of five of these (E2, E5, E6, E7, E8) matched the temporal gene expression signature (Fig. 4H). These data suggest that E2, E5, E6, E7, and E8 all contribute to VEGFA activation of DUSP1 by gaining H3K27ac, increasing promoter contact frequency, and driving DUSP1 promoter activity. While the DUSP1 SE creates a primed environment for rapid DUSP1 up-regulation, dynamic changes at component enhancers rapidly create a functional enhancer that drives productive gene transcription. Together, these results demonstrate that dynamics of component enhancers within SEs is an important and novel epigenetic feature that regulates transcription in response to extracellular cues such as VEGFA.

    Dynamic TF chromatin occupancy coordinates epigenetic and transcriptional changes

    TFs recognize specific sequence motifs to direct the shaping of the chromatin landscape and transcriptome. To further elucidate the TF network in endothelial cells in regulation of dynamic enhancers and transcription responsive to VEGFA, we set out to map the genome-wide chromatin occupancy of key endothelial TFs. We first scanned enhancer regions for overrepresented TF motifs. Out of 2470 motifs compiled from multiple sources (see Methods) that included TF heterodimer motifs (Jolma et al. 2015), we found 74 nonredundant motifs that were highly overrepresented (P-value < 10−15 in regions from at least one ChIP-seq sample). Based on an analysis of the co-occurrence of these motifs within the same regions, we were able to cluster them into six main groups (Fig. 5A). Many of these motifs contained the recognition sites of ETS factors, which highlights the dominant role of ETS family TFs in endothelial cells (De Val et al. 2008; De Val and Black 2009; Zhou et al. 2017). Most of these ETS-related motifs were heterodimers in which one partner was an ETS family member and the other belonged to other TF families. In addition to the recognized FOX-ETS motif, additional heterodimer motifs suggested interaction between ETS factors and members of the TEAD, HOX, SOX, BHLH, and TCF families. Moreover, we observed different clusters of co-occurring ETS heterodimer motifs, which points to functional differences between them.

    Figure 5.
    View larger version:
    Figure 5.

    TF dynamics are linked to VEGFA-induced changes in the chromatin landscape and transcriptome. (A) Co-association of motifs at enhancers. AEs were scored for presence of motifs, and the frequency that regions contained two different motifs was calculated. Six major groups of co-associated TFs were found (green rectangles). The members of these groups are labeled on the right. (B) Chromatin occupancy of six TFs and two epigenetic enzymes. Regions were grouped based on their temporal pattern of factor binding, from binding at all time points (“stable,” top) to binding at some but not all time points (“partially_shared”), to binding at a single time point (“unique”). VEGFA stimulated dramatic changes in chromatin occupancy by these eight proteins, as demonstrated by the large majority of regions being bound at a subset of time points (“partially_shared” and “unique”). (C) Linkages between ETS1, JUN, EP300, and EZH2 peaks and AEs visualized in Circos plots (Krzywinski et al. 2009). TFS and AE clusters were defined as in Figures 4A and 6A. Each line represents an overlap between a TF and an AE cluster and is colored by the associated TF cluster. The bar plots on the inner side of the circle indicate the enrichment value of TF clusters on AE or RE clusters. (D) The enrichment of tested transcriptional regulators at SEs and typical enhancers (TEs). The size of the dot represents the TF enrichment on both TEs and SEs. The enrichment score was calculated as the sum of peaks on TEs and SEs normalized by the total number of peaks across the whole genome. The dots above and below the diagonal line indicate preferential enrichment on TEs or SEs, respectively. (E) Enrichment of DEGs among genes’ neighboring regions bound by the indicated transcriptional regulator and bound in the indicated temporal pattern. DEGs were mapped to bound regions within 2, 10, 50, or 100 kb. Temporal clusters are named by time points with regulator binding, for example, H0-1-4 indicates that a peak was called on member regions at H0, H1, and H4. Enrichment P-value was calculated by Fisher's exact test. Multiple temporally dynamic TF clusters, particularly the time unique clusters, were overrepresented for DEGs. (F) Enrichment of DEG clusters among genes’ neighboring regions bound by the indicated transcriptional regulator and bound in the indicated temporal pattern. This panel is similar to D, except that DEGs are separated by the temporal expression pattern. Rectangles highlight enrichment of DEG clusters that are referred to in the text.

    To test the hypothesis that altered TF occupancy contributed to enhancers dynamics, we selected six overrepresented TFs in Figure 5A (ETS1, ERG, FLI, GATA2, JUN, and RBPJ) with known functions in angiogenesis and profiled their genomic occupancy in HUVEC by ChIP-seq at 0, 1, 4, and 12 h after VEGFA stimulation. At the same time, we also performed ChIP-seq on two epigenetic modulators, EZH2 and EP300 (Fig. 5B; Supplemental Fig. S5A). Specificity of ChIP antibodies was validated by western blot (Supplemental Fig. S5A). As a result, the protein levels of several of these factors (EZH2, GATA2, and RBPJ) were altered by VEGFA (Supplemental Fig. S5A). Overall, we obtained 75,739 TF binding peaks (Supplemental Fig. S5B; Supplemental Table S3). Motif scanning of these peaks illustrated a high enrichment of their cognate DNA-binding motifs (Supplemental Fig. S5C), reflecting the high quality of these ChIP-seq data.

    To reveal how VEGFA stimulation influences the binding pattern of these TFs over time, we grouped TF-bound regions based on TF occupancy at each time point (Fig. 5B). The majority (83.91%–97.20%) of TF-bound regions changed TF occupancy at least at one time point. These observations were confirmed by ChIP-qPCR, which showed that 70%–90% of unique sites changed TF occupancy, whereas stable sites did not (Supplemental Fig. S5D). The dynamic binding at these sites was further manifested by low correlations between time points in unique clusters versus high correlations in stable clusters (Supplemental Fig. S5E). These data indicate that VEGFA altered global TF activity, which is similar to epigenetic modifications. To analyze the relationship of TF binding to the epigenetic landscape, we compared TF binding clusters to AE clusters. Several TFs, most notably JUN, FLI, and ETS1, as well as the histone acetyltransferase EP300, were enriched at both dynamic and stable AEs (Fig. 5C; Supplemental Fig. S6A,B). Dynamic TF activity was linked to significant changes of H3K27ac and H3K4me1 (Supplemental Fig. S6C). Except for EZH2, these factors were also dynamically enriched at both typical AEs and SEs (Fig. 5D). Some factors, such as GATA2, had greater enrichment in typical enhancers (TEs) compared to SEs, whereas others, such as ERG1 and EP300, had greater enrichment in SEs compared to TEs. Unlike the other factors profiled, the repressive factor EZH2 was not enriched in either TEs or SEs (Fig. 5D). Rather, EZH2 was enriched at static REs, as 86% of EZH2 bound regions were contained within static REs (Fig. 5C; Supplemental Fig. S6B). These results together suggested that dynamics of TF and epigenetic enzyme chromatin occupancy contribute to VEGFA-induced changes in the chromatin landscape, particularly at AE and SE sites.

    To test the possibility that changes in TF occupancy and epigenetic state contribute to the VEGFA transcriptional response, we interrogated the relationship between TF binding clusters and DEGs (Fig. 5E). We found that DEGs were highly enriched adjacent to regions with dynamic binding by sequence-specific TFs and EP300 (but not EZH2). The enrichment was strongest in clusters with TF binding at a single time point (H0, H1, H4, and H12) (Fig. 5E). We looked deeper into this relationship by determining the enrichment of DEG clusters G1–G7 among genes adjacent to temporal clusters of each TF (Fig. 5F). We found enrichment of specific, dynamic TF binding patterns with specific DEG clusters. For example, in line with our previous study that ETS1 positively regulates global transcription, ETS1 mainly correlated with up-regulated genes (Chen et al. 2017). Expression cluster G4, which featured transient up-regulation of genes at 1 h, was associated with transient occupancy of JUN, FLI, or EP300 at 1 h (Fig. 5F, white boxes; Supplemental Fig. S6D). The G4 cluster was also enriched among regions that bound GATA2 at hour 0 and lost GATA2 at subsequent time points (Fig. 5F, yellow box), consistent with a repressive role of GATA2 (Obara et al. 2008). On the other hand, regions that lost GATA2 binding were also overrepresented for expression clusters G2 and G3, which were down-regulated after VEGFA (Fig. 5F, dotted yellow box; Supplemental Fig. S6D), indicating GATA2-dependent gene activation at these sites. RBPJ sites present exclusively at either hour 0 or hour 12 were also overrepresented for expression cluster G4 (Fig. 5F, blue box; Supplemental Fig. S6D), consistent with a model in which RBPJ represses transcription in the absence of Notch activation (Kao et al. 1998). Hence, this study pinpointed divergent roles of the tested TFs in regulating gene expression downstream from VEGFA signaling.

    MAFs are master regulators of VEGFA transcriptional responses

    Transcriptional regulatory networks play essential roles in many physical and pathological processes. Within these networks, master regulators (MRs) play dominant roles in establishing the epigenetic landscape and regulating the transcription of key functional genes (Lefebvre et al. 2010). To identify MRs that control VEGFA-induced responses in endothelial cells, we integrated DEGs, dynamic TF-AE linkages, motif-dynamic AE linkages, AE-DEG associations, and protein–protein interactions from the STRING database (https://string-db.org/) to construct a transcriptional network (see Methods) (Fig. 6A). This network contained six major modules with functions relevant to cell proliferation, migration, and oxidative stress that have been proven to be essential processes of angiogenesis (Supplemental Fig. S7A).

    Figure 6.
    View larger version:
    Figure 6.

    MAF factors are master regulators (MRs) of endothelial cell VEGFA responses and angiogenesis. (A) The VEGFA transcriptional gene regulatory network. Nodes are DEGs, expressed TFs whose motifs are significantly enriched at AEs, or proteins that interact with the aforementioned factors. MR nodes identified by MARINa are highlighted in magenta, whereas other nodes are colored green. Edges represent TF–target gene interactions, as determined by ChIP-seq or inferred from enriched AE motifs. Edges with a positive expression correlation are colored red, and those without are colored green. Light blue edges represent protein–protein interactions from the STRING database. Node size represents the number of connected edges. Transparent blue ellipses circled different modules, and MRs in each module are highlighted nearby in black. (B) MRs of VEGFA transcriptional responses in endothelial cells inferred by MARINa. The blue heatmap (left) represents the relative expression of the indicated TFs over the time course. The orange heatmap (right) displays the relative enrichment score (M-score) for TF targets among DEGs at hours 1, 4, and 12 compared with hour 0. Positive and negative M-scores indicate positive or negative correlation between a TF and its target DEGs at the labeled time point. All MAF family members (MAFF, MAFG, and MAFK) were identified as MRs. (C) Overlap between DEGs and predicted MAF targets, inferred by the presence of the MAF motif in adjoining AEs. Enrichment P-value over genome-wide expectation was calculated using a Fisher's exact test. (D) Enriched Gene Ontology terms of predicted MAF target genes (adjacent to AEs containing the MAF motif). The top 15 most significant GO terms are shown. (E) Inhibition of MAFF or MAFG, but not MAFK, suppressed HUVEC cell migration induced by VEGFA. Cell migration was measured using the transwell assay. (F) Overexpression of MAFF, MAFG, or MAFK augmented VEGFA-induced migration, as measured by the transwell assay. (G) Inhibition of MAFF or MAFG, but not MAFK, dampened HUVEC cell proliferation. Cell number was measured using a tetrazolium dye reduction assay. (H) Lentiviral overexpression of MAFF, MAFG, or MAFK promoted HUVEC cell proliferation. (I) Matrigel plug assay to assess role of MAF proteins in angiogenesis in vivo. (I–M) Matrigel plug assay to assess role of MAF proteins in angiogenesis in vivo. HUVECs with knockdown or overexpression of MAFF and MAFG were mixed with stromal cells and Matrigel and injected subcutaneously into nude mice. One week later, vascular density (J), the fraction of MKI67 (also known as Ki-67) positive human endothelial cells (K), and caspase 3-positive human endothelial cells (L) were measured. PECAM1 staining was used to evaluate the invaded blood vessels from host mice (M). n = 5. Student's t-test: (**) P-value <0.01; (***) P-value <0.001. (I) Representative images of immunostained Matrigel plug sections.

    We further interrogated the relationship between expressed TFs and DEGs using MARINa, a MR analysis algorithm (Lefebvre et al. 2010). This analysis identified 41 MRs (Fig. 6B) with predicted enriched regulation of DEGs, measured as the M-score. Among these are TFs that are known to be central regulators of endothelial cells, including ETS1, MEF2C, FOXO1, JUN, ERG, and MAX (De Val et al. 2004; Birdsey et al. 2012; Wilhelm et al. 2016; Chen et al. 2017). Other predicted MRs do not have well-described roles in ECs or angiogenesis, such as MAFF, MAFK, MAFG, NR5A2, GABPA, and TGIF2. Some MRs were dynamically expressed in ECs after VEGFA treatment, and their M-score varied over time, suggesting temporally regulated MR activity (Fig. 6B). Nearly all of the endothelial MRs from MARINa formed highly interconnected nodes in the network, which independently validated the MARINa predications and strengthened the central role of these MRs in orchestrating VEGFA responses in endothelial cells.

    Both MARINa and the network analysis identified MAFK, MAFF, and MAFG, which compose the small MAF family, as highly interconnected MRs. These basic leucine zipper domain TFs bind to the Maf recognition element (MARE) and are known to regulate antioxidant responses, metabolism, aging, and tumorigenesis (Katsuoka and Yamamoto 2016). However, their roles in endothelial cell biology and angiogenesis have not been described previously. MAFF and MAFK were also targeted by SEs (Fig. 4A; Supplemental Fig. S4A). Moreover, all three MAFs were expressed in aortic, venous, and microvascular ECs and up-regulated by VEGFA (Supplemental Fig. S7B–D). The genes neighboring MAF motif–positive AEs were preferentially DEGs (χ2 test P-value < 2.2 × 10−16) (Fig. 6C) and functionally related to angiogenesis (Fig. 6D; Supplemental Table S2). Depleting MAFs using siRNAs altered the VEGFA response of six tested target genes (Supplemental Fig. S8A,B). Together these data indicate that the small MAFs are indeed MRs central to VEGFA signaling and may have essential functions in angiogenesis.

    To further test the function of small MAFs in angiogenesis, we measured the effect of their gain and loss of function in multiple angiogenesis assays. Individual knockdown of MAFF and MAFG, but not MAFK, by siRNA disrupted VEGFA-induced migration of HUVECs (Fig. 6E; Supplemental Fig. S9A). When all three MAFs were depleted, VEGFA-stimulated cell migration was almost completely abolished. On the other hand, up-regulation of each MAF stimulated cell migration (Fig. 6F; Supplemental Fig. S9B–D). These data indicate that the small MAFs have partial functional redundancy and collectively are critical for VEGFA-driven cell migration. Proliferation assays showed that MAFF and MAFG, but not MAFK, are similarly essential for VEGFA-stimulated endothelial cell proliferation and that their overexpression is sufficient to enhance this process (Fig. 6G,H). MAFs were reported to protect various type of cells from oxidative stress by up-regulating the expression of antioxidant genes (Katsuoka and Yamamoto 2016). In line with this effect, we also found that overexpression of small MAFs protected HUVECs from H2O2 damage (Supplemental Fig. S9E). To further probe if the proangiogenic effect of MAFs in VEGFA signaling required their transcriptional activity, we transfected endothelial cells with synthetic DNA decoys containing the MARE motif to determine if VEGFA activity requires MAF binding to DNA (Supplemental Fig. S9F). We found that MARE decoy DNA markedly reduced of endothelial cell proliferation compared to control DNA, indicating that the angiogenic activity of MAFs depended on their chromatin occupancy and subsequent transcriptional activity (Supplemental Fig. S9G). The proangiogenic activity of MAFs could be also mediated by enhanced expression of angiocrine factors secreted by endothelial cells, which recently has been revealed to be important in neovascularization (Rafii et al. 2016). To test this possibility, we profiled commonly studied angiogenic factors in MAF-overexpressing ECs. We found that three proangiogenic angiocrine factors (VEGFA, IL6, and ANGPT1) were down-regulated following the up-regulation of MAFs, whereas three factors (TGFB1, BMP2, and JAG1) with negative effects on angiogenesis were down-regulated (Supplemental Fig. S9H). These findings suggest that MAFs are unlikely to promote angiogenesis by stimulating angiocrine factor expression but rather by stimulating EC migration and proliferation directly.

    We used an in vitro tube formation assay to further examine if MAFs promote the formation of vascular networks. HUVEC tube formation was significantly attenuated by knocking down MAFF or all three MAFs simultaneously and was enhanced by overexpression of MAFF or MAFG (Supplemental Fig. S10A,B). To validate the proangiogenic roles of MAFs in vivo, we used the Matrigel plug assay. Consistent with in vitro experiments, depletion of MAFF, MAFG, or both in endothelial cells reduced cell proliferation, increased cell death, and dampened vascular tube formation within the transplanted Matrigel (Fig. 6J–L). In contrast, overexpression of either MAFF or MAFG dramatically enhanced new vessel formation by promoting cell growth and reducing apoptosis. Moreover, HUVECs overexpressing MAFs also stimulated host endothelial cell recruitment into the Matrigel plug, and MAF knockdown in HUVECs reduced host endothelial cell recruitment (Fig. 6I,M), further demonstrating the proangiogenic role of MAFs.

    Together, the data in this section demonstrate that the small MAFs are novel MRs of VEGFA signaling that augment angiogenesis.

    Previous SectionNext Section

    Discussion

    Here, we used multiple genome-wide approaches to illuminate the global alterations in epigenetic and transcriptional landscapes induced by VEGFA stimulation. Our study highlights the interplay between TFs and histone modifications at distal regions that underlies rapid transcriptional responses to extracellular cues. Furthermore, integrating these multiple data types yielded a core VEGFA-driven transcriptional network that identified a handful of master endothelial cell transcriptional regulators. Among these were the small MAF proteins, which we found to be critical for angiogenesis. Together, our results shed new mechanistic insights into transcriptional mechanisms responsible for VEGFA-stimulated angiogenesis.

    A dynamic epigenetic and TF program controls angiogenesis

    The impact of transient extracellular stimuli on the chromatin landscape and the functional importance of these epigenetic changes are critical understudied areas of gene regulation. We and others have shown that extracellular stimulation induces dynamic changes of H3K27ac that mark functional transcriptional enhancers (Ostuni et al. 2013; Zhang et al. 2013). However, whether other epigenetic modifications are also altered, where these alterations occur, and how they contribute to signal-induced transcriptional responses are key remaining questions. By comprehensively measuring endothelial cell responses to VEGFA, we addressed these questions by defining acute signal-induced changes in chromatin occupancy by TFs, epigenetic regulators, and histone marks and by delineating how these changes are related to signal-induced alterations in the transcriptome.

    Our study shows that although the epigenetic landscape provides a permissive environment for transcription, the dynamics of these epigenetic marks are not evidently associated with the temporal expression. Rather, we found that changes in TF binding and epigenetic regulation at distal regions were significantly linked to VEGFA responses with ample evidence: More enhancers were associated with DEGs, dynamic AEs and PEs overlapped with majority of DEGs, the temporal pattern of dynamic AEs is concordant with DEGs, and the dynamic REs destabilized their neighbor genes than did consistent REs.

    The relationship between distal regions and signal-responsive gene transcription was further supported by our SE study, which showed that SEs were substantially changed by VEGFA and that SEs were highly enriched for VEGFA-responsive genes. Our study of the DUSP1 SE, which was present at all time points, illustrates that this group of SEs also underwent significant signal-induced changes in component enhancer–promoter contacts, which correlated with changes in enhancer activity and gene transcription. These data suggest that SEs provide a platform for rapid, VEGFA-induced transcriptional changes and that stimulation of productive gene transcription is governed by local changes in chromatin conformation that regulate enhancer–promoter contact.

    Recent studies have shown that functional distal regions are key for the intricate regulation of transcription in lineage specification and disease states (Visel et al. 2009; Hnisz et al. 2013). Our study suggests that distal regions are likewise the major drivers of signal-responsive gene regulation.

    Small MAF factors: new regulators of angiogenesis

    MAFG, MAFK, and MAFF, composing the small MAF family of bZIP TFs, have been shown to play critical roles in development (Katsuoka and Yamamoto 2016). Mice lacking these three MAFs died at E13.5 with multiple organ deficiencies, but the effect on developmental angiogenesis was not evaluated (Onodera et al. 2000). Our integrative transcriptional network analysis and subsequent validation experiments demonstrate that small MAFs have critical functions in angiogenesis by promoting endothelial cell migration, proliferation, and tube formation. This activity appears to be shared with other extended MAF family members, because MAFB, a large MAF factor, had been also recently found to activate vascular sprouting in rodents (Jeong et al. 2017).

    MAFs bind DNA but do not themselves activate transcription; instead, transcriptional activation functions are conferred by interacting proteins, such as NRF2 or BACH1/2 (Igarashi et al. 1994; Sun et al. 2002). NRF2 is a key factor that regulates expression of antioxidant genes (DeNicola et al. 2011), and we showed that MAF factor overexpression increased HUVEC cell resistance to oxidant stress and prevented apoptosis within Matrigel plugs. BACH1 has been recently reported to promote angiogenesis (Jiang et al. 2015). Determining whether these factors, or other MAF-interacting factors, are essential for the angiogenic activity of MAFs will be a valuable direction for future studies.

    Excessive angiogenesis is critical for the pathogenesis of many forms of cancer and retinopathy, whereas other diseases such as stroke and ischemic cardiomyopathy result from inadequate organ perfusion. As a result, therapeutic strategies are needed to increase or inhibit angiogenesis. Thus, our finding that the small MAF proteins are MRs of endothelial cell responses to VEGFA has the potential to lead to new approaches to therapeutically modulate angiogenesis.

    In summary, our study identified temporal transcriptional and epigenetic changes at distal regions that control the response of cells to VEGFA stimulation, laying the foundation to decipher the chromatin code that controls developmental and pathological angiogenesis.

    Previous SectionNext Section

    Methods

    Cell culture

    Primary HUVECs (CC-2519), HPAECs (CC-2530), and HMVECs (CC-2930) were purchased from Lonza and grown in EGM2 full medium (CC-3162, Lonza). For VEGFA stimulation, HUVEC cells at three to six passages were grown overnight in low serum (EBM2 with 1% FBS) and then treated with 50 ng/mL VEGFA. HUVEC cells were collected at 0 (untreated), 1, 4, and 12 h after addition of VEGFA for downstream analyses. Where indicated, flavopiridol (100 nM) and JQ1 (500 nM) were added 1 h before VEGFA stimulation.

    The siRNAs listed in Supplemental Table S4 were used at 10 ng/mL and transfected with RNAiMAX (Thermo Fisher Scientific) following the manufacturer's instructions. Lentivirus was prepared by PEG6000 precipitation and used to treat HUVECs at 2–5 MOI, in the presence of 8 µg/mL hexadimethrine bromide (107689, Sigma-Aldrich). Migration and proliferation experiments were performed 24 h after virus treatment.

    Mice

    All experiments with mice were performed under protocols approved by the institutional animal care and use committee of Shanghai Jiao Tong University. Matrigel plug assay was performed as described previously (Chen et al. 2017) using 6- to 8-wk-old male nude mice (Si Lai Ke Experimental LLC). In brief, 1 × 106 HUVEC cells at passages 4–6 and 2 × 106 human mesenchymal stem cells (Lonza) were mixed with 200 µL ice-cold Matrigel (Corning 356237) and injected subcutaneously. After a week, the Matrigel plugs were dissected out, fixed with 4% paraformaldehyde overnight, and paraffin embedded. HUVECs were stained using Ulex europaeus agglutinin I (1:200, Vector Labs), a lectin that recognizes human but not mouse endothelial cells. Blood vessels from host mice were stained with PECAM1 (1:200, BD Biosciences). Proliferating cells were stained using Ki-67 antibody (1:200, Abcam), and apoptotic cells were stained using antiactivated caspase 3 antibody (1:200, Millipore). Images were acquired with a Nikon A1Si confocal microscope. The percentage of proliferating or apoptotic endothelial cells and the total vascular area were measured using ImageJ.

    ChIP-seq

    ChIP-seqs were performed as previously described with minor modifications (Zhang et al. 2013). In brief, HUVECs (5 × 107 cells for TF ChIP, 1 × 107 cells for histone ChIP) were crosslinked with 1% formaldehyde for 10 min at RT and neutralized with 0.125 M glycine. The nuclei were purified with 5 mL hypotonic buffer (20 mM HEPES at pH 7.5, 10 mM KCl, 1 mM EDTA, 0.2% NP-40, 10% Glycerol, 1 × Protease Inhibitor Cocktail [PI; Roche]) and then sonicated with a Misonix Sonicator 3000 (pulse on: 7 min; amplitude: 70) in sonication buffer (20 mM Tris Cl at pH 8.0, 2 mM EDTA, 150 mM NaCl, 1% NP-40, 0.1% SDS, 1 × PI). Sheared DNA was precleared with 5 µL of Protein G Dynabeads (Thermo Fisher Scientific) for 1 h, incubated with 5–10 µg ChIP antibody (Supplemental Table S5) overnight at 4°C, and then rinsed three to five times with RIPA buffer (50 mM HEPES at pH 8.0, 1 mM EDTA, 1% NP-40, 0.7% sodium deoxycholate, 1% Triton X-100) containing 1 × PI. Bead-bound DNA was decrosslinked in 50 mM Tris-Cl (pH 8.0), 1% SDS for 12–16 h at 65°C and then purified with QIAquick PCR purification columns. For ChIP-qPCR, the purified ChIP DNA was quantified by real-time PCR using primers listed in Supplemental Table S4.

    For ChIP-seq, ChIP DNA was converted into Illumina sequencing libraries using the NEBNext ChIP-seq library prep master mix set for Illumina. The barcoding primers were custom-synthesized (Supplemental Table S4). The libraries were sequenced on Illumina HiSeq 2000 and 2500 instruments to yield 50-nt single-end reads. Two biological replicates of EP300 were sequenced to verify reproducibility (Zhang et al. 2013). One biological replicate was used for the other chromatin features (Supplemental Table S1).

    RNA expression profiling

    Total RNA was recovered from HUVECs using the RNeasy mini kit (Qiagen) with on-column DNase I digestion. Polyadenylated RNA was purified from 5 µg total RNA using the Dynabeads mRNA DIRECT purification kit (Thermo Fisher Scientific). The oligo(dT)-based purification step was performed twice to minimize rRNA contamination. RNA-seq libraries were synthesized using the ScriptSeq V2 RNA-Seq library preparation kit (Illumina) and custom-synthesized barcoding primers (Supplemental Table S4). The libraries were sequenced on Illumina HiSeq 2500 instruments to yield 50-nt paired-end reads. Two biological replicates were sequenced at each time point.

    Transwell migration assay

    HUVECs were transduced with lentivirus or transfected with siRNA for 48 h, trypsinized, and resuspended into migration medium (EBM2 + 1% FBS). We placed 1 × 105 HUVECs in the upper side of transwell chambers (Corning), and 30 ng/mL VEGFA was added to lower wells to induce transmembrane migration. After 12 h, migrated cells were stained with crystal violet (0.09% crystal violet in 10% ethanol) for 1 h and rinsed with water. Five randomly selected regions were imaged with a Nikon SMZ800N microscope and quantified using the the “particle analysis” function of ImageJ.

    In vitro proliferation assay

    Three thousand HUVECs were added to each well of a 96-well plate 24 h after lentivirus transduction. At the indicated time point, the cell number was quantified using a tetrazolium dye reduction assay (Beyotime C0042). OD values were converted to cell number using a standard curve drawn with serially increasing cell numbers.

    Tube formation assay

    The tube formation assay was performed as previously described (Zhang et al. 2009). Briefly, 100 µL of ice-cold Matrigel (Corning 356237) was mixed with 2 volumes of EGM2 and then seeded in 24-well plates. After gelatinization, 2 × 105 MAF-transfected HUVEC cells were plated onto the surface of the gel for the formation of tubes. The tubes were imaged under a Nikon SMZ800N microscope 6–16 h after seeding, and tube length was measured with ImageJ software.

    Protein quantification

    Total protein was extracted from HUVECs treated with VEGFA (50 ng/mL) using RIPA buffer containing 1 mM PMSF and 10 mM NaF. Total protein concentration was determined with the BCA protein assay (Pierce 23225). Equal amounts of protein were loaded into each well of 4%–12% gradient SDS-PAGE gels and then transferred to PVDF membranes (Millipore). The membranes were probed with specific primary antibodies (Supplemental Table S5) and appropriate HRP-conjugated secondary antibodies and then developed using Immobilon Western Chemiluminescent HRP substrate (Millipore).

    Chromatin confirmation capture (3C)

    Chromatin conformation capture (3C) was performed as previously described (Zhang et al. 2013). We crosslinked 1 × 107 HUVECs with 2% formaldehyde for 10 min at room temperature. After inactivating formaldehyde with 0.125 M glycine, nuclei were extracted with nuclei extraction buffer (10 mM Tris-Cl at pH 8.0, 10 mM NaCl, 5 mM MgCl2, 0.2% NP-40, 0.1 mM EGTA), resuspended in NEB restriction enzyme buffer, and digested with 1500U HindIII overnight at 37°C with constant agitation. After inactivation of HindIII with 1.6% SDS at 65°C, the digested DNA was ligated with 100 Weiss units of T4 DNA ligase for 4 h at 16°C in 7 mL 1 × ligation buffer. The ligated DNA was decrosslinked by incubating overnight with 300 µg/mL Proteinase K at 65°C. DNA was then purified by pheno/chloroform extraction and ethanol precipitation.

    BAC clones containing genomic loci of DUSP1 (Thermo Fisher Scientific, clone 3091L13) and VWF (CHORI, clone RP11-1137J12) were digested with HindIII and religated with T4 DNA ligase to generate control templates for quantification of results across the restriction sites. The chromatin conformation contact frequency was measured by quantitative TaqMan PCR using primers and probes listed in Supplemental Table S4 and calculated as described by Zhang et al. (2013).

    Transcriptional activity assay

    The promoter of DUSP1 was amplified from a BAC clone (Thermo Fisher Scientific, clone 3091L13) and inserted into pGL4-Luc2P (Promega) to generate a pGL4-DuspPro vector. The component SE regions of DUSP1 were further inserted upstream of the DUSP1 promoter to produce pGL4-SE-DuspPro. Equal amounts of pGL4-SE-DuspPro and pGL4-CMV-hRluc were transfected into HUVECs with Lipofectamine LTX (Life Technologies 15338-100). Two days after transfection, HUVECs were serum starved overnight and then stimulated with 50 ng/mL VEGFA for 0, 1, or 4 h. Luciferase activity was measured using the Luciferase dual assay system (Promega E1910). Firefly luciferase activity was normalized to the activity of cotransfected Renilla luciferase.

    Computational analyses

    Many of our high-throughput tracks were generated and analyzed before 2013, when the GRCh38 released. After testing realignment of TF sequencing files to GRCh38, extremely few ChIP-seq peaks (0.08%–0.4%) called with GRCh38 differed from those called with hg19, and these rare differences did not alter the conclusions. Therefore, we used the hg19 as the genome reference in this study.

    ChIP-seq

    ChIP-seq reads were aligned to hg19 reference genome using Bowtie 2 (Langmead and Salzberg 2012) with default options. Duplicates were marked with Picard v1.129 (http://broadinstitute.github.io/picard). Low-quality (Q < 15) reads and duplicates were removed. MACS 1.4 (Zhang et al. 2008) was deployed to call ChIP peaks. We subtracted peaks that fell inside the ChIP-seq blacklist regions (https://sites.google.com/site/anshulkundaje/projects/blacklists).

    ChIP-seqs of histone markers (H3K27ac, H3K4me1, H3K4me3, and H3K27ac) were downloaded from ENCODE (The ENCODE Project Consortium 2012; www.encodeproject.org) under accession number ENCSR289LSP. The same protocol was applied on these data. Then, ChIP-seq correlations were calculated using deepTools2 (Ramirez et al. 2016).

    RNA-seq

    FASTQ files were aligned to hg19 reference genome using TopHat v2.1.0 (Trapnell et al. 2012) without discovery of novel splice junctions or novel indels. Gene abundances were quantified with Cuffquant with options “--multi-read-correct --frag-bias-correct.” Differential gene expression was assessed with Cuffdiff with options “--library-norm-method geometric --dispersion-method pooled --compatible-hits-norm --multi-read-correct --min-alignment-count 10 --FDR 0.1 --frag-bias-correct”. A total of 29,722 genes were expressed in our data set (FPKM ≥ 1 in at least one time point). Differential expression was assessed from pairwise comparisons to hour 0 with a fold change ≥2. DEGs were clustered with the partitioning around medoids method in the cluster R library (https://cran.r-project.org/web/packages/cluster/index. html) with Pearson correlation distance. The number of clusters was optimized by silhouette analysis.

    Differentially expressed lncRNAs were obtained from RNA-seq by using GENCODE gene annotations. We associated lncRNAs to neighboring genes within ±100 kb by applying the “Basal plus extension” mode of GREAT (McLean et al. 2010). The Pearson correlation coefficient (PCC) between expression of lncRNAs and neighboring genes expression was calculated. The PCC distributions of differentially expressed lncRNAs to their neighboring genes or randomly selected genes (neighboring genes excluded) were compared using the Kolmogorov–Smirnov test.

    Promoter chromatin feature analysis

    This analysis encompassed the 29,722 genes detectably expressed by RNA-seq. For genes with multiple potential transcripts, the ones with the highest H3K4me3 coverage within ±2 kb of annotated TSSs were used. For genes without significant H3K4me3 coverage, the TSSs of the longest transcripts were used. Coverage for H3K27ac, H3K4me1/2/3, H3K27me3, H3K36me3, and RNAP II was calculated in the 2-kb promoter region centered on the TSS by normalizing to sequencing depth and subtracting input signal. Promoters were grouped into nine clusters by the pattern of these chromatin features across the four time points using the partitioning around medoids method in the cluster R library with the PCC distances.

    To examine whether signal strength and change of epigenetic markers at promoters contributed to gene expression dynamics, we calculated the mean coverage and coefficient of variance of H3K4me2, H3K4me3, H3K27ac, and H3K27me3 at the promoters. A logistic regression binary classifier was applied to estimate the probability of a gene being differentially expressed based on the eight predicting variables (strength and variance of four epigenetic marks). By using 10-fold cross-validation, the performance of the model was measured by determining the area under the receiver operating characteristic (ROC) curve, that is, the plot of the true-positive rate (TPR) against the false-positive rate (FPR) at various threshold settings.

    Distal regulatory element analysis

    Distal regulatory elements (DREs) were defined by peak regions of H3K27ac, H3K4me1, and H3K27me3, after removing peaks whose centers were within ±2 kb of a TSS. The function MergePeaks (options: -d given) in HOMER (Heinz et al. 2010) was deployed to group the DREs by their patterns of peak overlap: AEs were defined as all H3K27ac peaks, PEs were all H3K4me1 peaks that did not overlap AEs, and REs as all H3K27me3 peaks that did not overlap AEs or PEs.

    To evaluate the correlation between DREs and DEGs, genes were assigned to enhancers by GREAT using distance thresholds of 10, 50, or 100 kb to their TSS. The enrichment of DEGs within genes assigned to each type of DRE was calculated using Fisher's exact test. To examine if dynamic DREs possibly regulate their target genes, we analyzed genes that were targeted by only one type of DRE. Change of gene expression associated with dynamic DREs was calculated as (FPKM before or after DRE presence)/(FPKM on DRE presence). Gene expression changes associated with dynamic DREs present at H0, H1, H12, H1+H12, and H0+H1+H12 were compared to those associated stable DREs present throughout the time course using the Wilcoxon rank-sum test.

    SE analysis

    SEs were called using the ROSE method and algorithm-based on H3K27ac signal and peaks (FDR = 5%) at each time point (Whyte et al. 2013). The overlap of SEs across the time points was calculated using BEDTools (intersected bed) with at least 1-kb overlap. The metagene plot of H3K27ac was generated with ngsplot (ngsplot.r –G hg19 –R bed) (Shen et al. 2014). To characterize the preference of TF binding to SEs or TEs, we calculated the enrichment score of each TF on either SEs or TEs: enrichment score = (|overlapped TF peaks|/AE peak length)/(|all TF peaks|/genome length). The peak length was defined by the MACS.

    TF analysis

    TF peaks were called using MACS 1.4 (default options) from six TFs plus two epigenetic regulators (EZH2 and EP300) at the four time points. TF peaks at two time points were considered overlapping if their summits were closer than 1 kb, as evaluated by HOMER MergePeaks. TF-bound regions were clustered by the temporal pattern of their peak overlaps across the four time points.

    Motif scanning was performed with HOMER (findMotifs Genome.pl hg19 –size given –mask –mknown <motif file>) using a motif list containing the default vertebrate motifs included with HOMER plus 618 recently defined heterodimer motifs (Jolma et al. 2015). Motifs with a P-value ≤10−15 were considered to be enriched. The redundant motifs were calculated and removed by STAMP (http://www.benoslab.pitt.edu/stamp/).

    Target genes were assigned to TF peaks by GREAT using 2-, 10-, 50-, and 100-kb distance thresholds. The enrichment of DEG genes among TF target genes was examined using Fisher's exact test. The analogous analysis was performed to evaluate the correlation between DRE clusters with DEG clusters.

    Global regulatory network construction

    The nodes of the gene regulatory network were seeded with all DEGs plus the expressed TFs, including MRs defined in MARINa. Putative MRs were colored red, whereas the rest were colored blue. Node radius was defined by the number of connections it has to other nodes in the network. Edges connecting nodes originated from three sources. First, edges were drawn between the eight transcriptional regulators (six TFs plus EZH2 and EP300) and their target genes. Second, edges were drawn between TFs and target genes based on the enrichment of the TF cognate motif within the target genes’ active distal enhancers. In these first two classes, edges between any two nodes whose absolute expression correlation was lower than 0.8 were eliminated, and the remaining were colored red (positive correlation) or green (negative correlation). Third, edges were drawn based on protein–protein interactions from the STRING database (Szklarczyk et al. 2017) and colored gray.

    Master regulator analysis

    We used MARINa (http://wiki.c2b2.columbia.edu/workbench/ index.php/MRA-FET), a master regulatory analysis algorithm (Lefebvre et al. 2010), to identify MRs within the gene regulatory network that governs VEGFA responses. The input for MRA-FET analysis is gene expression, a list of DEGs, and dynamic AE-target gene pairs (defined by GREAT). The enrichment threshold (FET P-value) was set to 0.01, and the P-value was further adjusted using the Bonferroni correction. The analysis was run in performed using a two-run FET calculation, which allows the differential activity of each TF to be examined. The output of MRA-FET is a list of potential MRs, each with a P-value, mode of activity, and other regulatory information.

    Previous SectionNext Section

    Data access

    All raw and processed sequencing data in this study have been submitted to the NCBI Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) under accession number GSE109626.

    Previous SectionNext Section

    Acknowledgments

    We thank Dr. Emery H. Bresnick (University of Wisconsin-Madison) for providing the GATA2 ChIP-seq antibody and Dr. Ruei-zeng Lin and Dr. Juan M. Melero-Martin for providing the ECFC cells and the training for in vivo Matrigel plug assays. G.C.Y. is supported by National Institutes of Health R01HG009663 and R01HL119099. W.T.P. is supported by an American Heart Association established investigator award, NIH 2UM1HL098166, and R01 HL138571 and by charitable donations to the Boston Children's Hospital Department of Cardiology. B.Z. is funded by the National Natural Science Foundation of China (91539109, 31671503), a Thousand Young Talent Award (16Z127060017), the Innovation program of Shanghai Municipal Education Commission (2017-01-07-00-01-E0028), National Key Research and Development Program of China (2018YFC1312702, 2018YFC1002400), and an AHA Scientist Development Grant (14SDG20380866).

    Author contributions: S.Y.W., X.D.L., F.Z., P.Y.Y., H.J.Y., S.M.S., Z.X.L., H.Y.L., and W.W.T. performed the experiments and analyzed the data. J.H.C. and S.P.G. performed computational analysis; J.H.C., G.C.Y., W.T.P., and B.Z. designed the study and wrote the manuscript. D.S.D., W.J.F., Z.G.H., X.S.L., and P.J.P. assisted with computational analysis and data processing. Y.Z. assisted with proliferation assay and ChIP-seq. Z.J.L. and K.S. were involved in the study design and manuscript revision. G.C.Y., W.T.P., and B.Z. oversaw all data analysis and study design.

    Previous SectionNext Section

    Footnotes

    • [Supplemental material is available for this article.]

    • Article published online before print. Article, supplemental material, and publication date are at http://www.genome.org/cgi/doi/10.1101/gr.239053.118.

    • Freely available online through the Genome Research Open Access option.

    • Received May 2, 2018.
    • Accepted December 12, 2018.

    This article, published in Genome Research, is available under a Creative Commons License (Attribution-NonCommercial 4.0 International), as described at http://creativecommons.org/licenses/by-nc/4.0/.

    Previous Section
     

    References

    1. ↵
      Bernstein BE, Mikkelsen TS, Xie X, Kamal M, Huebert DJ, Cuff J, Fry B, Meissner A, Wernig M, Plath K, et al. 2006. A bivalent chromatin structure marks key developmental genes in embryonic stem cells. Cell 125: 315–326. doi:10.1016/j.cell.2006.02.041
      CrossRefMedlineGoogle Scholar
    2. ↵
      Birdsey GM, Dryden NH, Shah AV, Hannah R, Hall MD, Haskard DO, Parsons M, Mason JC, Zvelebil M, Gottgens B, et al. 2012. The transcription factor Erg regulates expression of histone deacetylase 6 and multiple pathways involved in endothelial cell migration and angiogenesis. Blood 119: 894–903. doi:10.1182/blood-2011-04-350025
      Abstract/FREE Full Text
    3. ↵
      Carmeliet P. 2005. Angiogenesis in life, disease and medicine. Nature 438: 932–936. doi:10.1038/nature04478
      CrossRefMedlineGoogle Scholar
    4. ↵
      Chang S, Young BD, Li S, Qi X, Richardson JA, Olson EN. 2006. Histone deacetylase 7 maintains vascular integrity by repressing matrix metalloproteinase 10. Cell 126: 321–334. doi:10.1016/j.cell.2006.05.040
      CrossRefMedlineGoogle Scholar
    5. ↵
      Chen J, Fu Y, Day DS, Sun Y, Wang S, Liang X, Gu F, Zhang F, Stevens SM, Zhou P, et al. 2017. VEGF amplifies transcription through ETS1 acetylation to enable angiogenesis. Nat Commun 8: 383. doi:10.1038/s41467-017-00405-x
      CrossRefGoogle Scholar
    6. ↵
      De Val S, Black BL. 2009. Transcriptional control of endothelial cell development. Dev Cell 16: 180–195. doi:10.1016/j.devcel.2009.01.014
      CrossRefMedlineGoogle Scholar
    7. ↵
      De Val S, Anderson JP, Heidt AB, Khiem D, Xu SM, Black BL. 2004. Mef2c is activated directly by Ets transcription factors through an evolutionarily conserved endothelial cell-specific enhancer. Dev Biol 275: 424–434. doi:10.1016/j.ydbio.2004.08.016
      CrossRefMedlineGoogle Scholar
    8. ↵
      De Val S, Chi NC, Meadows SM, Minovitsky S, Anderson JP, Harris IS, Ehlers ML, Agarwal P, Visel A, Xu SM, et al. 2008. Combinatorial regulation of endothelial gene expression by Ets and forkhead transcription factors. Cell 135: 1053–1064. doi:10.1016/j.cell.2008.10.049
      CrossRefMedlineGoogle Scholar
    9. ↵
      DeNicola GM, Karreth FA, Humpton TJ, Gopinathan A, Wei C, Frese K, Mangal D, Yu KH, Yeo CJ, Calhoun ES, et al. 2011. Oncogene-induced Nrf2 transcription promotes ROS detoxification and tumorigenesis. Nature 475: 106–109. doi:10.1038/nature10189
      CrossRefMedlineGoogle Scholar
    10. ↵
      The ENCODE Project Consortium. 2012. An integrated encyclopedia of DNA elements in the human genome. Nature 489: 57–74. doi:10.1038/nature11247
      CrossRefMedlineGoogle Scholar
    11. ↵
      Goldberg AD, Allis CD, Bernstein E. 2007. Epigenetics: a landscape takes shape. Cell 128: 635–638. doi:10.1016/j.cell.2007.02.006
      CrossRefMedlineGoogle Scholar
    12. ↵
      Hake SB, Xiao A, Allis CD. 2004. Linking the epigenetic ‘language’ of covalent histone modifications to cancer. Br J Cancer 90: 761–769. doi:10.1038/sj.bjc.6601575
      CrossRefMedlineGoogle Scholar
    13. ↵
      Hasan SS, Tsaryk R, Lange M, Wisniewski L, Moore JC, Lawson ND, Wojciechowska K, Schnittler H, Siekmann AF. 2017. Endothelial Notch signalling limits angiogenesis via control of artery formation. Nat Cell Biol 19: 928–940. doi:10.1038/ncb3574
      CrossRefGoogle Scholar
    14. ↵
      Heintzman ND, Stuart RK, Hon G, Fu Y, Ching CW, Hawkins RD, Barrera LO, Van Calcar S, Qu C, Ching KA, et al. 2007. Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. Nat Genet 39: 311–318. doi:10.1038/ng1966
      CrossRefMedlineGoogle Scholar
    15. ↵
      Heinz S, Benner C, Spann N, Bertolino E, Lin YC, Laslo P, Cheng JX, Murre C, Singh H, Glass CK. 2010. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol Cell 38: 576–589. doi:10.1016/j.molcel.2010.05.004
      CrossRefMedlineGoogle Scholar
    16. ↵
      Hnisz D, Abraham BJ, Lee TI, Lau A, Saint-André V, Sigova AA, Hoke HA, Young RA. 2013. Super-enhancers in the control of cell identity and disease. Cell 155: 934–947. doi:10.1016/j.cell.2013.09.053
      CrossRefMedlineGoogle Scholar
    17. ↵
      Igarashi K, Kataoka K, Itoh K, Hayashi N, Nishizawa M, Yamamoto M. 1994. Regulation of transcription by dimerization of erythroid factor NF-E2 p45 with small Maf proteins. Nature 367: 568–572. doi:10.1038/367568a0
      CrossRefMedlineGoogle Scholar
    18. ↵
      Jeong HW, Hernandez-Rodriguez B, Kim J, Kim KP, Enriquez-Gasca R, Yoon J, Adams S, Scholer HR, Vaquerizas JM, Adams RH. 2017. Transcriptional regulation of endothelial cell behavior during sprouting angiogenesis. Nat Commun 8: 726. doi:10.1038/s41467-017-00738-7
      CrossRefGoogle Scholar
    19. ↵
      Jiang L, Yin M, Wei X, Liu J, Wang X, Niu C, Kang X, Xu J, Zhou Z, Sun S, et al. 2015. Bach1 represses Wnt/β-catenin signaling and angiogenesis. Circ Res 117: 364–375. doi:10.1161/CIRCRESAHA.115.306829
      Abstract/FREE Full Text
    20. ↵
      Jolma A, Yin Y, Nitta KR, Dave K, Popov A, Taipale M, Enge M, Kivioja T, Morgunova E, Taipale J. 2015. DNA-dependent formation of transcription factor pairs alters their binding specificity. Nature 527: 384–388. doi:10.1038/nature15518
      CrossRefMedlineGoogle Scholar
    21. ↵
      Kao HY, Ordentlich P, Koyano-Nakagawa N, Tang Z, Downes M, Kintner CR, Evans RM, Kadesch T. 1998. A histone deacetylase corepressor complex regulates the Notch signal transduction pathway. Genes Dev 12: 2269–2277. doi:10.1101/gad.12.15.2269
      Abstract/FREE Full Text
    22. ↵
      Katsuoka F, Yamamoto M. 2016. Small Maf proteins (MafF, MafG, MafK): history, structure and function. Gene 586: 197–205. doi:10.1016/j.gene.2016.03.058
      CrossRefGoogle Scholar
    23. ↵
      Krzywinski M, Schein J, Birol İ, Connors J, Gascoyne R, Horsman D, Jones SJ, Marra MA. 2009. Circos: an information aesthetic for comparative genomics. Genome Res 19: 1639–1645. doi:10.1101/gr.092759.109
      Abstract/FREE Full Text
    24. ↵
      Kumar MM, Goyal R. 2017. LncRNA as a therapeutic target for angiogenesis. Curr Top Med Chem 17: 1750–1757. doi:10.2174/1568026617666161116144744
      CrossRefGoogle Scholar
    25. ↵
      Langmead B, Salzberg SL. 2012. Fast gapped-read alignment with Bowtie 2. Nat Methods 9: 357–359. doi:10.1038/nmeth.1923
      CrossRefMedlineGoogle Scholar
    26. ↵
      Lefebvre C, Rajbhandari P, Alvarez MJ, Bandaru P, Lim WK, Sato M, Wang K, Sumazin P, Kustagi M, Bisikirska BC, et al. 2010. A human B-cell interactome identifies MYB and FOXM1 as master regulators of proliferation in germinal centers. Mol Syst Biol 6: 377. doi:10.1038/msb.2010.31
      Abstract/FREE Full Text
    27. ↵
      Leung DW, Cachianes G, Kuang WJ, Goeddel DV, Ferrara N. 1989. Vascular endothelial growth factor is a secreted angiogenic mitogen. Science 246: 1306–1309. doi:10.1126/science.2479986
      Abstract/FREE Full Text
    28. ↵
      Lu C, Han HD, Mangala LS, Ali-Fehmi R, Newton CS, Ozbun L, Armaiz-Pena GN, Hu W, Stone RL, Munkarah A, et al. 2010. Regulation of tumor angiogenesis by EZH2. Cancer Cell 18: 185–197. doi:10.1016/j.ccr.2010.06.016
      CrossRefMedlineGoogle Scholar
    29. ↵
      McLean CY, Bristor D, Hiller M, Clarke SL, Schaar BT, Lowe CB, Wenger AM, Bejerano G. 2010. GREAT improves functional interpretation of cis-regulatory regions. Nat Biotechnol 28: 495–501. doi:10.1038/nbt.1630
      CrossRefMedlineGoogle Scholar
    30. ↵
      Michalik KM, You X, Manavski Y, Doddaballapur A, Zornig M, Braun T, John D, Ponomareva Y, Chen W, Uchida S, et al. 2014. Long noncoding RNA MALAT1 regulates endothelial cell function and vessel growth. Circ Res 114: 1389–1397. doi:10.1161/CIRCRESAHA.114.303265
      Abstract/FREE Full Text
    31. ↵
      Obara N, Suzuki N, Kim K, Nagasawa T, Imagawa S, Yamamoto M. 2007. Repression via the GATA box is essential for tissue-specific erythropoietin gene expression. Blood 111: 5223–5232. doi:10.1182/blood-2007-10-115857
      CrossRefGoogle Scholar
    32. ↵
      Onodera K, Shavit JA, Motohashi H, Yamamoto M, Engel JD. 2000. Perinatal synthetic lethality and hematopoietic defects in compound mafG::mafK mutant mice. EMBO J 19: 1335–1345. doi:10.1093/emboj/19.6.1335
      Abstract/FREE Full Text
    33. ↵
      Ostuni R, Piccolo V, Barozzi I, Polletti S, Termanini A, Bonifacio S, Curina A, Prosperini E, Ghisletti S, Natoli G. 2013. Latent enhancers activated by stimulation in differentiated cells. Cell 152: 157–171. doi:10.1016/j.cell.2012.12.018
      CrossRefMedlineGoogle Scholar
    34. ↵
      Prasitsak T, Nandar M, Okuhara S, Ichinose S, Ota MS, Iseki S. 2015. Foxc1 is required for early stage telencephalic vascular development. Dev Dyn 244: 703–711. doi:10.1002/dvdy.24269
      CrossRefGoogle Scholar
    35. ↵
      Rafii S, Butler JM, Ding BS. 2016. Angiocrine functions of organ-specific endothelial cells. Nature 529: 316–325. doi:10.1038/nature17040
      CrossRefMedlineGoogle Scholar
    36. ↵
      Ramirez F, Ryan DP, Gruning B, Bhardwaj V, Kilpert F, Richter AS, Heyne S, Dundar F, Manke T. 2016. deepTools2: a next generation web server for deep-sequencing data analysis. Nucleic Acids Res 44: W160–W165. doi:10.1093/nar/gkw257
      CrossRefMedlineGoogle Scholar
    37. ↵
      Sacilotto N, Chouliaras KM, Nikitenko LL, Lu YW, Fritzsche M, Wallace MD, Nornes S, Garcia-Moreno F, Payne S, Bridges E, et al. 2016. MEF2 transcription factors are key regulators of sprouting angiogenesis. Genes Dev 30: 2297–2309. doi:10.1101/gad.290619.116
      Abstract/FREE Full Text
    38. ↵
      Shen L, Shao N, Liu X, Nestler E. 2014. ngs.plot: quick mining and visualization of next-generation sequencing data by integrating genomic databases. BMC Genomics 15: 284. doi:10.1186/1471-2164-15-284
      CrossRefMedlineGoogle Scholar
    39. ↵
      Sigova AA, Mullen AC, Molinie B, Gupta S, Orlando DA, Guenther MG, Almada AE, Lin C, Sharp PA, Giallourakis CC, et al. 2013. Divergent transcription of long noncoding RNA/mRNA gene pairs in embryonic stem cells. Proc Natl Acad Sci 110: 2876–2881. doi:10.1073/pnas.1221904110
      Abstract/FREE Full Text
    40. ↵
      Strahl BD, Allis CD. 2000. The language of covalent histone modifications. Nature 403: 41–45. doi:10.1038/47412
      CrossRefMedlineGoogle Scholar
    41. ↵
      Sun J, Hoshino H, Takaku K, Nakajima O, Muto A, Suzuki H, Tashiro S, Takahashi S, Shibahara S, Alam J, et al. 2002. Hemoprotein Bach1 regulates enhancer availability of heme oxygenase-1 gene. EMBO J 21: 5216–5224. doi:10.1093/emboj/cdf516
      Abstract
    42. ↵
      Szklarczyk D, Morris JH, Cook H, Kuhn M, Wyder S, Simonovic M, Santos A, Doncheva NT, Roth A, Bork P, et al. 2017. The STRING database in 2017: quality-controlled protein–protein association networks, made broadly accessible. Nucleic Acids Res 45: D362–D368. doi:10.1093/nar/gkw937
      CrossRefMedlineGoogle Scholar
    43. ↵
      Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. 2012. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 7: 562–578. doi:10.1038/nprot.2012.016
      CrossRefMedlineGoogle Scholar
    44. ↵
      Visel A, Blow MJ, Li Z, Zhang T, Akiyama JA, Holt A, Plajzer-Frick I, Shoukry M, Wright C, Chen F, et al. 2009. ChIP-seq accurately predicts tissue-specific activity of enhancers. Nature 457: 854–858. doi:10.1038/nature07730
      CrossRefMedlineGoogle Scholar
    45. ↵
      Voellenkle C, Garcia-Manteiga JM, Pedrotti S, Perfetti A, De Toma I, Da Silva D, Maimone B, Greco S, Fasanaro P, Creo P, et al. 2016. Implication of Long noncoding RNAs in the endothelial cell response to hypoxia revealed by RNA-sequencing. Sci Rep 6: 24141. doi:10.1038/srep24141
      CrossRefGoogle Scholar
    46. ↵
      Whyte WA, Orlando DA, Hnisz D, Abraham BJ, Lin CY, Kagey MH, Rahl PB, Lee TI, Young RA. 2013. Master transcription factors and mediator establish super-enhancers at key cell identity genes. Cell 153: 307–319. doi:10.1016/j.cell.2013.03.035
      CrossRefMedlineGoogle Scholar
    47. ↵
      Wilhelm K, Happel K, Eelen G, Schoors S, Oellerich MF, Lim R, Zimmermann B, Aspalter IM, Franco CA, Boettger T, et al. 2016. FOXO1 couples metabolic activity and growth state in the vascular endothelium. Nature 529: 216–220. doi:10.1038/nature16498
      CrossRefMedlineGoogle Scholar
    48. ↵
      Yao TP, Oh SP, Fuchs M, Zhou ND, Ch'ng LE, Newsome D, Bronson RT, Li E, Livingston DM, Eckner R. 1998. Gene dosage–dependent embryonic development and proliferation defects in mice lacking the transcriptional integrator p300. Cell 93: 361–372. doi:10.1016/S0092-8674(00)81165-4
      CrossRefMedlineGoogle Scholar
    49. ↵
      Yu H, Xue Y, Wang P, Liu X, Ma J, Zheng J, Li Z, Li Z, Cai H, Liu Y. 2017a. Knockdown of long non-coding RNA XIST increases blood–tumor barrier permeability and inhibits glioma angiogenesis by targeting miR-137. Oncogenesis 6: e303. doi:10.1038/oncsis.2017.7
      CrossRefGoogle Scholar
    50. ↵
      Yu W, Zhang F, Wang S, Fu Y, Chen J, Liang X, Le H, Pu WT, Zhang B. 2017b. Depletion of polycomb repressive complex 2 core component EED impairs fetal hematopoiesis. Cell Death Dis 8: e2744. doi:10.1038/cddis.2017.163
      CrossRefGoogle Scholar
    51. ↵
      Zhang CL, McKinsey TA, Chang S, Antos CL, Hill JA, Olson EN. 2002. Class II histone deacetylases act as signal-responsive repressors of cardiac hypertrophy. Cell 110: 479–488. doi:10.1016/S0092-8674(02)00861-9
      CrossRefMedlineGoogle Scholar
    52. ↵
      Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, Nusbaum C, Myers RM, Brown M, Li W, et al. 2008. Model-based Analysis of ChIP-Seq (MACS). Genome Biol 9: R137. doi:10.1186/gb-2008-9-9-r137
      CrossRefMedlineGoogle Scholar
    53. ↵
      Zhang B, Dietrich UM, Geng JG, Bicknell R, Esko JD, Wang L. 2009. Repulsive axon guidance molecule Slit3 is a novel angiogenic factor. Blood 114: 4300–4309. doi:10.1182/blood-2008-12-193326
      Abstract/FREE Full Text
    54. ↵
      Zhang B, Day DS, Ho JW, Song L, Cao J, Christodoulou D, Seidman JG, Crawford GE, Park PJ, Pu WT. 2013. A dynamic H3K27ac signature identifies VEGFA-stimulated endothelial enhancers and requires EP300 activity. Genome Res 23: 917–927. doi:10.1101/gr.149674.112
      Abstract/FREE Full Text
    55. ↵
      Zhou P, Gu F, Zhang L, Akerberg BN, Ma Q, Li K, He A, Lin Z, Stevens SM, Zhou B, et al. 2017. Mapping cell type-specific transcriptional enhancers using high affinity, lineage-specific Ep300 bioChIP-seq. eLife 6: e22039. doi:10.7554/eLife.22039
      CrossRefGoogle Scholar
    | Table of Contents
    OPEN ACCESS ARTICLE

    This Article

    1. Published in Advance January 22, 2019, doi: 10.1101/gr.239053.118 Genome Res. 29: 193-207 © 2019 Wang et al.; Published by Cold Spring Harbor Laboratory Press

    Article Category

    Share

    Preprint Server



    Navigate This Article

    Current Issue

    1. March 2026, 36 (3)
    1. Current Issue
    1. Alert me to new issues of Genome Research