ETV2 Upregulation Marks the Specification of Early Cardiomyocytes and Endothelial Cells During Co-differentiation

Abstract The ability to differentiate human-induced pluripotent stem cells (hiPSCs) efficiently into defined cardiac lineages, such as cardiomyocytes and cardiac endothelial cells, is crucial to study human heart development and model cardiovascular diseases in vitro. The mechanisms underlying the specification of these cell types during human development are not well understood which limits fine-tuning and broader application of cardiac model systems. Here, we used the expression of ETV2, a master regulator of hematoendothelial specification in mice, to identify functionally distinct subpopulations during the co-differentiation of endothelial cells and cardiomyocytes from hiPSCs. Targeted analysis of single-cell RNA-sequencing data revealed differential ETV2 dynamics in the 2 lineages. A newly created fluorescent reporter line allowed us to identify early lineage-predisposed states and show that a transient ETV2-high-state initiates the specification of endothelial cells. We further demonstrated, unexpectedly, that functional cardiomyocytes can originate from progenitors expressing ETV2 at a low level. Our study thus sheds light on the in vitro differentiation dynamics of 2 important cardiac lineages.


Introduction
In vivo, cardiomyocytes (CMs) and endothelial cells (ECs) originate from Mesp1+ progenitors specified during gastrulation. In mice, these cells appear in the primitive streak and subsequently migrate toward the lateral plate mesoderm around E6.5.- [1][2][3][4] The timing of segregation of CMs and ECs from their common progenitor is still controversial. Singlecell RNA-seq (scRNA-seq) of mouse Mesp1+ progenitors collected at E6.75 and E7. 25 showed that these cells were already segregated into distinct cardiovascular lineages, including CMs and ECs. 5 However, other studies showed that multipotential progenitors were still present in Flk-1expressing lateral plate mesoderm. 6,7 These cells were the first to be recognized as multipotent cardiac progenitor cells (CPCs). 8 Studies in mouse and chick showed that CPCs come from 2 different sources 9,10 : the first-and the second-heart field (FHF, SHF). The FHF in the cardiac crescent contributes to the primitive heart tube, which serves as a scaffold into which SHF cells can migrate before heart chamber morphogenesis. It has been shown that cells from the SHF are patterned before migration to give rise to different parts of the heart. 3,11 CPCs from FHF and SHF can be distinguished by the expression of ISL1, which is specific to the SHF. 12 Nkx2-5-expressing CPCs in both FHF and SHF from E7.5 to E7.75 contribute to both CMs and ECs in the heart. 13 ETS Variant Transcription Factor 2 (Etv2) is a master regulator of endothelial and hematopoietic cell lineages during early development. 14 Etv2 functions downstream of BMP, WNT, and NOTCH signaling pathways 15 and regulates the expression of early EC-specific markers, such as Tal1, Gata2, Lmo2, Tek, Notch1, Notch4, and Cdh5.- [15][16][17][18] In mouse embryonic stem cells (ESCs), VEGF-FLK1 signaling upregulates ETV2 expression to induce hemangiogenic specification via an ETV2 threshold-dependent mechanism. 19 ETV2 expression was also found to direct the segregation of hemangioblasts and smooth muscle cells (SMCs) in mouse ESCs. 20 In human heart development, much less is known about the specification of endothelial and myocardial lineages from multipotent CPCs, both in terms of timing and gene regulatory mechanisms. More specifically, it is still unclear whether ETV2 also plays a role in the segregation of ECs and CMs from CPCs in humans. Overexpression of ETV2 converts human fibroblasts into endothelial-like cells 21 and ETV2 expression levels have been modified in several studies to drive hiPSCs toward ECs in 2D and 3D cultures.- [22][23][24][25][26][27][28] Paik et al performed scRNA-seq analysis of hiPSC-derived ECs (hiPSC-ECs), which made up less than 10% of the cells that expressed the cardiac maker TNNT2. The developmental dynamics of ECs and cardiac lineages as such were not further studied. 29 In an scRNA-seq study of hiPSC-ECs obtained using a different differentiation protocol, 30 ECs were collected at multiple time points. This study showed that endothelial and mesenchymal lineages have a common developmental origin in mesoderm cells but the identity and differentiation potential of these cells were not described.
Previously, we found that MESP1+ progenitors derived from human ESCs could give rise to CMs, ECs, and SMCs. 31,32 We also developed a co-differentiation system for ECs and CMs from hiPSCs through a common cardiac mesoderm precursor. 33 Here we performed scRNA-seq analysis of this co-differentiation system to elucidate the relationship between ETV2 expression and specification of ECs and CMs from cardiac mesoderm. ETV2 expression was observed principally as an initial "pulse" in the endothelial lineage but also in a subpopulation of the myocardial lineage. Using a newly generated ETV2 mCherry hiPSC reporter line, we purified 2 subpopulations of ETV2+ cells and revealed their derivative EC and CM expression characteristics by bulk RNA-seq. These sorted populations also showed distinct differentiation potentials toward CMs and ECs upon further differentiation with VEGF. In summary, this study detailed ETV2 dynamics during the segregation of human CMs and ECs differentiated from hiPSCs.

hiPSC Culture
The NCRM1 hiPSC line (NIH Center for Regenerative Medicine (NIH CRM), obtained from RUDCR Infinite Biologics at Rutgers University, hPSCreg number CRMi003-A) was used in this study, except for the single-cell RNA-seq, which was done with LUMC0020iCTRL06 (hPSCreg number LUMCi028-A). hiPSC control lines were cultured in TeSR-E8 on Vitronectin XF and routinely passaged once a week using Gentle Cell Dissociation Reagent (all from STEMCELL TECHNOLOGIES). Prior to targeting, NCRM1 hiPSCs were passaged as a bulk on feeders in hESC medium. 34 RevitaCell (Life Technologies) was added to the medium (1:200) after every passage to enhance viability after single-cell passaging with TrypLE (Life technologies).

Generation of hiPSC Reporter Line Using CRISPR/ Cas9
The p15a-cm-hETV2-P2A-NLS-mCherry-neo repair template plasmid was generated using overlap PCR and restriction-based cloning and ligation. The ETV2 homology arms were amplified from genomic DNA and the neomycin cassette flanked by 2 flippase recognition target (FRT) sites were amplified from the P15 backbone vector (kindly provided by Dr Konstantinos Anastassiadis, Technical University Dresden). P2A-NLS-mCherry double-stranded DNA fragment was ordered from IDT. The sgRNA/Cas9 plasmid was modified from SpCas9-2A-Puro V2.0 plasmid (Addgene, Feng Zang). NCRM1 hiPSCs were passaged at a 1:2 or 1:3 ratio into 60 mm dishes to reach 60%-70% confluence the next day for transfection. Twenty microliters of lipofectamine (Invitrogen), 8 µg of a repair template, and 8 µg of sgRNA/Cas9 plasmid were diluted in 600 µL of Opti-MEM and added to each 60 mm dish. After 18 h the medium was changed to hESC medium. After another 6 h, G-418 (50 µg/mL) selection was started and was continued for 1 week. Surviving cells were cultured in hESC medium and passage into 6-well plates for the transfection of Flp recombinase expression vector to remove the neomycin cassette. 35 A total of 300 µL of Opti-MEM containing 10 µL lipofectamine and 4 µg CAGGs-Flpo-IRESpuro plasmid was added per well for 18 h. Puromycin (0.5 µg/ mL) selection was started 24 h post-transfection for 2 days. Once recovered, cells were passaged into 96-well plates for clonal expansion via limiting dilution. Targeted clones were identified by PCR and sequencing. Primers outside the ETV2 homology arms and primers inside the targeting construct were used to confirm on-target integration. The absence of mutations within the inserted sequence and untargeted allele was confirmed by Sanger sequencing (BaseClear).

Immunofluorescence Staining and Imaging
Cultured cells were fixed in 4% paraformaldehyde for 15 minutes, permeabilized for 10 minutes with PBS containing 0.1% Triton X-100 (Sigma-Aldrich), and blocked for 1 h with PBS containing 5% BSA (Sigma-Aldrich). Then cells were stained with primary antibody overnight at 4 °C. The next day, cells were washed 3 times (20 minutes each time) with PBS. After that, cells were incubated with fluorochromeconjugated secondary antibodies for 1 h at room temperature and washed 3 times (20 minutes each time) with PBS. Then, cells were stained with DAPI (Life Technologies) for 10 minutes at room temperature and washed once with PBS for 10 minutes. Both primary and secondary antibodies were diluted in 5% BSA/PBS. Images were taken with the EVOS FL AUTO2 imaging system (Thermo Fischer Scientific) with a 20× objective or using the Incucyte system (Sartorius). Confocal imaging was done using a Leica SP8WLL confocal laser-scanning microscope using a 63× objective and z-stack acquisition. Details of all antibodies used are provided in Supplementary Table S1.

FACS Analysis
Cells were washed once with FACS buffer (PBS containing 0.5% BSA and 2 mM EDTA) and stained with FACS antibodies for 30 minutes at 4 °C. Samples were washed once with FACS buffer and analyzed on the MACSQuant VYB (Miltenyi Biotech) equipped with a violet (405 nm), blue (488 nm), and yellow (561 nm) laser. The results were analyzed using Flowjo v10 (FlowJo, LLC). Details of all fluorochromeconjugated FACS antibodies are provided in Supplementary  Table S1.

Quantitative Real-Time Polymerase Chain Reaction (qPCR)
Total RNA was extracted using the NucleoSpin RNA kit (Macherey-Nagel) according to the manufacturer's protocol. cDNA was synthesized using an iScript-cDNA Synthesis kit (Bio-Rad). iTaq Universal SYBR Green Supermixes (Bio-Rad) and the Bio-Rad CFX384 real-time system were used for the PCR reaction and detection. Relative gene expression was determined according to the standard ΔCT calculation and normalized to housekeeping genes (mean of HARP and RPL37A). Details of all primers used are provided in Supplementary Table S2.
Differentially expressed genes were identified using generalized linear models as implemented in edgeR (3.24.3). 42 P-values were adjusted using the Benjamini-Hochberg procedure and FDR ≤.05 was considered significant. Analyses were performed using R (version 3.5.2). The PCA plot was generated with the built-in R function prcomp using the transposed normalized RPKM matrix. Correlation among samples was calculated using the cor function with the Spearman method and the correlation heatmap was generated with a heatmap function (NMF package).
Gene clusters were calculated with the CancerSubtypes package. 43 The top 3000 most variable genes across all chosen samples were identified based on the Median Absolute Deviation (MAD) using the FSbyMAD function, then expression was standardized for each gene. K clusters were calculated using k-means clustering with Euclidean distance. Clustering was iterated 1000 times for each k in the range of 2-10. Heatmaps of genes in all clusters were generated using the base R heatmap function. Gene ontology enrichment for each cluster was performed using the compareCluster function of clusterProfiler package (v3.10.1) 44 and q ≤ .05 was considered significant.

Single-Cell RNA Sequencing and Analysis Sample Preparation and Sequencing
Cells were dissociated into single cells on day 6 of CMEC differentiation and loaded into the 10× Chromium Controller for library construction using the Single-Cell 3ʹ Library Kit. Next, indexed cDNA libraries were sequenced on the HiSeq4000 platform. Mean reads per cell were 28 499 in the first replicate and 29 388 in the second replicate.

Preprocessing
Both replicates of day 6 CMEC differentiation were merged into 1 dataset. The average number of detected genes was 2643 and the average total expression per cell was 10 382 ( Supplementary Fig. S1A, S1B). Then, undetected genes (>1 UMI count detected in less than 2 cells) and cells with a low number of transcripts were removed from further analysis ( Supplementary Fig. S1A, S1B). This resulted in 5107 cells for the first replicate, and 3743 cells for the second replicate and 13 243 genes. Expression profiles were normalized with the R package scran (V 1.10.2) using the method described in Ref. 45 The 5% most highly variable genes (HVGs) for each replicate were calculated with scran after excluding ribosomal genes (obtained from the HGNC website without any filtering for minimum gene expression), stress-related genes 46 , and mitochondrial genes. For downstream analysis, the top 5% HVGs were used after excluding proliferation 47 and cell cycle 48 -related genes.

Cell-Cycle Analysis
For each dataset, cell-cycle analysis was performed with the scran package using the cyclone function 49 on normalized counts. Cells with a G2/M score higher than 0.2 were considered to be in G2/M phase. Otherwise, they were classified as G1/S. Using this binary classifier as predictor, we regressed out cell cycle effects with the R package limma (V 3.42.2) 50 applied to log-transformed normalized counts. The 2 replicates were then batch corrected with the fast mutual nearest neighbors (MNN) correction method 51 on the cell cycle corrected counts, using the 30 first principal components and 20 nearest neighbors.

Clustering
Batch-corrected counts were standardized per gene and then used to create a shared nearest neighbor (SNN) graph with the scran R package (d = 30, k = 20). Louvain clustering was applied to the SNN graph using the igraph python package (V 0.7.1) with 0.4 as the resolution parameter. This resulted in 5 clusters. Two of these 5 clusters were excluded from further analysis based on the expression of pluripotency markers. 50

Dimensionality Reduction and Pseudotime
Dimensionality reduction was performed using the python scanpy pipeline (V 1.4.6). First, a 20 nearest-neighbors (knn, k = 20) graph was created from diffusion components of the batch corrected datasets. Diffusion components are the eigenvectors of the diffusion operator which is calculated from Euclidean distances and a Gaussian kernel. The aim is to find a lower dimensional embedding that considers cellular dynamics. The graph was projected into 2 dimensions with the default force-directed graph layout and starting positions obtained from the partition-based graph abstraction (PAGA) algorithm. 52 PAGA estimates connectivities between partitions and performs an improved version of diffusion pseudotime. Diffusion pseudotime 51,52 was calculated on these graphs with root cells selected from the "Cardiac Mesoderm" cluster.
Average gene expression trajectories were calculated by dividing the cells of each cluster into bins along pseudotime. Fifty bins were created for cardiac mesoderm and 30 bins each for ECs and CMs. The average log-expression per bin was then calculated. The value of the threshold indicated in Fig. 1D, 1E was determined by calculating the point in pseudotime where the average ETV2 expression was the lowest in the endothelial cell cluster before the peak expression, which corresponds to a value around 0.25.

Differential Expression Analysis and Identification of Cluster Maker Genes
The R package edgeR (V 3.24.3, 31) 42 was used to perform differential expression analysis. We used raw counts and a negative binomial distribution to fit the generalized linear model. The covariates were comprised of 6 binary dummy variables that indicate the 3 remaining clusters per replicate and a variable that corresponds to the total number of counts per cell. Finally, P-values for each cluster considering both replicates were obtained and adjusted for multiple hypothesis testing with the Benjamini-Hochberg method.

Comparison to Bulk RNA-Sequencing Data
The MNN approach was used to integrate the 2 single-cell replicates, using normalized counts and the 10% HVGs per replicate, and the bulk RNA-sequencing data, with d = 30 and k = 20. After batch correction, a diffusion map was calculated on the MNN-corrected values with default parameters.

Statistics
Statistical analysis was conducted with GraphPad Prism 7 software. Data are represented as mean ± SD. A Student's t test was used for the comparison of the 2 samples. Ordinary one-way ANOVA was used for multiple sample comparison, and uncorrected Fisher's LSD test was applied. Two-way ANOVA was used for multiple group comparison and uncorrected Fisher's LSD test was applied. P < .05 was considered significant.

ETV2 is Upregulated after Bifurcation of Progenitors into CMs and ECs
To characterize the expression of ETV2 during codifferentiation of ECs and CMs 33 (Fig. 1A), we collected scRNA-seq data on day 6 of differentiation (Fig. 1B). We identified 3 distinct clusters: cardiac mesoderm, CMs, and ECs (Fig. 1B, top right and Supplementary Table S3). Pseudotime analysis revealed cardiac mesoderm as the common developmental origin of CMs and ECs (Fig. 1B, bottom right). We found that ETV2 was highly expressed in the EC cluster, as well as in a small fraction of cells in the cardiac mesoderm and CM clusters (Fig. 1B, left). We next focused on ETV2 expression dynamics along the developmental path from cardiac mesoderm to CMs and ECs. Notably, ECs extended to larger pseudotimes (0.15-0.8) compared to CMs (0.15-0.3), which might indicate faster differentiation kinetics in the EC lineage (Fig. 1C, Supplementary Fig. S1A). After the bifurcation into ECs and CMs (around pseudotime 0.15), ETV2 increased only slightly in the CM lineage. In the EC lineage, however, it was initially strongly upregulated (until pseudotime 0.25) and subsequently declined to a similar level as in cardiac mesoderm (Fig. 1C). ETV2 downstream target genes, such as TAL1, GATA2, and LMO2, 18 were only slightly increased in the CM lineage (Fig. 1D), while in the EC lineage, they were highly induced and strongly correlated with ETV2 (Fig. 1E). Notably, TAL1, GATA2, and LMO2 only showed significant expression after ETV2 expression exceeded 0.25 in ECs, an expression level that was not reached in CMs (Supplementary Fig.  S1B, S1C). Endothelial-specific genes KDR, CD34, SOX17, CDH5, and PECAM1 increased on the path from cardiac mesoderm to ECs (Supplementary Fig. S1E). Most of these genes started to increase when ETV2 was already declining, as exemplified by CDH5 (Supplementary Fig. S1D). Genes related to cardiac or muscle function, like ACTC1, PDLIM5, HAND1, PKP2, and GATA4, most of which were already expressed in the cardiac mesoderm, were further increased in the CM lineage ( Supplementary Fig. S1F). Identification of genes that are differentially expressed between ETV2+ CMs and ECs showed enrichment in CM-and EC-specific genes, respectively (Fig. 1F, Supplementary Table S4). Taken together, these analyses confirmed the differentiation of cardiac mesoderm into CMs and ECs, which we had discovered previously. They also revealed the increase of ETV2 as a global genes that were differentially expressed between ETV2+ ECs and ETV2+ CMs in the scRNA-seq dataset. 128 and 136 genes were upregulated in ETV2+ ECs and CMs respectively (P adjusted < .05). A complete list of GO terms can be found in Supplementary Table S3. Color represents the P adjusted of the enrichment analysis and dot size represents the count of genes mapped to the GO term.
indicator of early lineage separation and a transient pulse of high ETV2 at the beginning of EC specification.

Generation and Characterization of an ETV2 mCherry Fluorescent hiPSC Reporter Line
To follow ETV2 expression in real-time, we introduced a fluorescent reporter for ETV2 in hiPSCs. P2A-mCherry with a nuclear localization signal (NLS) and a neomycin selection cassette was inserted into the endogenous ETV2 locus before the stop codon using CRISPR/Cas9-facilitated homologous recombination ( Fig. 2A; Supplementary Fig. S2A). After neomycin selection and excision of the selection cassette, targeted hiPSC clones were validated by PCR ( Supplementary   Fig. S2A-S2D) and Sanger sequencing (data not shown). The hiPSC clone with ETV2 mCherry in both alleles was further characterized by measuring pluripotency marker expression and G-banding karyotyping ( Supplementary Fig. S2E-S2H). Karyotype analysis revealed an additional duplication in the 1q32.1 locus (Supplementary Fig. S2H). This duplication occurs frequently in hPSCs possibly imposing positive natural selection. 53 However, this did not appear to affect the differentiation of the hPSCs to ECs.
ETV2 and mCherry mRNAs were highly expressed on days 4-5 of differentiation and downregulated from day 6 (Fig. 2B,  2C). ETV2 and mCherry protein appeared from day 4 and peaked on day 5. ETV2 protein was downregulated on day 6 and absent on day 8. mCherry was retained in a fraction of the cell population for somewhat longer because of its halflife (Fig. 2D, 2E; Supplementary Fig. S3A and Supplementary Online Video 1). Flow cytometry analysis at different stages of differentiation revealed upregulation of ETV2 (mCherry protein) as early as day 4 of differentiation followed by the upregulation of the EC-specific marker CD144 ( Fig. 2F;  Supplementary Fig. S3B). Quantification of the percentages of single positive (SP; ETV2 mCherry +CD144-) and doublepositive (DP; ETV2 mCherry +CD144+) cells on days 4, 5, 6, and 8 of differentiation from at least 3 independent experiments showed a decrease and an increase of SP and DP cells, respectively (Fig. 2G, 2H). mCherry protein remained present for a longer period than ETV2 protein and endogenous ETV2 and mCherry mRNA (Fig. 2B-2H), likely due to the relatively longer half-life of the mCherry protein. This explains the persistence of mCherry signal in both the DP and SP population (Fig. 2G, 2H), and offers the possibility to use mCherry as a lineage tracer, identifying cells that previously passed through a stage of being ETV2+.

The ETV2 mCherry Fluorescent Reporter Allows the Purification of Differentiating Cells with Lineage-Specific Expression Profiles
We next sorted DP and SP cells at different stages of differentiation (Fig. 2F) and performed bulk RNA-seq on at least 3 independent replicates. ETV2 mRNA showed similar trends in DP and SP cells ( Supplementary Fig. S4), consistent with the earlier qPCR result (Fig. 2B).
Principal component analysis (PCA) showed that DP and SP populations diverged progressively over time (Fig. 3A). Mapping of the bulk transcriptomes to the scRNA-seq data revealed that DP samples aligned on the EC branch and SP cells on the CM branch (Fig. 3B). Notably, SP and DP cells collected at later time points were further away from cardiac mesoderm, reflecting ongoing differentiation (Fig. 3B).
Taken together, time-resolved bulk RNA-seq of sorted SP and DP populations confirmed that ETV2-positive cells contained transcriptionally distinct subpopulations. DP cells were part of the EC lineage, while SP cells corresponded to the CM lineage.

ETV2+ Cells Contain Lineage-Predisposed Subpopulations
Next, we wanted to find out how the various subpopulations we identified differed in terms of their further differentiation potential. To this end, we sorted cells on the basis of ETV2 reporter levels shortly after the bifurcation (on day 5) and attempted to differentiate them further toward ECs by adding VEGF (Fig. 4A, 4B). After 5 days of additional differentiation, ETV2+ cells produced more than 90% CD144+CD31+ ECs, while ETV2-cells gave rise to only 10%-15% ECs (Fig. 4C,  4D; Supplementary Fig. S6A). Only cells derived from ETV2+ cells expressed endothelial-specific markers, as observed by qPCR and immunofluorescence (Fig. 4E-4H; Supplementary  Fig. S6C). These cells also upregulated pro-inflammatory markers, such as ICAM-1 and E-Selectin upon TNF-α stimulation ( Fig. 4I-4L; Supplementary Fig. S6B), as shown previously for hiPSC-derived ECs. 54 We thus concluded that the majority of ETV2+ cells on day 5 has a strong propensity to produce ECs.
Both the analysis of the scRNA-seq data and the time-resolved bulk RNA-seq of sorted cells identified a subpopulation of ETV2+ cells with CM characteristics. We strongly suspected that the differentiation of these cells would  be predisposed to the CM lineage. To test this hypothesis with our reporter line, we co-differentiated cells until day 7. We chose a later time point for this experiment because the majority of cells are past bifurcation at this point and it is therefore easier to identify the ETV2+ population that does not correspond to early ECs. We co-stained for CD144 and sorted the cells into DP, SP, and double-negative (DN) populations. These subpopulations were then further cultured in the presence of VEGF until day 18 (Fig. 4A, 4M). The majority (>80%) of DP cells differentiated into CD144+CD31+ ECs, in agreement with the previous experiment (Fig. 4N, 4O). In contrast, more than 50% of SP and DN cells differentiated into cTnT+ CMs while very few ECs were detected (Fig. 4N, 4O). Interestingly, CMs derived from SP cells seemed to proliferate more and formed a monolayer composed of a contracting cell sheet, while CMs from DN cells proliferated to a lesser extent and produced only a few, isolated clusters of contracting cells (Supplementary Online Video 2). Almost all DP cells on day 18 expressed the EC marker CD31, while only few cells derived from SP and DN cells were positive for CD31 ( Fig. 4P-4R). Most cells derived from SP and DN expressed CM-specific α-actinin and cTnT and showed typical sarcomeric structures (Fig. 4Q, 4R; Supplementary S6D, S6E). A small number of SP and DN-derived cells were also positive for the smooth muscle cell marker SM22, while negative for cardiac markers (data not shown). Furthermore, the α-actinin-positive CMs derived from the SP cell fraction were positive for SM22, possibly indicating their immaturity (Fig. 4R).
Taken together, the VEGF differentiation experiments showed that DP and SP cells are predisposed to the EC and CM lineages, respectively. DN cells were largely unable to give rise to EC but produced CMs, albeit with lower efficiency than SP cells. Entering a transient state characterized by high ETV2 expression, thus seems necessary to initiate EC specification.

Discussion
In this study, we characterized the dynamics of EC and CM co-differentiation from hiPSCs. 33 ETV2 was identified as an early indicator of lineage segregation and found to be strongly, but transiently, upregulated in ECs, in agreement with its essential role in hemangiogenic development. 55 Interestingly, ETV2 expression was also observed in a small population of cardiac mesoderm and CMs. This is reminiscent of a recent study where etv2 expression was observed in lateral plate mesoderm and the CM population in zebrafish. 56 In our experiments, expression of ETV2 target genes seemed to occur only above a threshold of ETV2 expression, although this observation could also be explained by a temporal delay between ETV2 upregulation and target gene expression. An ETV2 threshold in hiPSC differentiation would be in line with previous reports of an ETV2 threshold in hemangiogenic specification. 19,20 Our results thus support an ETV2 pulseand threshold-dependent specification of ECs.
With the ETV2 mCherry hiPSC reporter line, generated to track, isolate, and characterize ETV2+ cells, we showed that ETV2+ cells could give rise to both ECs and CMs. Over time, EC and CM precursors acquired more specific endothelial and myocardial identities, respectively, as well as downregulating cell cycle-related genes, which indicated exit from the progenitor state and further maturation.
In the DP subpopulation (EC precursors), several key angiogenesis and Notch signaling pathway genes, like LEPR, FOXO4, DLL4, NOTCH4, and EGF, strongly increased starting from day 4, indicating a specified EC fate but an immature state on day 4. These relatively late-expressed genes could potentially be used as markers to distinguish early and late ECs during development in vitro or in vivo. Genes involved in heart development and definitive hematopoiesis were also upregulated during EC development, suggesting a mixture of cardiac endothelial-and hemogenic endothelial identity of these ECs. A better characterization hematopoietic potential of these cells would be interesting but beyond the scope of this study. ECs that were further differentiated with VEGF showed a clear endothelial identity and were fully functional based on their inflammatory response upon TNFα stimulation. Notably, they also expressed a number of cardiac markers like MEOX2, GATA4, GATA6, and ISL1, suggesting a cardiac-specific EC identity. 33 The SP subpopulation (CM precursors) had already committed to a cardiac fate on day 4, as evidenced by the expression of cardiac genes HAND1, MYH10, NKX2-5, ISL1, TNNC1, MYOCD, and LMO4. However, some crucial CM genes were still absent, including MYH6 and TNNT2. MYH6 encodes the major CM thick filament protein MHC-α and TNNT2 is routinely used as a CM marker. Both genes are essential for CM contractility and started to be expressed only after day 4. Their relatively late expression could allow us to identify early and late cardiac progenitors during cardiac development in future studies. CMs were still early progenitors on day 6 of the differentiation as no functionally contracting CMs were observed yet at this stage. Pseudotime analysis also suggested that ECs had differentiated further compared to CMs on day 6. After additional VEGF differentiation, SP cells gave rise to contracting CM, which provided direct evidence they were CM precursors. More importantly, it demonstrated that both ECs and CMs could be derived from ETV2+ progenitors, confirming the presence of a common precursor implied by our earlier studies. 33 Notably, ETV2-cells (DN population) also gave rise to contracting CMs after VEGF treatment, albeit less frequently than SP cells. This difference could be due to either the different cell growth rates or their different developmental origins (FHF vs. SHF). More work is needed to establish the identity of CMs from SP and DN populations in the future. . T test (D, K, L) and uncorrected Fisher's LSD test (E-F, O) were used. ns = non-significant, *P < .05, **P < .01, ***P < .001, ****P < .0001.

Conclusion
Bulk-and single-cell transcriptomic analysis in this study provide insights into the differentiation dynamics of cardiomyocytes and cardiac endothelial cells, 2 important human cardiac lineages. This rich dataset is now available for comparison with in vivo data. The ETV2 fluorescent reporter we generated in hiPSCs allowed the identification of a new subpopulation of early CM precursors that expressed ETV2.