Profiling of transcriptional and epigenetic changes during directed endothelial differentiation of human embryonic stem cells identifies FOXA2 as a marker of early mesoderm commitment

Introduction Differentiation of vascular endothelial cells (ECs) in clinically relevant numbers for injection into ischaemic areas could offer therapeutic potential in the treatment of cardiovascular conditions, including myocardial infarction, peripheral vascular disease and stroke. While we and others have demonstrated successful generation of functional endothelial-like cells from human embryonic stem cells (hESCs), little is understood regarding the complex transcriptional and epigenetic changes that occur during differentiation, in particular during early commitment to a mesodermal lineage. Methods We performed the first gene expression microarray study of hESCs undergoing directed differentiation to ECs using a monolayer-based, feeder-free and serum-free protocol. Microarray results were confirmed by quantitative RT-PCR and immunocytochemistry, and chromatin immunoprecipitation (ChIP)-PCR analysis was utilised to determine the bivalent status of differentially expressed genes. Results We identified 22 transcription factors specific to early mesoderm commitment. Among these factors, FOXA2 was observed to be the most significantly differentially expressed at the hESC–EC day 2 timepoint. ChIP-PCR analysis revealed that the FOXA2 transcription start site is bivalently marked with histone modifications for both gene activation (H3K4me3) and repression (H3K27me3) in hESCs, suggesting the transcription factor may be a key regulator of hESC differentiation. Conclusion This enhanced knowledge of the lineage commitment process will help improve the design of directed differentiation protocols, increasing the yield of endothelial-like cells for regenerative medicine therapies in cardiovascular disease.


Introduction
The directed differentiation of human embryonic stem cells (hESCs) towards endothelial cell (EC) lineages offers therapeutic potential in the treatment of myocardial infarction, peripheral vascular disease and stroke [1]. While successful derivation of endothelial-like cells from hESCs has been demonstrated [2][3][4][5], relatively little is understood regarding early commitment to mesoderm and subsequent specification [6].
To commit to a specified lineage, pluripotent cells must undergo radical transcriptional change [2], partly regulated by miRNAs [2,3]. Recent studies suggest epigenetic influences also play a significant role in the determination of cell fate [7]. Indeed, the poised transcriptional state classically associated with pluripotent hESCs is maintained via effects at the chromatin level with gene expression determined by post-translational modification of histones [8]. While some histone marks, such as trimethylation of lysine 4 of histone H3 (H3K4me3), are associated with gene activation, others, such as trimethylation of lysine 27 of histone H3 (H3K27me3), are associated with repression. However, around 3,000 developmental regulatory genes in embryonic stem cells are labelled with both H3K4me3 and H3K27me3 [7], allowing the genes to be rapidly activated upon differentiation or to remain silenced during commitment to lineages not requiring their expression. Unsurprisingly, these bivalently marked genes are proposed to be master regulators of the differentiation process, although their role in endothelial differentiation has not been extensively investigated to date.
The study described herein was designed to profile transcriptional and epigenetic changes during early hESC commitment to a mesodermal and endothelial-like fate with a view to improving understanding of this process and to optimise the generation of ECs for regenerative medicine purposes.

Methods
Cell culture hESC lines SA461 (Cellartis, Dundee, UK), H1 and H9 (WiCell Research Institute, Madison, WI, USA) and RC10 (Roslin Cells Ltd, Edinburgh, UK) were cultured in a monolayer-based, serum-free and feeder-free system. Pluripotency was maintained and endothelial differentiation induced as previously described [2]. Primary human saphenous vein endothelial cells (HSVECs) were isolated on the day of surgery by standard collagenase digestion based on a modified version of the protocol described by Jaffe and colleagues [9]. HSVECs were then cultured as previously described [2]. All participants gave written informed consent. The study was approved by the West of Scotland Ethics Committee (06/S0703/110) and complies with the principles of the Declaration of Helsinki.

Microarray analysis
SA461 hESCs subjected to directed endothelial differentiation were harvested at days 0, 2, 4 and 10, and RNA was isolated using the miRNAeasy Mini Kit (Qiagen, Crawley, UK). HSVEC RNA was harvested to act as a positive control. Complementary RNA was prepared for hybridisation with Illumina® HT-12 v3 Expression BeadChip microarrays (Illumina®, Saffron Walden, UK) and the gene expression assessed. Microarray data were deposited in the ArrayExpress public repository [ArrayExpress: E-MTAB-1510]. Figure 1A outlines the experimental design.
Data were quantile normalised and background subtracted in Genome Studio (Illumina®) and were then exported to Partek® Genomics Suite™ (Partek® Inc., St. Louis, MO, USA) and visualised using principle component analysis [10]. An error-weighted analysis of variance with a false discovery rate multiple testing adjustment threshold of 0.05 was used to identify differentially expressed probe sets [11], which were then uploaded to Ingenuity Pathway Analysis software (2009; Ingenuity® Systems, Redwood City, CA, USA) and analysed to identify dynamic gene expression.

In silico prediction of bivalency
H3K27me3 and H3K4me3 data from the H9 chromatin immunoprecipitation (ChIP)-sequencing dataset of Ku and colleagues [12] were mined and integrated with microarray data to give an in silico prediction of the bivalent status of genes.
SICER software (http://home.gwu.edu/~wpeng/Software. htm) was used to determine enriched domains for each histone modification using a stringent E value of 0.1, a window size of 200 and gap sizes of 200, 400, 600, 800 and 1,000 base pairs [13]. For each gene within the Ensembl gene set (NCBI36.1), the transcription start site (TSS) was determined and an interval spanning ±2 kb around the TSS was defined.
For a specific gap size, if the TSS contained enriched domains for both H3K27me3 and H3K4me3 marks, then this gene was classified as bivalent. The bivalency score was defined as the count of gap sizes where a gene was considered bivalent. A score of 0 would thus indicate that the gene was never classified as bivalent, while a score of 5 indicated that the gene was found to be bivalent regardless of the gap size chosen.

Chromatin immunoprecipitation and PCR identification
ChIP assays were performed based on a modified version of the method of Rai and colleagues [15]. Chromatin was prepared from pluripotent H9s and SA461s, and H3K4me3 and H3K27me3 were immunoprecipitated using Dynabeads® M-280 sheep anti-rabbit IgG (Invitrogen) and H3K4me3 (AB8580; Abcam, Cambridge, UK) and H3K27me3 (C36B11; Cell Signaling Technology, Beverly, MA, USA) specific antibodies. Immunoprecipitations with total H3 (AB1791; Abcam) and control IgG (M7023; Sigma-Aldrich Company Ltd., Dorset, UK) were included as positive and negative controls, respectively.
Using the University of California Santa Cruz Genome Browser, primer pairs were designed to span the FOXA2 TSS (Additional files 1, 2 and 3) and DyNAmo™ SYBR® Green quantitative PCR (Thermo Fisher Scientific UK Ltd., Loughborough, UK) was performed on immunoprecipitation eluates, in addition to 2% chromatin input not subjected to immunoprecipitation. Quantitative PCR data were normalised to IgG negative control and displayed as fold enrichment.

Statistical analyses
Values are presented as mean ± standard error of the mean. Data from multiple groups were analysed using repeatedmeasures analysis of variance. Significant differences were determined by Tukey post-hoc testing and P <0.05 (twotailed) was considered significant.

Gene expression analysis of endothelial differentiation
Principle component analysis of global transcription data, designed to detect early transcriptional changes during directed differentiation ( Figure 1A), revealed sufficient separation of cell groups. As expected, HSVECs were clearly distinct from all hESC-derived cells, and hESC-EC on day 4 and day 10 were divergent compared with day 0 and day 2 timepoints (Figure 1B). A large number of significantly differentially expressed probe sets was observed at each of the three differentiation timepoints, as compared with day 0 ( Figure 1C). Comparison of samples from day 10 hESC-EC with HSVEC samples revealed differential expression of 6,133 different probe-sets, defining markedly different cells ( Figure 1C).
Further investigation of the hESC-EC day 2 timepoint was carried out and 223 significantly differentially expressed, unique probe sets were identified ( Figure 1C). We reasoned that corresponding genes might be involved in early differentiation. The Ingenuity Pathway Analysis software revealed that these probe sets correspond to 178 genes (Additional file 4), a relatively high proportion (>12%) of which encode for transcription factors (Table 1).
To validate this differential expression of transcription factors, quantitative RT-PCR was performed. Of the 22 genes listed, 12 were confirmed as being differentially expressedincluding FOXA2, which demonstrated the most significant and rapid upregulation of all genes examined (Figure 2A). Upregulation of FOXA2 at hESC-EC day 2 was also confirmed in additional H1 and RC10 hESC lines ( Figure 2B). Again, upregulation was transient as analysis of day 4, day 7 and day 10 differentiations failed to detect expression of FOXA2 mRNA (data not shown).
To confirm day 2 differential expression of FOXA2 at the protein level, SA461s, H1s and RC10s were stained using anti-FOXA2 antibodies. Expression of the pluripotency marker OCT4 was observed in day 0 pluripotent cells but diminished by hESC-EC day 2 ( Figure 2C). FOXA2, although present in a subset of pluripotent cells, was abundantly expressed at day 2. Interestingly, FOXA2 did not co-localise with OCT4 ( Figure 2C).

Epigenetic analysis of FOXA2
On aligning in silico data with differentially expressed genes identified via microarray analysis, a large proportion (2,684/3,883) achieved a bivalency score of 5 and was therefore predicted to be bivalently marked.
With epigenetics known to influence cell fate, the predicted bivalent status of genes encoding hESC-EC day 2 differentially expressed transcription factors was of particular interest. FOXA2 was identified as having both H3K4me3 and H3K27me3 modifications at its TSS (Table 1) and ChIP-PCR confirmed bivalency of the FOXA2 TSS in pluripotent H9s and SA461s (Figure 3).

Discussion
To the best of our knowledge, this study describes the first global gene expression profiling of hESCs at intermediate timepoints during their directed differentiation to an endothelial-like lineage using a monolayer-based, feeder-free and serum-free protocol. The study is also the first to identify transcription factors potentially specific to early mesoderm commitment, including FOXA2 which we show to be bivalently marked at the TSS. In addition, we report that cells at day 10 of the differentiation process remain genetically distinct from mature endothelial cells, supporting our previous findings that this timepoint yields progenitor-type cells, capable of further differentiation in vivo [2].
Rapid and transient FOXA2 mRNA upregulation was identified and confirmed at the protein level, suggesting the hepatocyte nuclear factor may represent a marker of early mesoderm lineage commitment. A member of the forkhead class of DNA binding proteins, FOXA2 has been shown to bind and open compacted chromatin at histones H3 and H4 [16] and is one of the first transcription factors to bind target genes on differentiation [17]. Traditionally endoderm associated, FOXA2 has not, as far as we are aware, previously been linked to hESC-EC differentiation. However, with cells used for microarray analyses likely to represent a somewhat heterogeneous population containing precursors for both endoderm and mesoderm germ layers [18,19], further investigation utilising purified cell populations will be required to confirm the importance of FOXA2 in early EC commitment.
Epigenetic control of transcription is now known to play a major role in cell fate determination and, as stated, we have identified the FOXA2 TSS as being bivalently marked in pluripotent hESCs, displaying histone modifications associated with both gene activation and repression. Conceivably, FOXA2 may therefore be a marker of early mesoderm lineage commitment and a potential regulator of hESC-EC commitment. Such enhanced knowledge of the complex commitment process is likely to prove important for optimising the design and development of directed endothelial differentiation protocols.

Conclusion
This study is the first to generate global gene expression profiles for hESCs undergoing directed differentiation to ECs using a monolayer-based, serum-free and feederfree system, and is the first to identify transcription factors, including FOXA2, likely to be specific for early mesoderm lineage commitment.

Additional files
Additional file 1: A figure showing the University of California Santa Cruz Genome Browser visualisation of H3K4me3 and H3K27me3 ChIP sequencing performed on pluripotent H9 hESCs. Genome browser output from FOXA2 genomic location showing Figure 3 In vitro validation of predicted FOXA2 bivalency. Immunoprecipitation of H3K4me3, H3K27me3 and total H3 from chromatin of pluripotent SA461s and H9s revealed enrichment of the histone modifications at the FOXA2 transcriptional start site, as identified by quantitative PCR using primers for three regions, covering around 2.5 kb in total. *P <0.05 versus IgG control.
presence of H3K4me3 and H3K27me3 at the transcriptional start site. Output is based on the pluripotent H9 hESC ChIP sequencing data of Ku and colleagues [12].