Comprehensive transcriptomic profiling reveals SOX7 as an early regulator of angiogenesis in hypoxic human endothelial cells

Endothelial cells (ECs) lining the vasculature of vertebrates respond to low oxygen (hypoxia) by maintaining vascular homeostasis and initiating adaptive growth of new vasculature through angiogenesis. Previous studies have uncovered the molecular underpinnings of the hypoxic response in ECs; however, there is a need for comprehensive temporal analysis of the transcriptome during hypoxia. Here, we sought to investigate the early transcriptional programs of hypoxic ECs by using RNA-Seq of primary cultured human umbilical vein ECs exposed to progressively increasing severity and duration of hypoxia. We observed that hypoxia modulates the expression levels of approximately one-third of the EC transcriptome. Intriguingly, expression of the gene encoding the developmental transcription factor SOX7 (SRY-box transcription factor 7) rapidly and transiently increased during hypoxia. Transcriptomic and functional analyses of ECs following SOX7 depletion established its critical role in regulating hypoxia-induced angiogenesis. We also observed that depletion of the hypoxia-inducible factor (HIF) genes, HIF1A (encoding HIF-1α) and endothelial PAS domain protein 1 (EPAS1 encoding HIF-2α), inhibited both distinct and overlapping transcriptional programs. Our results indicated a role for HIF-1α in down-regulating mitochondrial metabolism while concomitantly up-regulating glycolytic genes, whereas HIF-2α primarily up-regulated the angiogenesis transcriptional program. These results identify the concentration and time dependence of the endothelial transcriptomic response to hypoxia and an early key role for SOX7 in mediating angiogenesis.

The ability of cells to survive and adapt to low-oxygen conditions is one of the most ancient and complex cellular processes and a critical feature for survival of multicellular organisms. The expression of thousands of genes was shown to be responsive to low oxygen levels, or hypoxia, but only a fraction of these genes are expressed within any given cell (reviewed in Ref. 1). Consequently, the hypoxic response varies significantly depending on the cell type and other factors, such as the level of hypoxia and its duration (2,3). The vast majority of hypoxia up-regulated gene expression changes in cells have been attributed to the hypoxia-inducible factor (HIF) 4 proteins, HIF-1␣ and HIF-2␣. These proteins are degraded in the proteasome during normoxic conditions, whereas they form heterodimeric complexes with HIF-1␣ protein during hypoxic stress. The HIF-␣-HIF-␤ complex binds to conserved hypoxia-response elements (consensus sequence 5Ј-RCGTG-3Ј) in promoters and enhancers of multiple genes and drives the gene expression regulating a wide range of functions including glycolytic metabolism and angiogenesis (reviewed in Ref. 4). However, the presence of consensus HIF-binding sites in gene promoter regions cannot accurately predict HIF-driven transcriptional responses because they are required but not sufficient for activity (5).
Endothelial cells (ECs) play important physiological roles in tissue development and disease progression, and these roles critically depend on their ability to sense and adapt to oxygen availability and generate new blood vessels. Thus, understanding the hypoxic controls of ECs and mechanisms of hypoxiainduced angiogenesis is important for emerging therapeutics such as specific HIF-2␣ inhibitors for renal cancer (6,7). Earlier studies uncovered hypoxic transcriptional programs related to glycolysis, angiogenesis, metastasis, and proliferation (reviewed in Ref. 8). Studies more recently showed that regulation of programs, such as glycolysis via the enzyme PFKFB3, can drive angiogenic EC behavior (9,10). Current thinking suggests that these programs are potentially dependent on the suite of active and available transcription factors in ECs, although recent research has also addressed the importance of chromatin remodeling, as defined by H3K4me, H3K27Ac, and RNA polymerase II occupancy, in determining HIF-regulated gene expression (3,11). Transcriptional studies of hypoxic ECs were conducted with early gene expression arrays (2,12), with microarrays with time series (13), and using CoCl 2 as a hypoxia mimic (14). However, the scope of these studies was limited in the numbers of genes detected, in the resolution and duration of time series, and by the inadequacy of CoCl 2 to faithfully recapitulate the hypoxic transcriptional response (15,16).
Although it is clear that the HIF-␣ subunits have redundant roles in regulation of some genes, it has become equally evident that the HIF-␣ subunits can differentially regulate genes with roles in development, disease, and metabolism (17)(18)(19). Thus, distinguishing HIF-1␣ from HIF-2␣ functions is an area of intense interest because these roles have manifested as differential regulation of intravasation/extravasation in tumor metastasis (20), the identification of HIF-2␣-independent clear cell kidney tumors (6), and their distinct roles in vascular disease (see also reviews in Refs. 1 and 21). The reasons for these differences are complex because HIF-1␣ and HIF-2␣ possess numerous differences in post-translational modifications, expression levels, and cofactor binding (reviewed in Ref. 8).
Here we describe detailed transcriptional characterization of the hypoxic EC response via RNA-Seq under multiple oxygen levels from 0.1% to normoxic (21%) conditions and from 1 to 48 h duration of hypoxia. We demonstrated that most gene activity is initiated within 8 h of hypoxic exposure by examining the time when individual genes first exhibit changes in expression. However, the point at which genes reach maximum expression varied considerably. Using these comparisons, we uncovered the previously unrecognized early and transient upregulation of the developmental angiogenesis regulator SRYrelated high-mobility-group box gene SOX7 by hypoxia. Gene depletion experiments showed that SOX7 is required for hypoxia angiogenic expression programs. Our broad survey revealed 1,889 up-regulated and 2,142 down-regulated genes that encompassed ϳ35% of all detected genes in HUVECs, a larger number than previously observed using CoCl 2 to mimic hypoxia (14), microarrays with only single hypoxia oxygen levels (2,12), and cancer cell lines such as MCF7 (22). Lastly, using HIF-1␣, HIF-2␣, and combined knockdowns in normoxia and hypoxia, we ascribed gene expression changes to HIF-1␣, HIF-2␣, either HIF isoform, or both HIFs. We identified SOX7 as a rapid HIF-2␣-regulated gene and decomposed the hypoxic response into HIF-specific programs that enhance understanding of the complex interplay of HIFs in homeostatically regulating endothelial hypoxic response.

A survey of the hypoxic transcriptional response of ECs
A growing body of research showed that signaling initiated by hypoxia is cell type-specific, time-dependent, and exquisitely complex (2,8,12,23,24). To better understand the genetic programs of hypoxia in ECs, we sought to precisely define the hypoxic transcriptional response in the context of a widely used EC model system, primary cultured HUVECs. Briefly, we conducted stranded RNA-Seq of low-passage HUVECs subjected to five different oxygen levels, ranging from 0.1 to 21%, and 10 different time points, from 1 to 48 h of hypoxia (see "Experimental procedures"). We used the normoxic condition (21% O 2 ) at these time points to control for potentially confounded gene expression changes caused by cell cycling, serum, or other time-dependent influences.
We sequenced to a depth of 16.7 Ϯ 1.1 million reads, which provided a sufficient balance for gene detection, sampling density, and ability to detect potential PCR biases. We detected expression of 11,469 protein-encoding genes and 1,635 noncoding genes (Ͼ30 reads in at least one sample and Ͼ100 reads across all samples; human genome version GRCh38.87). We observed that although the feature detection rate correlated with sequencing depth across samples (from 10 to 26 million reads), the average increase was at a rate of approximately three features per additional 1 million reads, suggesting that this design was sufficiently robust to avoid missing transcripts with moderate expression (Fig. S1).
We detected expression changes within the first 1-2 h of hypoxia exposure in several genes (Table 1). This included many established HIF target genes such as BHLHE40 (basic helix-loop-helix family member e40), known for its role as a hypoxia-regulated EMT gene; HK2 (hexokinase 2), an initiator role in glucose metabolism; and ADM (adrenomedullin), a potent hypotensive peptide. Activation of genes with less wellknown roles in hypoxia was also observed, such as DUSP6 (dual specificity phosphatase 6; but see Ref. 12), a regulator of mitogen-activated protein kinase activity, as well as SOX7 (SRY-box 7), a developmental transcriptional regulator. SOX7 is required for vasculogenesis and angiogenesis during development (25, 26) but has not previously been shown to be important for the hypoxic response. This analysis also revealed an increase in the hypoxia-responsive long noncoding RNA (lncRNA), MIR210HG, earlier than 24 h as shown (27). In fact, transcriptional changes for a variety of genes under hypoxia were noticed at the earliest 1-h time point, rivaling the speed of some early cellular responses to hypoxia, including changes in mitochon-

SOX7 regulates hypoxia-dependent angiogenesis
drial respiratory activity (28). In contrast, the well-characterized HIF-1␣ target gene encoding vascular endothelial growth factor A, VEGFA, increased gradually in expression during 48 h of hypoxia (Fig. 1A). Although transcript decay may play a role in gene expression levels during hypoxia, our examination of recent work with HUVECs supports a strong concordance between nascent transcripts (labeled with 4-thiouridine) as compared with total transcripts for genes with large expression changes in hypoxia (29) (GSE89831; Fig. S2). This finding shows that at least for the most differentially expressed transcripts, our approach approximates changes in transcription rather than degradation. We identified 4,023 protein-coding and 284 noncoding genes (35 and 17% of detected genes, respectively) that showed hypoxia-responsive expression changes during 48 h of hypoxia exposure (FDR-adjusted p Ͻ 0.01). Based on the most significant changes within two consecutive time points, this corresponded to 1,884 up-regulated and 2,139 down-regulated protein-coding genes along with 176 up-regulated and 108 down-regulated ncRNA transcripts (Fig. 1, B and C, and Table  S1). The greater number of differentially regulated genes reported here is greater than those reported for HUVECs or other cell types under hypoxic stress (2,12,14,22). It is likely due to a combination of factors, including the detection of transient changes and transcripts more amenable to stranded library RNA-Seq than microarrays. Although we cannot rule out that apoptotic expression programs may have contributed to this difference, others have shown that HUVECs are minimally apoptotic within 24 h of severe (Ͻ1% O 2 ) hypoxia (up to 3% terminal deoxynucleotidyl transferase dUTP nick end labeling-positive cells) (30). Moreover, we observed neither detached cells within 24 h of hypoxia nor up-regulation of genes in the apoptotic transcriptional signature, HALLMARK_APO-PTOSIS (MsigDBv6.2).
To resolve domains of expression activity for hypoxia-responsive genes, we compared expression profiles across the entire series of oxygen and time points. Briefly, we used affinity propagation clustering (31) of gene expression during 0.1% O 2 hypoxia to assign genes into groups based on similar changes in expression. We identified 29 clusters using this method, ranging in size from 40 to 291 genes, with a median of 151 genes (Fig.  1C). We conducted enrichment analysis for Gene Ontology (GO) terms (Molecular Function and Biological Process) for these clusters to better understand the unique domains of transcriptional activity. Given the established link between hypoxia and metabolism (reviewed in Ref. 4), many of the terms enriched in up-regulated genes were related to metabolic functions. Interestingly, some genes with transient increases in expression may be related to developmental and physiological processes, as is evident in the first three gene clusters (GO terms: multicellular organismal process, thyroid gland development, and embryonic morphogenesis). Moreover, a group of genes up-regulated at 48 h included those related to a stress response (GO term: ATF6-mediated unfolded protein response), as well as a metabolic change (GO term: cellular amino acid metabolic process), further supporting distinctions between the cellular response to acute (Ͻ24 h) and longer duration (Ͼ24 h) hypoxia.
Mitochondria play a central role in cellular metabolism, and hypoxia has been shown to ameliorate mitochondrial defects by altering metabolic pathways (32). As with other cell types, ECs down-regulated mitochondrial-encoded gene transcription in hypoxia (Fig. S3A). In addition to changes in glycolytic regulator gene expression, mitochondrial gene expression patterns were consistent with a metabolic shift; e.g. the mitochondrial fission gene, MTFP1, was increased concordantly with decreased expression of mitochondrial fusion genes, MFN1/2, during the first 5-6 h ( Fig. S3B) (33). We also noticed a concurrent increase in PDK1 expression during hypoxia, an enzyme driven by HIF-1␣ and necessary for the switch from oxidative phosphorylation to glycolysis (34) (reviewed in Ref. 35). Of note, the transient increase in MTFP1 and PDK1 also coincided with other metabolic genes such as HK2, mentioned above (Table S1).
ECs of different vessel types have both shared and unique roles in vascular function. Recent work has uncovered unique molecular attributes of different endothelial subtypes (36). To determine whether EC subtypes exhibited differences in their hypoxic transcriptional programs, we conducted additional RNA-Seq of primary cultured human pulmonary arterial endothelial cells (HPAECs) and human pulmonary microvascular endothelial cells (HPMVECs). We compared expression at 7 h of 1% hypoxia versus normoxic conditions and found that the changes in gene expression were in general agreement among the three EC subtypes, notably in the up-regulation of genes PPFIA4, AK4, EGLN3, PFKFB4, SLC2A1, and down-regulation of CYP1A1 and SPSB2 ( Fig. S4 and Table S2), many of which are known to play a role in metabolism. Although a few genes were found to have relatively subtype-specific increases with the 7-h hypoxia treatment, including FAM13A in HPMVECs and HUVECs and FAM162A in HPAECs and HUVECs, the strong similarities in hypoxic profiles indicated that the impact of hypoxia on HUVEC transcription is generally applicable to a variety of low passage primary EC subtype lines.

Temporal activation of the hypoxic program
Few studies have examined transcriptional program activity of ECs during hypoxia simultaneously with levels of oxygen. Our goal was therefore to determine whether different transcriptional programs were engaged at different oxygen levels. In comparing expression changes in moderate hypoxia (5% O 2 ) to severe hypoxia (0.1% O 2 ), we found that differences were primarily in the degree of expression change rather than the specific genes involved ( Fig. 1C and Fig. S5, A and B). Although the protein levels and cellular activities may vary depending on thresholds of gene expression and other more complex feedback mechanisms, this observed gradation in transcriptional response suggested that similar cis-regulatory mechanisms are engaged across levels of hypoxia, such that HIFs exert their influences in a predictable, titratable manner.
To elucidate the dynamic controls of the EC hypoxic response, we assigned times to the onset of gene expression changes based on when genes in the hypoxic treatment (1% O 2 ) deviated from normoxic conditions. Briefly, activation (or inactivation) for each gene was determined by the first point at which the hypoxic expression Loess regression trend differed from normoxia (Fig. S6). This approach has the advantage that it does not distinguish continued expression changes or transient expression changes and is therefore more relevant for understanding the factors required for initiation of an expression change. Using these data to pinpoint times for the onset of expression changes of each gene, we found the expression of most hypoxia up-regulated genes was noticeably increased within 9 h of hypoxia (78%) with a mode of 3 h and most downregulated genes noticeably decreased within 9 h of hypoxia (69%) with a mode of 4 h ( Fig. 1D and Table S3). However, expression changes for many of these genes continued throughout the 48-h experiment. ϳ28 and 32% of genes increasing or decreasing, respectively, had either reached a maximal change at 48 h or were continuing to change at 48 h ( Fig. 1D and Table  S3). This observation suggests that differing factors drive initial activation versus sustained gene expression activity in hypoxia, possibly relating to availability of transcription binding partners and chromatin state.
Although the HIFs are primarily regulated via protein degradation in the proteasome, we observed distinct and complementary transcriptional activities where HIF1A expression declined over time and EPAS1 expression increased between 3 and 12 h (Fig. S7A). Interestingly, HIF3A, a potential negative transcriptional regulator of HIF-1␣ and HIF-2␣ targets, steadily increased in expression during 48 h of (0.1% O 2 ) hypoxia. We also found an almost 8-fold (log2) increase in expression of the lncRNA HIF1A-AS2 from ϳ3 to 24 h, in opposition to HIF1A expression (Fig. S7A). This observation may have been obscured in previous studies without temporal sampling or strand-specific RNA-Seq libraries. There is little

SOX7 regulates hypoxia-dependent angiogenesis
known about the HIF1A-AS2 lncRNA, although it is important for mesenchymal stem cell maintenance in brain tumors (37). We found that HIF1A-AS2 contains five consensus HIF-binding sites at its promoter, an indication that its expression may be driven by HIF (Fig. S7B). To better understand the relationship between HIF1A and HIF1A-AS2, we examined ChIP sequencing for HIF-1␣ and active chromatin marks H3K27Ac and H3K4me3 in HUVECs after 24 h of hypoxia (1% O 2 ) (GSE39089) (11); this corresponded to the time of low HIF1A expression and the highest HIF1A-AS2 expression. We observed that HIF-1␣ bound at the promoter of HIF1A-AS2 (adjusted p value ϭ 5e-6) but not its own promoter or that of HIF1A-AS1, whereas both regions had consensus HIF-binding sites and active chromatin marks (Fig. S7). Combined with our expression profiles, these observations suggest a feedback regulatory role for HIF1A-AS2 and/or EPAS1 in HIF1A expression as seen in periodontal ligament cells (38) and renal cell carcinoma and glioblastoma cells (39).

SOX7 is a key hypoxia-inducible angiogenic regulator
We found 100 transcription factors, 37 cotranscription factors, and 20 chromatin modifiers whose expression increased during hypoxia (of 1,515 detected from AnimalTFDB v2.0). Several of these gene expression regulators showed increasing trends in expression during the 48-h time course (Fig. 2A). Interestingly, the gene with relatively largest and earliest peak in expression was the SoxF family member transcription factor SOX7, showing an increase in expression within 1 h and a peak

SOX7 regulates hypoxia-dependent angiogenesis
in expression by 5 h of hypoxia (Fig. 2, A and B, and Table 1). A similar hypoxia-dependent change in expression was not observed for the other SoxF family members, SOX17 and SOX18, even though all three were found expressed in HUVECs at appreciable levels (SOX7 5.8 Ϯ 0.4, SOX17 4.2 Ϯ 0.6, and SOX18 6.2 Ϯ 0.4; mean RPKM Ϯ S.D.).
The SOX7 locus only contains two consensus HIF-binding sites within 2 kb of the TSS, whereas SOX17 and SOX18 contain three and four, respectively, supporting the observation that consensus HIF-binding sites are only one part of several regulatory mechanisms for HIFs (reviewed in Ref. (8). ETS family transcription factor, ETV2, was shown to directly regulate SOX7 expression during development (40) but was below the detection threshold in our RNA-Seq experiments. Although SOX7 was expressed at a higher level in hypoxic HUVECs (7.0 Ϯ 0.1) at 7 h (1% O 2 ), it was also expressed in the other two EC types studied: HPAEC, 6.0 Ϯ 0.7; and HPMVEC, 6.1 Ϯ 0.3. We confirmed that SOX7 protein also exhibited a transient increase at 8 h of hypoxia, similar to SOX7 mRNA that was decreased by 24 h (Fig. 2, C and D).
SOX7 has been shown to regulate important EC genes such as CDH5 (25); thus, we addressed how SOX7 contributed to the angiogenic transcriptional program of hypoxic ECs. We thus conducted SOX7 knockdowns in HUVECs and subjected them to normoxic or hypoxic (1% O 2 ) conditions followed by RNA-Seq. We confirmed loss of SOX7 expression at the 3-h time point (Fig. S8) and compared the expression profiles of SOX7 versus scramble siRNA-treated cells following 3-6 h of hypoxia. We identified several hypoxia-induced genes that were dependent on SOX7 expression (Fig. 3A and Table S4). We showed SOX7-dependent regulation of EC identity genes at 3-6 h of hypoxia, such as CDH5 (logFC 0.9, FDR-adjusted p ϭ 2e-7), PECAM1 (logFC 1.7, FDR-adjusted p ϭ 2e-15), FLT1 (logFC 1.0, FDR-adjusted p ϭ 6e-7), and EMCN (logFC 0.8, FDR-adjusted p ϭ 1e-5) but did not see changes in the expression of other EC genes such as TEK and PTPRB (both FDRadjusted p Ͼ 0.05). Two prominent hypoxia-induced angiogenesis related genes, NOTCH4 and DLL4, were also reduced in expression following SOX7 depletion ( Fig. 2A), consistent with a link between hypoxia, Notch signaling, and angiogenesis (41). Gene Ontology analysis of genes regulated by hypoxia within 3-7 h (͉logFC͉ Ͼ0.5, FDR-adjusted p Ͻ 0.05) and genes regulated by SOX7 during normoxia or 3-6 h of hypoxia (1% O 2 ) showed that SOX7 also induced hypoxia-dependent genes  Fig. 1C). Filled data points are those with significant expression changes (FDR-adjusted p Ͻ 0.05). Labeled points are genes implicated in angiogenesis, contained in one or more of the gene sets: HALLMARK_ANGIOGENESIS, GO_ANGIOGENESIS, GO_SPROUTING_ANGIOGENESIS, and GO_CELL_MIGRATION_IN-VOLVED_IN_SPROUTING_ANGIOGENESIS. B, intersection of differentially expressed genes from A followed by GO analyses. Hypoxia-dependent genes were limited to ͉log2 FC Ͼ 0.5͉ and FDR-adjusted p Ͻ 0.05, and SOX7-dependent genes were limited to log2 FC Ͼ 1 and FDR-adjusted p Ͻ 0.01. Unambiguous Gene Ontology terms were chosen for representation of full list (black text), and unlabeled regions had associations with FDR-adjusted p Ͼ 0.01, fewer than five genes in the gene set, or fewer than three DE genes (available as Table S5). C, representative images and quantification for tube length in a hypoxic tube-forming assays with SOX7 knockdown (Scr n ϭ 6, SOX7 siRNA n ϭ 5). Scale bar, 250 m. pval, p value.
To establish whether SOX7 influenced angiogenesis in hypoxic conditions, we used tube formation assays. SOX7 KD in HUVECs reduced tube length by 30% under hypoxic (1% O 2 ) conditions (Fig. 3C). However, neither vessel number nor branch number were changed (both p Ͼ 0.05). Thus, SOX7 was required for aspects of angiogenesis such as vessel elongation in hypoxic conditions, similar to its role in regulating vasculogenesis and angiogenesis during development (40).

HIF1 and HIF2 programs
There are over one million sequences matching the consensus HIF-binding site sequence RCGTG in the human genome with fewer than 1% of these with HIF binding in cells under hypoxia (42). Differences in HIF-1␣ and HIF-2␣ targets appear to be in part driven by cofactors (43,44), and thus discerning HIF-1␣-and HIF-2␣-dependent programs based on promoter motifs is not currently possible. Although ChIP-seq experiments shed light on differences between HIF-1␣-and HIF-2␣-driven hypoxia transcriptional programs in breast cancer cells (42,45), it is unknown whether these differences are present in ECs. Moreover, only ϳ8% of the genes bound by HIF-1␣ are similar between breast cancer cells and HUVECs in hypoxia (11). Recently, it was shown using CRISPR-Cas9 knockout of HIFs in a transformed renal cell line that HIF-1␣ and HIF-2␣ did not compete for binding sites in hypoxia, supporting the importance of establishing specific transcriptional activities for HIFs (46). Studies have identified gene expression changes shared between ECs with HIF-1␣ overexpression and 24 h of hypoxia (12), but the role of HIF-2␣-mediated regulation is not clear. Therefore, to understand the functional program differences driven by HIF-1␣ versus HIF-2␣ in ECs under hypoxia, we conducted HIF1A (HIF-1␣) and EPAS1 (HIF-2␣) knockdown experiments in HUVECs exposed to normoxic (21% O 2 ) and hypoxic (1% O 2 ) conditions followed by RNA-Seq.
Knockdown was efficient at both the protein and mRNA levels with either of two unique siRNAs targeting each HIF, with 86 -89% efficiency for HIF1A and 90 -91% efficiency for EPAS1 at the mRNA level (Fig. 4, A and B). Following RNA-Seq, we used RUVseq for normalization across treatments and Scr-siRNA controls. Differential expression comparisons between siRNA-treated and Scr-siRNA-treated cells showed 425 and 290 genes activated (down in siRNA-treated) and inactivated (up in siRNA-treated) by HIF-1␣ in hypoxia, respectively (FDR-adjusted p Ͻ 0.05 and log2 FC Ͼ 1). Similarly, we found 359 and 338 genes activated and inactivated by HIF-2␣ in hypoxia, respectively ( Table 2). We found that numerous glycolytic genes, such as PFKFB4, LDHA, and HK2, were up-regulated by HIF-1␣ and either not affected or down-regulated by HIF-2␣ ( Fig. 4C and Table S6). PFKFB4, the glycolytic metabolite regulator required for ATP generation via activation of PFK-1, is a rate-limiting step in glycolysis and enables cancer cell survival during hypoxia (47). We found that PFKFB4 expression in HUVECs was positively regulated by HIF-1␣ but negatively regulated by HIF-2␣. However, PFKFB3, a gene encoding a homologous protein that drives vessel sprouting (10), was not altered by HIF-2␣ levels and depended only on HIF-1␣ for its hypoxia-dependent expression increase. Some genes, such as PPARG and INSR, appeared to be dependent on both HIF-1␣ and HIF-2␣ for their expression in hypoxia, although this observation did not distinguish between the requirement for each independent HIF. We also identified numerous genes for which hypoxia-induced gene expression was abrogated only when both HIFs were depleted, indicating that either HIF isoform was sufficient for mediating hypoxiainduced transcriptomic changes. One such example is ANG-PTL4, a gene encoding a hypoxia-regulated protein; thus, these data extend prior observations that ANGPTL4 is regulated by HIF-1␣ (12) by demonstrating that it can also be regulated by HIF-2␣.
HIFs have many similar and opposing roles in a variety of contexts (reviewed in Ref. 8), likely because of their different binding sites across the genome (46). We compared GO annotations for genes that fell into categories relating to whether their hypoxic expression depended on HIF-1␣, HIF-2␣, both HIFs, or either HIF and whether they were also up-regulated in hypoxia (refer to "Experimental procedures" and Fig. 4D). HIF-1␣-dependent genes increasing in hypoxia were involved in glycolysis and pyruvate metabolism (Table S7 and Fig. 4D). In contrast, HIF-2␣-dependent genes increasing in hypoxia were involved with angiogenesis. Interestingly, HIF-2␣dependent genes that did not increase in hypoxia were associated with development (Fig. 4D). Although specific programs likely encompass numerous activities beyond a simple GO analysis, these comparisons show that during hypoxia, genes are mostly dependent on one HIF or other (HIF-1␣, n ϭ 122; HIF-2␣, n ϭ 73) rather than either HIF (n ϭ 22) or both HIFs (n ϭ 28).
Because SOX7 plays a critical angiogenic role during development and appears to regulate angiogenic programs in hypoxia, we addressed whether its expression depended on HIF-1␣, HIF-2␣, or both. Examination of the HIF KD RNA-Seq data showed as late as 7 h of hypoxia exposure that expression of SOX7 depended on HIF-2␣ (logFC 1.4, FDR-adjusted p ϭ 4e-8) but minimally on HIF-1␣ (logFC 0.4, FDR-adjusted p ϭ 8e-3). To determine whether HIF-2␣ interacted directly with the SOX7 promoter during hypoxia, we conducted ChIP-PCR for two regions of the SOX7 promoter with HIF-2␣ antibody. Despite large variability in enrichment, one region of the SOX7 5Ј-UTR was associated with HIF-2␣ during hypoxia (Fig. S9), consistent with the role of HIF-2␣ in regulating SOX7 expression in hypoxia.
A fundamental question regarding the hypoxic response of cells is whether and to what extent HIF-independent oxygensensing mechanisms drive transcriptional activity. Reports such as in Drosophila melanogaster showed that a homolog of the estrogen-related receptor acted both with and independently of the HIF homolog sima to drive metabolic changes (48). Thus, to identify HIF-1␣/HIF-2␣-independent changes during hypoxia, we analyzed genes up-regulated in hypoxia, as well as following combined HIF-1␣/HIF-2␣ knockdown. However, no genes were detected as being strongly up-regulated by hypoxia in both the presence and absence of HIFs (Fig. S10), indicating that HIF-independent transcriptional changes were not important in ECs.

Discussion
We showed by sampling across various time points and oxygen levels that the hypoxic transcriptional response of primary ECs was far greater than expected, engaging more than onethird (4,023/11,469) of all genes detected within RNA-Seq experiments. A major caveat in translating the transcriptional response in EC cultures to physiologically meaningful hypoxia is that oxygen tension varies across tissues of the body. Thus, ECs adapted to relatively hypoxic environments may be less sensitive to mild hypoxia as compared with other ECs. Hence there is a need to understand responses to both different levels and durations of hypoxia and different types of EC. Our observation that many gene expression changes occurred similarly in different ECs and in a gradual manner across oxygen levels from 0.1 to 21% (Fig. S5) justifies the value of primary EC culture as a meaningful system for delineating features of the hypoxic transcriptional response.
Hypoxia-driven gene expression changes observed in the RNA-Seq experiments may be result of either 1) an increasing proportion of cells exhibiting a large change or 2) all cells exhib- The scores are the product of RMS scaled significance (Ϫlog10 FDR-adjusted p) and log2 fold-change values for HIF1A and EPAS1 knockdown comparisons. D, heat map showing relative expression of genes that retained expression (within 90% of control) or decreased with HIF siRNA treatment (less than 70% of control). Genes are included based on average expression across siRNA versions and replicates and grouped according to whether they are increased in hypoxia (right side, logFC Ͼ 0, FDR-adjusted p Ͻ 0.01). Gene Ontology terms were assigned to groups using hypergeometric test for gene enrichment (top hit listed). Relative expression is thresholded for visualization ([-2,2] log2 FC). Refer to Table S7.

SOX7 regulates hypoxia-dependent angiogenesis
iting a small change. By utilizing primary cultures of HUVEC monolayers, we were more confident that the cells are of homogenous origin and were exposed to the same oxygen levels. This would argue that the observed changes occurred in every cell. With this caveat, hypoxia produced a graded transcriptional response all the way down to virtually anoxic condition (0.1% O 2 ). The implication is that this fine-tuned hypoxiaactivated transcriptional response enables cells to adapt certain traits, such as their metabolic programs or degree of angiogenic signaling, to a specific hypoxic condition. Such a model leads to implicit gradients of angiogenic factors, such as VEGF, that induce graded angiogenesis.
We observed that the EC hypoxic transcriptional program began within 1 h of exposure. We detected expression changes for most genes regulated by hypoxia within the first 8 h of exposure. Interestingly, many genes showed transient changes in expression (Fig. 1C) that would not have been evident without time-course sampling. Many of these transient gene expression changes were associated with developmental GO terms (Fig.  1C) and included the developmental angiogenesis regulator gene SOX7. Although Sox7 deletion in mice is embryonic lethal because of cardiovascular defects (49), a role for SOX7 in the cellular response to hypoxia has not been previously observed. Our finding that SOX7 expression increased abruptly in response to hypoxia and contributed to vessel growth as implicated by GO comparisons and tube-forming assays indicates a potential mechanism by which SOX7 controls angiogenesis during development under hypoxic condition. Moreover, our observation that the global hypoxic transcriptional response is similar between primary cultures of human arterial, venous, and microvascular endothelial cells (Fig. S4), combined with the similar expression levels for SOX7 in these cells, suggests SOX7 also functions in a nondevelopmental context such as after vascular injury.
Studies have identified important and overlapping roles for HIF-1␣ and HIF-2␣ in development (50,51) and disease (18,(52)(53)(54)(55)(56). Thus, understanding the individual roles of HIF-1␣ and HIF-2␣ in the vascular system is becoming increasingly important as we seek to pharmacodynamically modulate these transcriptional controls in pathological situations such as clear cell renal cell carcinoma with circulating HIF-2␣ antagonists (6,7). The individual contributions of HIF-1␣ and HIF-2␣ to the hypoxic control of ECs is poorly understood in part because of previous technological limitations. Here we showed that HIF-1␣ drives a shift to glycolytic metabolism, whereas HIF-2␣ appears to drive angiogenic programs, partly through SOX7.
Interestingly, we did not find that the other SoxF family members SOX17 and SOX18 shared the early-response hypoxic profile of SOX7. Several publications have highlighted the redundant roles of these three transcription factors in development (57)(58)(59), although whether they can substitute for one another as part of the hypoxic response is not yet known. Given that all three SoxF family members are expressed in HUVECs and contain consensus HIF-binding sequences in their promoters, the differences in cis-regulatory elements and cofactors could be critical for understanding why SOX7 is one of the earliest transcription factors up-regulated in response to hypoxia in ECs. Although it is unknown how many of the hypoxic changes are directly regulated by the HIFs, the large numbers of genes that are dependent on one or the other HIF suggest that binding partners and cotranscriptional activators are as important as the presence of hypoxia-response elements in shaping the hypoxic response.
Our analyses with stranded RNA-Seq support a further complexity of HIF regulation (reviewed in Ref. 8), whereby expression of each HIF is regulated in a time-dependent manner and the level of HIF1A mRNA is opposed by levels of HIF1A-AS2 lncRNA. Consistent with previous reports, we found ECs to show a variety of unique and shared roles for HIF-1␣ and HIF-2␣ in our gene knockdown comparisons. Many glycolytic gene expression changes, such as regulation of HK2, PFKFB3, PFKFB4, LDHA, ENO2, SLC1A1, etc., in hypoxia were found to be dependent on HIF-1␣ alone, whereas HIF-2␣ appeared to independently control a variety of genes pertaining to signaling and apoptotic processes.
Finally, the observation that most expression changes begin within 1 h to just a few hours following hypoxic exposure and that some of these changes are transient is not just an important caution for research. It suggests that the transcriptional programs of the cellular hypoxic response are in continual and rapid flux, enabling cells to rapidly adapt their metabolic programs to changing hypoxic conditions. The transient nature of some gene expression patterns (e.g. SOX7) suggests feedback mechanisms exist for some genes, whereas others appear to be unrestricted (e.g. VEGFA), enabling differential acute and chronic hypoxic responses.

Cell culture
All cell lines were obtained from Lonza (Walkersville, MD). Single donor primary HUVECs (P4, catalog no. C2517A) were grown in EGM-2 (catalog nos. CC-3156 and CC-4147). P6 HUVECs were split into 6-well plates and grown for 5 days with fresh medium every 24 h and then subjected to hypoxic or normoxic conditions. The medium was changed 1 h prior to the beginning of experiment and at 24 h for 48-h samples. Primary human pulmonary microvascular endothelial cells (HPMVECs, catalog no. CC-2527) were grown in EGM-2 (catalog no. CC-3156 and CC-4147) until P6. Primary human pulmonary arterial endothelial cells (HPAECs, catalog no. CC-2530) were grown in EGM-2 until P7. Hypoxia was controlled via the ProOx chamber (model 110; BioSpherix, Lacona, NY) and 10% FBS was added to all cells.

RNA-Seq
Cells in hypoxic and normoxic conditions were immediately rinsed with PBS and lysed with TRIzol (ThermoFisher, Waltham, MA) at specified time points. Total RNA was isolated, and samples were submitted to the University of Chicago Genomics Facility (Chicago, IL). RNA quality was assessed via Bioanalyzer (Agilent, Santa Clara, CA) and purified using the oligo(dT) method. cDNA was prepared with the TruSeq stranded mRNA library kit (Illumina, San Diego, CA), and samples were multiplexed for 50-bp single-end read sequencing on a HiSeq 2500 sequencer (Illumina, San Diego, CA).

RNA-Seq data pre-processing
R base (60) and Bioconductor (61) software packages were used for data processing. The TrimGalore (62) wrapper for Cutadapt (63) and FastQC (64) was used to remove adapters, low-quality sequences (Ͻ 28), and resulting reads shorter than 35 bp. Sequence reads were uniquely mapped to the hg38 human genome using the strand-specific seed-and-vote methodology implemented in the Rsubread package (65). File indexing and manipulations were conducted using Rsamtools (66). Feature annotations were established with GRCh38.87, and quantitation was conducted using reads fully overlapping with features via GenomicFeatures and GenomicAlignments (67). The large number of HUVEC samples enabled minimization of potential PCR duplication bias. This was accomplished by modeling read depth as a polynomial function of uniquely mapping reads for each gene. The read counts for each outlier gene (residual Ͼ 1.5 times the interquartile range above the upper quartile residual) were then iteratively shrunken toward the line of best fit for all samples. Alignment for SOX7 knockdown experiment data were conducted with Star (68) and quantitation of reads with Salmon (69).

RNA-Seq data modeling
ncRNA was considered for initial comparisons, although subsequent analyses were limited to protein coding genes. The genes were filtered such that they were represented by at least 100 read counts across all samples in the initial three sequencing batches and separately in the SOX7 KD experiments. Further filtering included the requirement of 30 read counts in at least one sample within each of the three initial sequencing batches and at least five samples with 20 read counts in the SOX7 KD experiments. EdgeR (70) was used for data handling and trimmed mean of M values normalization, as well as differential expression with the quasi-likelihood F test. Dispersion was estimated using log2 transformed time (h) and oxygen levels (%). Batch adjustments for normalized data to control for sequencing run batch and other artifacts were made using Combat in SVA (71) for most experiments and RUVSeq (72) for HIF gene KD experiments. SOX7 KD experiment samples were sequenced simultaneously, and therefore no batch correction was used.

Temporal RNA-Seq modeling
Loess regression was used to interpolate relative gene expression changes (from normoxia, 1 h) for the 0.1% O 2 hypoxic time series. One hundred time points were interpolated from the model whereby relative expression (cpm) was modeled as a function of time (h). Confidence bands (95%) were estimated from regression curves. Onset of activation was defined by three criteria: 1) the first significant difference, 2) in the direction of the largest change in the 48-h time series, and 3) persistence for at least 15% of the log2 time series in cases where it was a transient change (see RELA in Fig. S6).
To identify the strongest, specific changes associated with HIF1A and EPAS1 knockdowns, several criteria were used. For HIF1A, these genes had lower expression in the HIF1A and combo HIF1A/EPAS1 knockdowns than the control, had lower expression in the HIF1A than the EPAS1 knockdown (all logFC Ͻ Ϫ0.5, FDR-adjusted p Ͻ 0.05), and were not differentially expressed in EPAS1 knockdown versus control (͉logFC͉ Ͻ 0.5, FDR-adjusted p Ͼ 0.2). The reciprocal filtering was conducted for EPAS1. For genes that were dependent on both HIFs for their expression, the criteria of being decreased in both HIF1A and EPAS1 knockdowns (logFC Ͻ Ϫ0.5, FDR-adjusted p Ͻ 0.05) and not different between single versus combination knockdowns (͉logFC͉ Ͻ 0.5, FDR-adjusted p Ͼ 0.2) was used. Lastly, for the remaining genes to be considered as dependent on either HIF for their expression, they showed reduced expression in the combination knockdown (logFC Ͻ Ϫ0.5, FDR-adjusted p Ͻ 0.05) but not the individual HIF gene knockdowns (logFC Ͼ 0 or FDR-adjusted p Ͼ 0.2) nor between HIFs (FDRadjusted p Ͼ 0.05). Lastly, all of these genes were grouped into categories based on whether their expression increased in hypoxia Gene knockdowns P6 HUVECs were allowed to grow to ϳ85% confluency prior to siRNA transfection. Independent siRNAs targeting HIF1A (A: D-004018-05 and B: D-004018-03), EPAS1 (A: D-004814 -02 and B: D-004814 -03), and nontargeting control (Scr: D-001210-02) were purchased from Dharmacon (Lafayette, CO). DharmaFECT4 (catalog no. T-2004-01) was used with Optimem (ThermoFisher) for siRNA transfections. Effective total siRNA concentrations were determined to be 15 nM (7.5 nM for each HIF1A and EPAS1 in combined knockdown). HIF-1␣ and HIF-2␣ protein levels were measured immediately after 7 h of exposure to 21% oxygen or 1% oxygen. Whole cell lysates were prepared with radioimmune precipitation assay lysis buffer and sample buffer (2ϫ Laemmli ϩ 5% beta-mercaptoethanol), and proteins were detected with Western blotting. Antibodies were obtained from Cell Signaling Technology (Danvers, MA) (HIF-1␣: catalog no. 3716 and HIF-2␣: catalog no. 7096), Abcam (Cambridge, MA) (␤-actin: catalog no. ab8229), and R&D Systems (Minneapolis, MN) (SOX7: catalog no. AF2766). mRNA levels of HIF1A and EPAS1 following knockdowns were confirmed subsequently with RNA-Seq.

Angiogenesis assays
Low passage (P6) HUVECs were grown overnight and siRNA-treated (SOX7 or scramble) the following morning. The cells were washed and detached 24 h after transfection and resuspended in fresh medium containing no VEGF. The cells were seeded into wells of a 96-well plate coated with Matrigel SOX7 regulates hypoxia-dependent angiogenesis (10, 000 cells/well, n ϭ 3). The plates were incubated in hypoxia as described above using 1% oxygen or in 21% oxygen. The center of each well was imaged at 6 h and 24 h using 4ϫ magnification. Branching points, branch lengths, and number of branches were quantitated using ImageJ software. Two independently targeting SOX7 siRNAs were used (SOX7-1: D-019017-01-0002 and SOX7-2: D-019017-02-0002; Dharmacon).

ChIP-PCR
Low passage (P6) HUVECs were grown to confluency followed by 1% formaldehyde treatment (10-min incubation at room temperature). Cross-linking was quenched with 2.5 M glycine, and cells were scraped into a hypotonic buffer on ice. Following centrifugation to pellet nuclei, the pellets were resuspended in radioimmune precipitation assay buffer and protease inhibitor and then sheared. DNA shearing was examined on agarose gel to ensure 200 -600 bp fragments were present. Protein G magnetic beads were used to reduce nonspecific binding, and cleared lysates were incubated with 10 g HIF-2␣ antibody (ab199; Abcam) overnight. After washing and eluting immunoprecipitates, cross-linking was reversed with NaCl, and DNA was purified with a DNA spin column (Qiagen).