Introduction

Congenital heart disease (CHD) is the most common human birth defect, affecting nearly 1% of live births¹. CHD is generally understood to derive from dysfunction in the growth, morphogenesis, and/or differentiation of cardiac-fated progenitor cell derivatives; however, the genetic and molecular networks governing these processes remain under investigation. Treatment options are limited to surgical interventions that can save the life of the infant but do not correct the underlying genetics. The causative genetic variants may contribute to complications and morbidity in adulthood². Research aimed to expand our understanding of the molecular genetics regulating normal human cardiogenesis is thus needed to advance therapeutic strategies in the detection and treatment of CHD.

Studies of cardiac development in model organisms have led to the discovery of key molecular drivers of heart organogenesis, including the GATA family of zinc-finger transcription factors. GATA factors, named by their ability to bind to the consensus sequence (A/T)GATA(A/G)³, represent six family members with GATA4/5/6 expressed in the developing heart of early stage embryos^4,5. GATA4/5/6 transcription factors regulate multiple stages of heart development including the specification, proliferation, and differentiation of cardiac progenitor cells (CPCs)⁶. Notably, GATA6 knockout (KO) mice exhibit early embryonic arrest during primitive streak formation due to defects in visceral endoderm, a phenotype that has complicated the study of the early developmental function of GATA6 in vivo⁷. Mice with conditional deletion of GATA6 from early CPCs develop heart abnormalities including atrioventricular canal defects, ventricular septal defects, and partial loss of trabeculae^8,9. While these GATA6 conditional homozygous KO mice die at birth due to heart malformations, few cardiac phenotypes were reported in the GATA6 heterozygous mutants, except for recent reports of bicuspid aortic valve and heart rhythm abnormalities^10,11.

In humans, heterozygous mutations in GATA6 are associated with various forms of CHD including outflow tract (OFT) defects, septal defects, and Tetralogy of Fallot^12,13. GATA6 is required for the development of multiple mesoderm and endoderm derived organs and some CHD patients with heterozygous GATA6 mutations have comorbidities including pancreatic agenesis, neonatal diabetes, or congenital diaphragmatic hernia^14,15. The phenotypic diversity of birth defects in patients containing GATA6 mutations illustrates the complexity of human haploinsufficiency; there is likely a combination of variants in modifying or interacting genes that converge to influence the disease phenotype^16,17. Previous studies of the function of GATA6 in model organisms have not sufficiently modeled the variety of cardiac (and non-cardiac) defects observed in humans and most describe phenotypes only after homozygous inactivation of the gene. This raises an important question; why do CHDs in humans associate with GATA6 haploinsufficiency? Species-specific functions for GATA6 may differ between mouse and humans, for example, raising the question of how observations made in model organisms relate to human cardiac development. There is a clear need to investigate the function of GATA6 in the context of a human system to better understand the basis for haploinsufficiency.

Current in vitro cardiac-directed differentiation protocols for human pluripotent stem cells (hPSCs) provide powerful tools for studying the transcriptional regulatory networks relevant to CHD¹⁸. A previous study using human induced pluripotent stem cells (hiPSCs) with CRISPR-Cas9 induced GATA6 mutations showed that GATA6 is required at the stages of CPC and cardiomyocyte (CM) differentiation¹⁹. GATA6 is also required for definitive endoderm (DE) differentiation toward the pancreatic lineage as reported in several studies using hPSCs^20–23. However, the function of GATA6 during early mesoderm patterning of human cardiac differentiation has not been explored, and how in vitro cardiac and DE lineage phenotypes caused by GATA6 loss-of-function relate to one another is not understood.

To determine the function of GATA6 during the earliest stages of cardiac differentiation, we performed in vitro CM directed differentiation using GATA6 loss-of-function human embryonic stem cell (hESC) lines and explored regulatory functions for GATA6 during precardiac mesoderm patterning that is ultimately required for CM generation. Transcriptome analysis revealed that GATA6 mutations impaired expression of genes important for lateral and cardiac mesoderm patterning, and by integrating GATA6 CUT&RUN sequencing with RNA-seq we identified GATA6-regulated targets including those of the WNT and BMP networks, key morphogenetic signaling pathways required to induce anterior primitive streak formation^24,25. Furthermore, we found that GATA6 interacts with important developmental transcription factors and chromatin remodelers during early mesoderm patterning stages that may promote DNA accessibility required for cardiac lineage commitment. Taken together, these findings uncover multiple functions for GATA6 in regulating very early precardiac mesoderm patterning and the lineage specific gene networks required for human CM generation.

Results

GATA6 loss-of-function impairs CM differentiation

A monolayer cytokine-based in vitro cardiac differentiation protocol for hPSCs was adapted from previously established methods^26,27 (Figure 1A). GATA6 heterozygous mutant (GATA6^+/-), homozygous mutant (GATA6^-/-), and isogenic wildtype (WT) control (GATA6^+/+) H1-hESC lines were previously generated^20,28 (Supplemental Figure 1A, B). These lines were differentiated in parallel to determine the impact of GATA6 loss-of-function during human cardiogenesis. Western blotting at days 2 or 5 of cardiac differentiation confirmed the loss of GATA6 protein in GATA6^-/- cells and intermediate levels of protein in GATA6^+/- cells relative to WT controls (Supplemental Figure 1C). The percentage of cardiac troponin T positive (%cTnT⁺) CMs in individual differentiation assays remained stable starting around day 13 and thus CM differentiation efficiency was analyzed during a window of days 13 to 18 of cardiac differentiation. GATA6^-/- cells failed to generate beating CMs and yielded close to 0% cTnT⁺ CMs when examined by flow cytometry (Figure 1B, C). In contrast, GATA6^+/- hESCs did generate beating CMs (∼25% cTnT⁺) but at a significantly reduced efficiency compared to WT cells (∼55% cTnT⁺, Figure 1B, C). GATA6^+/- CMs obtained were comparable morphologically to WT CMs with no obvious sarcomere defects (Supplemental Figure 1D). No significant differences in %cTnT⁺ were observed between clones of the same genotype (Figure 1C), therefore data from lines with the same genotypes were combined in all subsequent comparisons.

Figure 1.

We also generated an iPSC line derived from a CHD patient presenting with an atrial septal defect and congenital diaphragmatic hernia; this patient was found to have a heterozygous frameshift mutation in GATA6 (c.1071delG)¹⁵. The mutant allele was corrected to WT sequence (GATA6^corr/+) using CRISPR-Cas9 gene editing (Supplemental Figure 1E, F). The iPSC colonies expressed normal levels of the pluripotency markers NANOG and SOX2 demonstrating efficient reprogramming (Supplemental Figure 1G). GATA6^1071delG/+ cells generated cTnT⁺ CMs less efficiently than GATA6^corr/+ cells, on average about 50% (Supplemental Figure 1H), indicating that the GATA6 c.1071delG mutation similarly disrupts cardiac differentiation as observed in the hESC heterozygous mutant lines. Likewise, the CMs obtained by differentiation of either GATA6^1071delG/+ or GATA6^corr/+ iPSCs appeared morphologically normal (Supplemental Figure 1I). The corrected iPSC line was shown to have a grossly normal karyotype and to restore higher levels of GATA6 protein as expected (Supplemental Figure 1J, K).

GATA6 is required for CPC specification

The impaired CM generation from GATA6 mutant hESCs may be due to deficient CPC specification. Therefore, RT-qPCR assays were performed on GATA6^+/-, GATA6^-/- and WT cells harvested at day 6 of cardiac differentiation to examine the expression levels of the CPC markers NKX2.5, TBX5, TBX20, and GATA4. CPC marker expression levels were significantly reduced in GATA6^+/- and GATA6^-/- cells relative to WT controls, indicating a defect in CPC specification (Figure 1D). At this time point here was a clear trend but not a statistically significant difference in the expression levels of CPC markers comparing the GATA6^+/- and GATA6^-/- genotypes. However, an RT-qPCR time course showed that CPC marker expression level was increased (although not to WT levels) in GATA6^+/- cells at later stages (days 8 to 11) while GATA6^-/- CPC gene expression levels remained “flat” in comparison with no increase over time (Supplemental Figure 2A).

GATA6^+/+, GATA6^+/-, and GATA6^-/- hESCs were differentiated toward the CM lineage using an alternative protocol based on chemical manipulation of the WNT pathway²⁹ (CHIR protocol, Supplemental Figure 2B). We observed the same phenotype in GATA6^-/- cells and GATA6^+/- cells using the CHIR protocol as with the cytokine-based protocol (data not shown), consistent with a previous report¹⁹. Bulk RNA-seq was performed at day 5 using the CHIR protocol (Supplemental Figure 2C-E). Principal component analysis (PCA) showed that the cells derived from each genotype were well separated from each other (Supplemental Figure 2C). Differential gene expression and gene ontology (GO) analyses revealed that GATA6^-/- cells had strongly reduced expression levels of cardiac development related gene sets including “OFT septum morphogenesis” and “heart development” compared to WT controls (Supplemental Figure 2D). Furthermore, while GATA6^+/- cells had relatively higher expression levels for many of these genes compared to GATA6^-/- cells, numerous CPC regulatory genes showed decreased levels compared to WT controls (including NKX2.5, TBX20, TBX5, MEF2C, and HAND2; Supplemental Figure 2E) confirming an impairment of CPC development in differentiating GATA6 mutant hESCs. GO analysis of genes with increased expression levels in GATA6^-/- cells relative to WT included gene sets such as “nervous system development” (Supplemental Figure 2D) suggesting the adoption of a non-cardiac genetic signature. Notably, gene sets related to retinoic acid (RA) signaling were significantly increased in GATA6^-/- cells relative to WT (Supplemental Figure 2D), and ALDH1A2 expression (the main enzyme that generates RA) was increased in GATA6^-/- cells relative to WT from days 4 to 6 using the cytokine-based protocol (Supplemental Figure 2F), consistent with a previous study of cardiac differentiation using GATA6 loss-of-function hiPSCs¹⁹. Proper RA signaling is crucial for cardiogenesis including OFT development and cardiac chamber morphogenesis³⁰, and inhibiting excessive RA signaling in GATA6 mutant hiPSCs was shown to partially rescue CPC marker expression¹⁹. Thus, upregulated RA signaling in GATA6^-/- hESCs during the CPC developmental stage may partially underly the CM differentiation defects.

GATA6 loss-of-function disrupts cardiac mesoderm patterning

Examining protein expression at multiple stages during cardiac differentiation revealed the earliest (and highest) expression level for GATA6 at day 2 of differentiation (Figure 2A). The impaired CPC development in GATA6^+/- and GATA6^-/- cells may therefore be due to defective cardiac mesoderm patterning during the early stages of cardiogenesis. We examined co-expression of the precardiac mesoderm markers KDR and PDGFRα (KP)²⁷ using flow cytometry from day 3 to 5 of cardiac differentiation (Figure 2B). At day 3, there was no significant difference in the %K⁺P⁺ cells between genotypes (∼45% K⁺P⁺), but by day 5 the %K⁺P⁺ increased to ∼72% in WT and GATA6^+/- cells, while the %K⁺P⁺ was markedly reduced in GATA6^-/- cells (∼22%, Figure 2B, C). This is consistent with previous reports where sorting and reseeding K⁺P⁺ cells at day 3 of cardiac differentiation generated both mesenchymal cells and CMs but day 5 sorted K⁺P⁺ cells gave rise predominantly to CMs, indicating that day 5 KP expression is relatively enriched for precardiac mesoderm^27,31. GATA6^+/- cells had a small but significant decrease in the %K⁺P⁺ cells at day 4 relative to WT (Supplemental Figure 3A, B) suggesting that the developing K⁺P⁺ cardiac mesoderm in GATA6^+/- cells is temporarily defective or delayed relative to WT controls. Notably, the decreased %K⁺P⁺ observed at day 4 (GATA6^+/- and GATA6^-/- cells) and day 5 (GATA6^-/- cells) relative to WT cells was due to reduced KDR but not PDGFRα expression (Figure 2D, E, Supplemental Figure 3A, B), indicating that KDR is (directly or indirectly) regulated by GATA6 during cardiac mesoderm patterning.

Figure 2.

We next sought to determine whether GATA6 loss-of-function impaired mesoderm induction by analyzing the expression of the pan-mesoderm marker BRACHYURY (T). Flow cytometry performed on WT, GATA6^+/-, or GATA6^-/- cells at day 2 or 3 of cardiac differentiation revealed no significant difference in abundance of BRACHYURY⁺ cells between genotypes (Figure 2F, G). Likewise, there was no significant change in T or the mesendoderm marker EOMES transcript levels at day 2 when assessed by RT-qPCR (Figure 2H, I). These results indicate that GATA6 mutant hESCs have the capacity to differentiate into a nascent mesendoderm/early mesoderm population that is deficient in patterning cardiac mesoderm.

Transcriptome analysis of GATA6 mutant hESCs during early mesoderm patterning

During gastrulation, sequentially lineage-restrictive mesodermal cell populations emerge along a developmental trajectory involving primitive streak to lateral mesoderm formation, followed by cardiac mesoderm patterning³². As GATA6 loss-of-function hESCs differentiated to an emergent mesendoderm/mesoderm stage but failed to properly generate cardiac mesdoderm, we hypothesized that GATA6 is required during the earliest stages of lateral and precardiac mesoderm patterning. GATA6-regulated genes important during these early mesoderm patterning stages were investigated by performing bulk RNA-seq at day 2 or day 3 of cardiac differentiation (Supplemental Figure 4A, B). Again, PCA showed the profiles representing each genotype were well separated. GO analysis and gene set enrichment analysis (GSEA) using the biological process (BP) subcollection for differentially expressed genes (DEGs; defined as p-adj<0.05 and log2(fold change)>0.5) comparing GATA6^-/- to WT samples revealed decreased enrichment for GO terms including “heart development” (e.g., HAND1) and “OFT septum morphogenesis” (e.g., NRP1) at days 2 and 3 suggesting a disruption of early cardiac signature genes (Figure 3A, C, D). Conversely, GO analysis of increased DEGs in GATA6^-/- samples relative to WT at days 2 and 3 included “extracellular matrix (ECM) organization”, “skeletal system development”, and “AP axis specification” (Figure 3B, D), suggesting that GATA6^-/- cells are biased to other lineages such as paraxial mesoderm or retain an early mesodermal phenotype by day 3.

Figure 3.

GATA6^+/- samples also had reduced expression levels for genes associated with cardiac GO terms such as “cardiac muscle contraction” and “OFT septum morphogenesis” when comparing to WT controls at day 3 (Supplemental Figure 4C, D). Indeed, there was overlap for many of the decreased DEGs (dDEGs) observed in GATA6^+/- and GATA6^-/- cells at day 2 or day 3 (relative to WT) including cardiac development genes such as ERBB4 (important for endocardial cushion and ventricular trabeculae development³³), the second heart field (SHF) related transcription factor HAND2, and the cardiac mesoderm marker LRRC32³⁴; although the decrease in expression levels for these genes relative to WT was not as severe in the GATA6^+/- genotype as for GATA6^-/- (Supplemental Figure 4D, E, F).

We observed dysregulation of genes related to the WNT and BMP networks in GATA6^-/- samples relative to WT at day 2 (Supplemental Figure 4G) and GO and GSEA analyses revealed decreased levels of “positive regulation of the BMP pathway” and increased levels of “negative regulation of canonical WNT pathway” at days 2 or 3 relative to WT (Figure 3A, B, D, E). Several dDEGs identified in GATA6^-/- cells at day 2 or day 3 relative to WT samples included negative regulators of WNT signaling, including SFRP1 and TLE1 (Supplemental Figure 4E, G), and antagonists of the BMP pathway, such as BMPER and SMAD6 (Figure 3A, E, Supplemental Figure 4E, G), indicating that the overall transcriptional response induced by the WNT and BMP pathways is dysregulated. GATA6^+/- cells also exhibited perturbed gene expression for several WNT- and BMP-related genes (at an intermediate level to those observed in GATA6^-/- cells) at days 2 or 3 relative to WT (Figure 3E, Supplemental Figure 4G). Mesoderm patterning is controlled by the input of WNT, BMP, and TGFβ signaling and dysfunction in these pathways cause severe developmental defects^35,36. Therefore the cardiac defects caused by GATA6 loss-of-function may be caused by these alterations in WNT and BMP signaling.

We examined activity of the BMP, TGFβ and WNT pathways by analyzing the effector proteins SMAD1/5/9 (BMP activated), SMAD2/3 (TGFβ activated), and β-catenin (WNT activated) by western blotting at day 2. There was no change in the levels of active phospho-SMAD1/5/9, phospho-SMAD2/3, or non-phospho-β-catenin when comparing GATA6^+/-, GATA6^-/- and WT cells (Figure 3F). These results indicate that GATA6 regulates the transcriptional response of genes activated by the BMP and WNT pathways downstream of canonical effector protein activation.

To further examine the impact of GATA6 loss-of-function on early mesoderm patterning, we correlated RNA-seq data from a study examining hESC differentiation toward cardiac mesoderm³⁷ to our day 2 or day 3 RNA-seq data. GSEA analysis using increased DEGs at day 2 (lateral mesoderm) or day 3 (precardiac mesoderm) from Koh et al.³⁷ revealed significantly reduced enrichment in our GATA6^-/- samples (relative to WT) on days 2 or 3, respectively (Figure 3G, H). Furthermore, GATA6^+/- cells had partially reduced expression for these lateral and precardiac mesoderm genes relative to WT cells (Figure 3H), suggesting that the impaired CPC and CM differentiation observed for GATA6^+/- hESCs originates as disrupted lateral and cardiac mesoderm patterning.

A previous study of GATA6 function during DE differentiation from hESCs reported an interaction between GATA6, EOMES, and SMAD2/3 to co-regulate genes important during mesendoderm and transition to DE development²². DE and cardiac mesoderm originate from a bipotent germlayer mesendoderm population³⁸ and likely share a gene expression signature at this early stage. We therefore performed GSEA analysis using the CHIP-seq dataset for triple-overlap of GATA6, EOMES, and SMAD2/3 during early DE differentiation from Chia et al.²² and found a significantly decreased enrichment in GATA6^-/- cells relative to WT at day 2 (Supplemental Figure 4H). EOMES is an important transcriptional regulator of anterior mesoderm development³⁹, and as GATA6 was shown to interact with EOMES during early DE differentiation^22,40 we hypothesized that an interaction is similarly necessary during early lateral mesoderm patterning. We therefore performed GSEA using another published dataset of EOMES targets identified during mesoderm and DE (ME) differentiation in mouse ESCs (mESCs)⁴¹. GSEA analysis of our day 2 RNA-seq data with the “EOMES ME direct activation dataset” human gene equivalents revealed a significant decrease in enrichment in GATA6^-/- cells relative to WT (Supplemental Figure 4H), supporting the notion that GATA6 cooperates with EOMES to co-regulate genes required for lateral mesoderm patterning.

GATA6 binds to putative regulatory regions associated with BMP, WNT, and EOMES related dDEGs

To determine the genomic localization of GATA6 occupancy during mesoderm patterning, we performed cleavage under targets and release using nuclease (CUT&RUN) sequencing for GATA6 at day 2 of differentiation using GATA6^+/+ hESCs. 1,273 significant differentially bound sites were identified for GATA6 relative to IgG control samples, associated with 953 transcriptional start sites (TSS) for unique genes (Figure 4A). Significant GATA6 binding sites predominately localized to intronic (46.03%) and distal intergenic (32.21%, >5kb from TSS) regions (Figure 4A), consistent with previous reports of GATA6 binding distribution^42,43. Transcription factor binding motif enrichment analysis revealed the strongest enrichment at GATA6 bound loci for “GATA” motifs, as well as “LHX-like” (LIM homeobox transcription factor; expressed in cardiac and other mesodermal derivatives⁴⁴) and “EOMES” motifs (Figure 4B). GO analysis of the gene list nearby significant GATA6 binding peaks showed striking similarity to GO analysis from RNA-seq dDEGs in GATA6^-/- cells (relative to WT), including “OFT septum morphogenesis”, “BMP signaling pathway” and “WNT signaling pathway” (Figure 4C). Additionally, we observed an overlap of CUT&RUN identified genes with day 2 (17.9%) and day 3 (16.8%) RNA-seq dDEGs in GATA6^-/- cells (relative to WT) representing high likelihood direct targets of GATA6 including the cardiac development related genes HAND1 and MYOCD (Figure 4D). Visualizing the localization of these GATA6 peaks that mapped to day 2 and day 3 dDEGs confirmed binding in intronic regions or at distal intergenic regions that overlapped with known enhancer-like signatures identified by the ENCODE Project⁴⁵ (Figure 4E). These included several WNT-related genes such as MCC and LGR5, as well as BMP-related genes including SMAD6 (Figure 4E).

Figure 4:

Day 2 GATA6 CUT&RUN data revealed an enrichment for binding to sequences containing EOMES motifs (Figure 4B). We compared our day 2 GATA6 CUT&RUN and RNA-seq data to previously published EOMES CHIP-seq data performed on hESCs at day 2 of DE differentiation⁴⁶. Comparing these datasets revealed 32 genes that had triple-overlap of day 2 GATA6 CUT&RUN identified genes (953 genes nearby significantly bound sites), day 2 RNA-seq dDEGs, and EOMES CHIP-seq targets⁴⁶ (Figure 4F). Among this list of triple-overlap genes are known cardiac development factors including LGR5, which was reported to be required for hESC-CM differentiation⁴⁷, and PRDM1, which is expressed in the SHF and involved in OFT morphogenesis⁴⁸. Visualizing the location of these triple-overlap genes confirmed that the GATA6 peaks directly overlapped with many EOMES CHIP-seq peaks (Supplemental Figure 5), suggesting that GATA6 and EOMES interact to co-regulate genes required for precardiac mesoderm patterning.

GATA6 interactome analysis during precardiac to cardiac mesoderm patterning stages

We next analyzed the GATA6 interactome using rapid immunoprecipitation of endogenous proteins (RIME) on WT cells at day 2 or day 4 of cardiac differentiation to define the chromatin associated protein interactions. 99 significant interacting proteins at day 2 and 52 unique proteins at day 4 of cardiac differentiation were identified bound to GATA6 compared to IgG control samples (Figure 5A). EOMES was one of the most enriched proteins bound to GATA6 at day 2 but was not identified in the day 4 GATA6-RIME analysis (Figure 5B). 81.8% of proteins identified by GATA6-RIME at day 2 were not encoded by significantly differentially expressed genes identified by day 2 RNA-seq (GATA6^-/- vs. WT) while 18.2% were, including TLE1, VRTN, and ZIC2 (Supplemental Figure 6A), indicating that GATA6 functions in parallel with most interacting proteins identified. Many of the significantly bound proteins were components of chromatin remodeling complexes such as members of the SWI/SNF complex (e.g., SMARCC1, SMARCA4, and ARID1A), which functions to mobilize nucleosomes and promote DNA accessibility⁴⁹, and the Nucleosome Remodeling and Deacetylase (NuRD) complex (e.g., KDM1A, MTA1/2/3, and GATAD2A), a chromatin remodeling complex associated with transcriptional repression⁵⁰ (Figure 5B). GATA6 protein interactions were also observed for other epigenetic and transcriptional regulators including those related to histone deacetylation (i.e., HDAC2) and mRNA splicing (i.e., SF1, Figure 5B, D). Interestingly, WNT/β-catenin transcriptional regulators were identified bound to GATA6, particularly at day 2 of cardiac differentiation, including the β-catenin transcriptional co-activator EP300⁵¹, and the TCF/LEF transcriptional co-repressors CTBP1/2^52,53 and TLE1⁵⁴ (Figure 5B), suggesting that GATA6 may directly modulate the transcriptional response induced by WNT signaling by binding to these proteins.

Figure 5.

Proteins identified by GATA6-RIME only at day 2 included the transcription factors EOMES and MIXL1 (both transcription factors involved in primitive streak formation and mesoderm patterning^34,55) as well at TLE1, EP300 and SMARCA5, while the unique day 4 proteins were all RNA binding proteins such as TAF15 (Figure 5C). The overlapping day 2 and day 4 GATA6-RIME results included many chromatin remodelers such as SMARCC1, SMARCA4 and KDM1A (Figure 5C), although the enrichment for these proteins was greater at day 2 compared to day 4 (Figure 5B, Supplemental Table 1). As day 2 of differentiation had the highest level of GATA6 expression and enrichment for RIME identified proteins, we performed GATA6 immunoprecipitation (IP) on GATA6^+/+ and GATA6^-/- samples at day 2 and confirmed a weak interaction with EOMES and SMARCC1 in the WT samples, with no observable bands in GATA6^-/- samples (Figure 5E). The EOMES and SMARCC1 protein bands identified in GATA6-IP samples were notably less than observed for the input controls, therefore it’s possible that the interaction between GATA6 and these proteins is indirect (i.e., GATA6 may be co-bound in a complex with EOMES and/or SMARCC1 among other proteins).

Previous studies performed GATA6-RIME at day 4 and EOMES-RIME at day 2 of DE differentiation^21,40, so we compared our day 2 and day 4 GATA6-RIME data to these previously published RIME datasets (Supplemental Figure 6B). 26 proteins were present in all 4 datasets (Supplemental Figure 6B, C), implying there is a common interaction between GATA6 and various remodeler proteins, the transcription factor SALL4, and the TCF/LEF co-repressors CTBP1/2 during mesendoderm patterning toward both cardiac mesoderm and DE lineages. Protein interactions unique to day 2 GATA6- and EOMES-RIME included MIXL1, SMARCC2, and ARID1B and may represent protein interactions specific to mesendoderm development (Supplemental Figure 6C). Notably, EOMES and the WNT-related proteins TLE1 and EP300 were not identified in our day 4 dataset but were present in our day 2 GATA6-RIME data and the published EOMES and GATA6-RIME datasets during early DE development^21,40. Therefore, an interaction between GATA6 and these proteins appears to be relinquished by the cardiac mesoderm stage (day 4) but persists during DE patterning. Cardiac differentiation requires a transient inhibition of the canonical WNT signaling pathway⁵⁶, and as day 4 of our cardiac differentiation protocol occurs after WNT inhibition these varying protein interactions may reflect the differing contributions of the WNT pathway during cardiac and DE differentiation.

Early manipulation of the WNT and BMP pathways partially rescues the CM defects in GATA6 mutant hESCs

LGR5, a component of the WNT pathway⁵⁷, was among the most significant dDEGs identified from day 2 or day 3 RNA-seq in GATA6^-/- cells (relative to WT) and was also found significantly bound by GATA6 (and in a previous study, EOMES⁴⁶) in WT cells (Supplemental Figure 5). As LGR5 is required for efficient CM differentiation in hPSCs^47,58, we hypothesized that infection with a doxycycline (DOX) inducible LGR5-expressing lentivirus might rescue the cardiac defects in GATA6^-/- hESCs. In GATA6^-/- hESCs that had been transduced with the inducible-LGR5 lentivirus (iLGR5), DOX was added from days 1 to 4 of cardiac differentiation (Figure 6A) to mimic the LGR5 expression pattern (Supplemental Figure 7A); we validated that this increased LGR5 expression in a DOX-dose-dependent manner (Supplemental Figure 7B). In all differentiation experiments using iLGR5 treated with DOX (250ng/mL), we did not observe generation of beating CMs. However, DOX treatment of iLGR5-transduced GATA6^-/- hESCs did yield a small but significant increase in the %K⁺P⁺ population at day 5 of differentiation relative to empty vector (EV) transduced controls, associated with increased KDR expression (Figure 6B). Thus, LGR5 expression can partially rescue a proportion of day 5 cardiac mesoderm in GATA6^-/- hESCs but cannot sufficiently improve subsequent CM differentiation.

Figure 6.

As LGR5 is a single component of the complex WNT pathway, we tested a broader pharmacological approach by using the WNT agonist CHIR to determine if early WNT activation could more efficiently rescue the cardiac defects in GATA6^-/- hESCs. We found that inclusion of CHIR (3μM) from days 0 to 2 in differentiating GATA6^-/- hESCs yielded by day 13 small clusters of beating CMs (Supplemental Video 1). When assessing quantitatively by flow cytometry the %cTnT⁺ CMs in CHIR-treated GATA6^-/- cells relative to vehicle treated controls, there was a variable trend for CM induction that was however not statistically significant (Supplemental Figure 7C). Similarly, while the %K⁺P⁺ cardiac mesoderm cells at day 5 was slightly increased upon CHIR treatment in GATA6^-/- cells, this was not a statistically significant improvement (Supplemental Figure 7D).

BMP-related gene expression was also dysregulated at day 2 in GATA6^-/- cells and thus manipulation of the BMP4 concentration in combination with CHIR treatment from days 0 to 2 was tested to further promote CM rescue. Interestingly, we found that decreasing the BMP4 concentration (5ng/mL) combined with CHIR treatment yielded small clusters of beating CMs (Supplemental Video 2) more consistently than increasing the BMP4 concentration (data not shown) or treatment with CHIR on its own. Quantification by flow cytometery for %cTnT⁺ CMs in GATA6^-/- cells treated with CHIR and the lower BMP4 concentration (CHIR LB) showed a small but significant increase (∼4% cTnT⁺) relative to GATA6^-/- controls indicating that CM differentiation is partially rescued upon CHIR LB treatment (Figure 6C). The %K⁺P⁺ and %K⁺ day 5 cardiac mesoderm population was also significantly increased with CHIR LB treatment in GATA6^-/- cells relative to untreated controls and was not significantly different than WT controls indicating a rescue of day 5 K⁺P⁺ cardiac mesoderm (Figure 6D). Taken together, these results indicate that modulating the early input of the WNT and BMP pathways in GATA6^-/- cells can rescue the cardiac mesoderm patterning defects and partially rescue subsequent CM differentiation.

GATA6^+/- cells also exhibited partial dysregulated gene expression related to WNT signaling via day 2 RNA-seq relative to WT controls (Supplemental Figure 4G). We therefore hypothesized that the addition of CHIR from days 0 to 2 in differentiating GATA6^+/- hESCs could rescue CM differentiation more potently than observed in GATA6^-/- hESCs. Indeed, we found that the %cTnT⁺ CMs was significantly increased in GATA6^+/- cells supplemented with CHIR relative to vehicle treated controls and was not significantly different from WT untreated controls, indicating an efficient rescue of CM differentiation (Figure 6E). CHIR addition to WT hESCs had a deleterious effect on CM differentiation efficiency, with the %cTnT⁺ CMs reduced from ∼60% (DMSO treated) to ∼25% (CHIR treated, Figure 6E), suggesting that there is a threshold “window” for the degree of WNT pathway activation required during early primitive streak to mesoderm patterning stages and that reducing or exceeding this threshold has an adverse effect on later CM differentiation capacity. Furthermore, this data indicates that the rescue in CM differentiation efficiency in GATA6^+/- cells treated with CHIR was not due to an overall improvement in the general cardiac differentiation protocol, but rather related to correcting the consequence of loss of one GATA6 allele.

Discussion

We found that GATA6 loss-of-function mutations severely disrupt CPC gene expression and impair CM differentiation, consistent with the results described by Sharma et al.¹⁹. Here, we discovered that GATA6 is required very early in this process and that the later cardiac phenotypes observed in GATA6^+/- and GATA6^-/- hESCs originate from an early mesoderm patterning defect in which lateral and cardiac mesoderm genes fail to be expressed at normal levels. Gene networks induced by the WNT and BMP pathways that are normally required for primitive streak organization and subsequent anterior lateral mesoderm emergence are also mis-expressed when GATA6 is deficient. GATA6 CUT&RUN analysis revealed binding at distal enhancer-like signatures nearby dDEGs identified by RNA-seq, and GATA6 binding was found to have a considerable overlap with previously published EOMES CHIP-seq data during early DE differentiation⁴⁶ nearby many of these genes. Day 2 GATA6 RIME analysis confirmed binding with EOMES and various chromatin remodelers including components of the SWI/SNF and NuRD remodeling complexes that are known to regulate chromatin accessibility of cardiac lineage-specific gene networks^59–63. Furthermore, we found that early manipulation of the WNT pathway through CHIR treatment from days 0 to 2 rescued CM defects observed in GATA6^+/- hESCs while CHIR treatment and reduced BMP4 dosage rescued K⁺P⁺ cardiac mesoderm and partially rescued CM differentiation in GATA6^-/- hESCs. Together, these results present multiple novel functions for GATA6 in regulating early precardiac mesoderm patterning during human cardiac differentiation.

WNT and BMP related genes were dysregulated in GATA6^+/- and GATA6^-/- cells at day 2 and day 3 of cardiac differentiation but β-catenin and SMAD protein activation was intact indicating that the transcriptional response of these pathways is defective rather than the upstream activation induced by morphogens. GATA6 RIME analysis revealed an interaction with several WNT transcriptional regulators including the β-catenin co-activator EP300 and the TCF/LEF transcriptional co-repressors TLE1 and CTBP1/2, and binding to these proteins could have positive or negative effects on their function. One potential mechanism is that GATA6 binds to these WNT regulators and prevents their repressive (or in the case for EP300, potentiating) function on TCF-mediated transcription, causing a net increase in the β-catenin transcriptional response. Indeed, there was an enrichment for negative regulation of the WNT signaling pathway at day 2 in GATA6^-/- cells according to RNA-seq data, suggesting that GATA6 loss-of-function increases feedback inhibition of TCF/LEF-mediated transcription. GATA6 CUT&RUN analysis also revealed binding near many genes related to WNT and BMP signaling implying direct participation in modulating appropriate β-catenin/SMAD transcriptional responses. However, many WNT/BMP response genes dysregulated in GATA6^+/- and GATA6^-/- cells were not identified as targets by GATA6 CUT&RUN, suggesting that another method for regulation is indirect, perhaps related to the interaction with chromatin remodelers affecting DNA accessibility.

RIME analysis revealed GATA6 bound to many chromatin remodelers including components of the SWI/SNF (e.g., SMARCA4) and NuRD (e.g., KDM1A) complexes at days 2 and 4 of cardiac differentiation. GATA6 is a pioneer transcription factor with specialized function to bind and promote accessibility of chromatin during CPC development¹⁹, and GATA6, GATA4, and GATA3 have previously been shown to interact with various members of SWI/SNF complexes such as SMARCA4^19,21,64,65. However, the molecular mechanism for how GATA6 modifies the accessibility of chromatin and mobilizes nucleosomes is unclear. The AP-1 pioneer factor interacts with SWI/SNF complexes to promote their localization to non-permissive regions of chromatin⁶⁶. GATA6 may function similarly, for example, by binding to SWI/SNF complex proteins and recruiting them to regions of closed chromatin containing lineage-specific enhancers where they then function to promote genomic accessibility.

The function of GATA6 during early mesoderm patterning presented in this study has unique implications toward understanding the pathogenesis of CHD and extra-cardiac phenotypes associated with human GATA6 mutations. GATA6 homozygous loss-of-function mutations have not been described in humans, most likely due to early embryonic lethality causing non-viable embryogenesis. This is consistent with the early mesoderm patterning defects of GATA6^-/- hESCs described in this study, previous reports of DE defects in GATA6-null hPSC lines^20–23, as well as the early lethality reported in GATA6 KO mice⁷. In the case of haploinsufficiency, CHD presenting as OFT or septal defects associated with heterozygous GATA6 mutations may be due to defects that occur during early mesoderm patterning that diminish the later developmental pool of CPCs and CMs. Quite remarkably, profiles of cells at day 2 of differentiation already identify OFT morphogenesis as gene sets disturbed by loss of GATA6. Relatively minor changes in early mesoderm cell population numbers have recently been reported to associate with later heart defects in a mouse model of Cornelia de Lange syndrome using NPBL1^+/- mice⁶⁷, and thus GATA6 heterozygous mutations may similarly have a subtle negative effect on developing mesodermal cell populations to cause CHD. Human mutations in GATA6 commonly associate with OFT defects, a primarily second heart field (SHF) derived structure, and recent work has established that SHF lineage progenitors are derived from early defined mesodermal subpopulations⁶⁸. Therefore, GATA6 mutations might impair specific heart field progenitor populations based on the timing and positioning of these developing mesodermal populations in the primitive streak, long before the progenitors are associated with SHF. GATA6 haploinsufficiency in humans also associates with endoderm defects such as pancreatic agenesis, and DE and early mesoderm co-develop from a common mesendoderm population adjacent to one another in the anterior primitive streak during gastrulation⁶⁹. Defects in the regulatory functions of GATA6 during mesendoderm patterning, such as interactions with EOMES and specific chromatin remodelers, may cause downstream organ defects in both mesoderm and endoderm derived organ systems in some patients with GATA6 mutations.

Materials and Methods

hPSC Culture

hPSC lines were cultured on Matrigel (Corning BD354277) coated tissue culture plates using StemFlex (Thermo Fisher Scientific A3349401) or mTESR Plus (Stem Cell Technologies 1000276) medium at 37°C and 5% CO2. Cells were passaged using Accutase Cell Detachment Solution (VWR 10761-312) and supplemented with 10 μM ROCK inhibitor (RI) Y-27632 (Selleckchem S1049) for one day. hPSCs were routinely tested for mycoplasma contamination using the LookOut Mycoplasma PCR Detection Kit (Sigma-Aldrich MP0035).

iPSC generation

A skin biopsy was obtained from a patient with an atrial septal defect and congenital diaphragmatic hernia associated with a heterozygous frameshift mutation in GATA6 (c.1071delG)¹⁵. Human dermal fibroblasts (HDF) derived from this biopsy were transduced with 4 Sendai viruses expressing OCT3/4, SOX2, KLF4, and c-MYC using the CytoTune-iPS 2.0 Sendai Reprogramming Kit (Invitrogen A16517). Transduced HDFs were reprogrammed on CF1 Mouse Embryonic Fibroblasts (MEFs) (Thermo Fisher Scientific A34180) with iPSC medium consisting of DMEM/F12 (VWR 16777-255) supplemented with 20% KnockOut serum replacement (KOSR, Thermo Fisher Scientific 10828010), 1% non-essential amino acids (Thermo Fisher Scientific), 1% L-glutamine (Corning 25-005-Cl), 1% Penicillin:Streptomycin Solution (Corning 30-002-CI), 0.0007% 2-mercaptoethanol (Thermo Fisher Scientific 21985023), 10 ng/mL bFGF (R&D Systems 233-FB-025). Once iPSC colonies emerged, colonies were manually isolated and expanded into feeder-free conditions by culturing with mTESR1 (Stem Cell Technologies 85850) on Matrigel coated plates.

CRISPR Gene Editing

The GATA6 mutant H1 (WA01, NIH #0043) hESC cell lines were established in a previous study²⁰ using the H1-iCas9 cell line and CRISPR gene editing²⁸ targeting the C-terminal ZnF domain (Supplemental Figure 1A, B). The GATA6 c.1071delG iPSC mutant allele was corrected to WT sequence to create the isogenic GATA6^corr/+ line (as well as unmodified GATA6^1071delG+ CRISPR-control clones) using CRISPR technology with Nickase-Cas9 according to previously established methods⁷⁰. A pair of gRNA sequences were designed targeting opposing DNA strands surrounding the G deletion site. The gRNA 1 sequence was specifically designed to overlap the G deletion site to prevent Nickase-Cas9 activity on the WT GATA6 allele. A 91 nucleotide (nt) length single-stranded DNA (ssDNA) repair template (Integrated DNA Technologies) was designed with 45 nt homology arms flanking the G mutation site, including a G to A silent mutation in the right homology arm and gRNA 2 recognition sequence to prevent Nickase-Cas9 activity after homology directed repair (HDR). gRNAs 1 and 2 (Supplemental Table 2) were cloned into the pSpCas9n(BB)-2A-Puro (PX462) V2.0 plasmid (Addgene 62987) as described previously⁷⁰. 2 million dissociated GATA6 c.1071delG/+ iPSCs were electroporated with 4 μg each gRNA-PX462 vector with 3 µL of 100 µmol/L ssDNA donor oligonucleotide using the Neon Transfection System (2 pulses at 1450V for 20 milliseconds, Invitrogen NEON1S) and plated on Matrigel coated dishes with mTESR Plus supplemented with 10 μM RI. A T7 endonuclease I assay (New England Biolabs M0302S) was performed to quantify indel frequency of gRNA pairs following PCR amplification (Phusion Hot Start Flex DNA Polymerase, New England Biolabs M0535L) and agarose gel electrophoresis from genomic DNA isolated 2-3 days after electroporation using the DNeasy Blood & Tissue Kit (Qiagen 69506). 24 hr after electroporation, the culture medium was replaced with mTESR Plus supplemented with 0.5 μg/mL Puromycin (Sigma-Aldrich P8833) for antibiotic resistance selection for 48 hr. iPSCs were then dissociated with Accutase and plated at a low density (∼3000 iPSCs per 10cm plate) on Matrigel coated dishes with mTESR Plus medium (with 10μM RI for the first 24 hr) and allowed to grow for about 2 weeks before iPSC colonies were manually isolated and plated into Matrigel coated 96 well plates and cultured with mTESR Plus. DNA was isolated from iPSC clones using QuickExtract DNA Extraction Solution (Biosearch Technologies QE09050), PCR amplification was performed using primers surrounding the mutation site (Supplemental Table 2) and Sanger sequencing was performed to verify gene editing (GENEWIZ). Karyotyping was performed by WiCell Research Institute (Madison, WI).

hPSC Cardiomyocyte Differentiation

A monolayer cytokine-based cardiac differentiation protocol was adapted from previously established methods (Figure 1A)^26,27. Briefly, a “base medium” for cardiac differentiation was prepared using RPMI 1640 (Thermo Fisher Scientific 22400089) supplemented with 0.5x B27 (Thermo Fisher Scientific 17504044), 0.5 mM 1-Thioglycerol (Sigma-Aldrich M6145), 50 μg/mL L-Ascorbic acid (Sigma-Aldrich A4544), and 1x GlutaMAX (Thermo Fisher Scientific 35050061). The hPSCs were plated at ∼4x10⁴ cells/cm² and allowed to grow for 2 days before changing with “Day 0 medium” consisting of the “base medium” plus 20 ng/mL BMP4 (R&D Systems 314-BP-050), 30 ng/mL Activin A (R&D Systems 338-AC-050), 5 ng/mL bFGF (R&D Systems 233-FB-025), and 5 μL/mL Transferrin (Roche 10652202001). BMP4 was used at 10 ng/ml for iPSC differentiation. Approximately 55 hr later, the medium was changed with “Day 2 medium” consisting of “base medium” supplemented with 5 μM XAV-939 (Sigma-Aldrich X3004), 5 ng/mL VEGF (R&D Systems 293-VE-050), and 5 μL/mL Transferrin. The culture medium was changed on days 4 and 6 of cardiac differentiation with “base medium” supplemented with 5 ng/mL VEGF. From day 9 and later the culture medium was changed with “base medium” every 2-3 days. For a subset of experiments, 3 μM CHIR99021 (Stem Cell Technologies 72054) was included in the “Day 0 medium” and/or the concentration of BMP4 was reduced to 5 ng/mL (Figure 6C-E, Supplemental Figure 7C, D). An alternative CM differentiation protocol based on chemical manipulation of the WNT pathway²⁹ was used for a subset of experiments (CHIR protocol, Supplemental Figure B-E). For the CHIR protocol, hESCs were plated at ∼4x10⁴ cells/cm² and allowed to grow for 2 days in mTESR Plus before changing medium to RPMI 1640 with 1x B27 minus insulin (Thermo Fisher Scientific A1895601) and 6 μM CHIR99021 (Stem Cell Technologies 72054) for 24 hr. Medium was then changed to RPMI 1640 with 1x B27 minus insulin for another 24 hr, before changing medium on day 3 to RPMI 1640 with 1x B27 minus insulin and 5 mM IWP2 (Tocris 3533). On day 5, culture medium was changed to RPMI 1640 with 1x B27 minus insulin once more and starting on day 7 medium was changed to RPMI with 1x B27 every other day. CMs were purified using a lactate-selection protocol to prepare for immunocytochemistry by changing the differentiation medium on day 15 of the cytokine-based protocol to a lactate-selection medium consisting of RPMI 1640 lacking glucose (Thermo Fisher Scientific 11879020) supplemented with 5 mM Sodium L-lactate (Sigma-Aldrich L7022), 213 μg/mL L-ascorbic acid, and 500 μg/mL Human Albumin (Sigma-Aldrich A9731). Lactate-selection medium was refreshed on day 17 of differentiation, after which the medium was restored to “base medium” on day 19 and changed every 2 to 3 days prior to fixation.

Flow Cytometry

Cells were detached from culture plates using 0.5% Trypsin / 0.5mM EDTA (VWR 45000-664) or TrypLE Express Enzyme (Thermo Fisher Scientific 12605010) and resuspended in a FACS buffer consisting of PBS with 10% FBS (for cTnT staining) or 0.5% BSA (for live cell staining). Cells were fixed by pelleting at 300g and resuspended in 0.2% paraformaldehyde (PFA) for 25 min at RT before staining with primary antibody for cTNT (Thermo Fisher Scientific MA5-12960) in FACS buffer with 0.5% Saponin (Sigma-Aldrich S7900-25G) for 1 hr at RT, washed 3x in FACS buffer, and then stained for 1 hr at RT with FACS buffer plus Saponin containing secondary antibody (Invitrogen A21236). Live cells were used for staining the cell surface proteins KP by incubating with FACS buffer containing PE-conjugated KDR antibody (R&D Systems FAB357P-100) and APC-conjugated PDGFRa antibody (R&D Systems FAB1264A) or IgG control antibodies (R&D Systems IC002P, IC002A) for 30 min at RT, and then stained with SYTOX Green Dead Cell Stain (Thermo Fisher Scientific S34860) in FACS buffer for 5 min at RT. BRACHYURY flow cytometry was performed using the Foxp3 Transcription Factor Staining Buffer Set (eBioscience 00-5523-00) according to the manufacturer’s protocol using a PE-conjugated BRACHYURY antibody (R&D Systems IC2085P) or IgG control antibody (R&D Systems IC108P). Flow cytometry was performed and analyzed using an Attune NxT Flow Cytometer (ThermoFisher Scientific) for fixed cells and a BD-Accuri C6 Flow Cytometer (BD Biosciences) for live cells.

RT-qPCR

RNA was isolated using the RNeasy Mini Kit (QIAGEN 74106) and reverse-transcribed using the SuperScript VILO cDNA Synthesis Kit (Invitrogen 11754250) according to the manufacturer’s protocols. The qPCR reaction was prepared by mixing cDNA with LightCycler 480 SYBR Green Master mix (Roche 04887-352-001) and primers listed in Supplemental Table 2. RT-qPCR was performed using a LightCycler 480 Instrument II (Roche). Gene expression data was normalized to endogenous HPRT expression levels as a control for each sample.

Western Blotting

Whole cell lysates were prepared using RIPA Lysis Buffer (Thermo Fisher Scientific 89900) supplemented with 1x Protease/Phosphatase Inhibitor Cocktail (Cell Signaling Technology 5872S) and 1x Benzonase (Millipore 70664-3). Protein was quantified using the Pierce BCA Protein Assay Kit (Thermo Fisher Scientific 23225) and electrophoresed on NuPAGE 4 to 12% Bis-Tris Gels (Invitrogen NP0335BOX) before transferring to PVDF membranes (Bio-Rad 162-0177). Blots were washed with 1x TBS with 0.1% Tween-20 (TBST) and blocked using 5% IgG-free BSA (VWR RLBSA50). Primary antibodies (Supplemental Table 2) were added to blocking buffer and stained overnight at 4°C. Secondary antibody staining was performed for 1 hr at RT using HRP-conjugated antibodies (Supplemental Table 2) in blocking buffer. Pierce™ ECL Western Blotting Substrate (Thermo Fisher Scientific 32106) used used to visualize proteins, and imaging was performed using a C-DiGit Blot Scanner (LICOR) with ImageStudio (v5.2) and Adobe Photoshop (v25.9.1) software.

Co-Immunoprecipitation

Co-immunoprecipitation (Co-IP) was performed using the Pierce™ Crosslink Magnetic IP/Co-IP Kit (Thermo Fisher Scientific 88805) according to the manufacturer’s protocol using cells at day 2 of cardiac differentiation. 7μg of the GATA6 antibody (Cell Signaling Technology 5851) was used to cross-link to the Protein A/G Magnetic beads and 1 mg of protein was used for the IP reaction performed overnight at 4°C on a rotating platform. GATA6 Co-IP samples were electrophoresed on a gel and imaged as described for Western Blotting using 30 μg of protein lysate from the Co-IP used as an input control. GATA6^-/- cell samples served as a negative control.

Immunocytochemistry and microscopy

Cells were fixed in tissue culture dishes with 4% PFA for 20 min at RT. Cells were washed in PBS and then blocked for 1 hr at RT with blocking buffer consisting of PBS with 10% FBS, 0.1% IgG-free BSA (Rockland Immunochemicals BSA-10), 0.1% Saponin, and 10% goat serum (Jackson ImmunoResearch Laboratories Inc. 005-000-121). Cells were then stained with primary antibodies (Supplemental Table 2) diluted in blocking buffer (without goat serum) overnight at 4°C. Secondary antibody staining was performed in blocking buffer (without goat serum) for 1 hr at RT. Cells were stained with DAPI to visualize nuclei and imaged using a Zeiss LSM800 laser scanning confocal microscope with ZEN imaging software or a Zeiss Epifluorescence microscope with AxioVision software, and edited using Adobe Photoshop software (v25.9.1). Videos were taken using an Infinity5 Microscope Camera (Teledyne Lumenera) with Infinity Analyze software.

Lentivirus production and infection

LGR5 full length cDNA was first cloned into a DOX-inducible Lentiviral vector containing a GFP selection marker (CS-TRE-PRE-Ubc-tTA-I2G⁷¹) by performing PCR for LGR5 using WT cDNA samples from day 2 of cardiac differentiation and primers designed to incorporate the BsiWI and MfeI restriction enzyme recognition sites flanking LGR5 (Supplemental Table 2) to allow digestion and insertion into the CS-TRE-PRE-Ubc-tTA-I2G vector. LGR5 sequence was validated by Sanger sequencing (GENEWIZ) using sequencing primers listed in Supplemental Table 2. Lentivirus production was performed by transfecting 293T cells with the lentiviral expression vectors (LGR5 or EV control) and packaging plasmids using PEI. Transfected 293T cells were then allowed to produce lentivirus for 48 hr before the medium was collected, and the virus was concentrated using Lenti-X Concentrator (Takara 631231) according to the manufacturer’s protocol and resuspended in StemFlex medium. GATA6^-/- hESCs were infected with lentivirus by replacing medium with iLGR5 or EV control⁷² lentivirus concentrate in StemFlex for 16 hr before changing to normal StemFlex and seeded for 2 to 3 days. Lentivirus infected hESCs were then sorted based on GFP expression using FACS through the Weill Cornell Medicine Flow Cytometry Core Facility. Doxycline (DOX, Sigma-Aldrich D9891) was added from days 1 to 4 to induce expression of iLGR5 or EV infected GATA6^-/- hESCs.

RNA-seq

Total RNA was extracted from hESCs approximately 52 hr (day 2) or 74 hr (day 3) after starting the cardiac differentiation protocol using the RNeasy Mini Kit (QIAGEN 74106) from at least two independent biological replicates per sample. Total RNA integrity was verified using a 2100 Bioanalyzer (Agilent Technologies), and RNA concentrations were measured using the NanoDrop system (Thermo Fisher Scientific). Libraries were generated from day 2 RNA using the TruSeq RNA Library Prep Kit v2 (Illumina RS-122-2001) and day 3 RNA using the TruSeq Stranded mRNA Library Prep Kit (Illumina 20020594) according to the manufacturer’s instructions. Normalized cDNA libraries were pooled and sequenced on Illumina HiSeq4000 (Day 2, single read) or Illumina NovaSeq 6000 (Day 3, paired-end read) sequencing systems for 50 cycles through the Weill Cornell Medicine Genomics Resources Core Facility. The raw sequencing reads in BCL format were processed through bcl2fastq 2.19 (Illumina) for FASTQ conversion and demultiplexing. Day 2 RNA-seq alignment (GRChg19) and differential gene expression analysis was performed as described⁷³. Day 3 RNA-seq reads were aligned and mapped to the GRCh38 human reference genome by STAR (Version 2.5.2⁷⁴) and transcriptome resconstruction was performed by Cufflinks (Version 2.1.1, http://cole-trapnell-lab.github.io/cufflinks/) after trimming adaptors with cutadapt (Version 1.18, https://cutadapt.readthedocs.io/en/v1.18/). The abundance of transcripts was measured with Cufflinks in Fragments Per Kilobase of exon model per Million mapped reads (FPKM)^75,76. Raw read counts per gene were extracted using HTSeq-count v0.11.2⁷⁷. Gene expression profiles were constructed for differential expression, cluster, and principle component analyses with the DESeq2 package⁷⁸. For differential expression analysis, pairwise comparisons between two or more groups using parametric tests where read-counts follow a negative binomial distribution with a gene-specific dispersion parameter were performed. Corrected p-values were calculated based on the Benjamini-Hochberg method to adjust for multiple testing. GSEA was performed using GSEA 4.2.3 software (Broad Institute, Inc. and Regents of the University of California)⁷⁹ on RNA-seq gene lists ranked by log2(Fold Change) using the subcollection of biological process (BP) gene sets from the GO collection in the Human Molecular Signatures Database (MSigDB)⁸⁰, or using custom gene lists from published sequencing datasets as described in this manuscript^22,37,41, and gene sets with a FDR (False Discovery Rate) q-value < 0.25 and a nominal p-value < 0.05 were considered to be enriched significantly. GO analysis was performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID, david.ncifcrf.gov)^81,82 using the BP subcollection. Heatmaps were generated from normalized counts from selected gene lists using SRplot (bioinformatics.com.cn/srplot)⁸³. Differentially expressed genes (DEGs) were defined as p-value adjusted (p-adj) <0.05 and log2(fold change)>0.5, genes with a count below 50 for all samples were excluded.

CUT&RUN

500,000 H1-GATA6^+/+ cells were harvested approximately 52 hr after starting the cardiac differentiation protocol and CUT&RUN was performed using the CUTANA ChIC/CUT&RUN Kit (Epicypher 14-1048) according to the manufacturer’s protocol. Anti-GATA6 (n=3 independent biological replicates) (Cell Signaling Technology 5851) and Rabbit Anti-IgG (n=2 independent biological replicates) negative control (Epicypher 13-0042) antibodies were used at a 1:50 dilution (∼0.5μg). DNA libraries were generated using the 2S Library Kit (Swift) and sequenced using a NextSeq500 Sequencing System (Illumina) through the Weill Cornell Medicine Epigenomics and Genomics Core Facilities. FASTQ files were aligned to the GENCODE GRChg38 build of the human genome using bowtie2 with the –dovetail parameter. The resulting BAM files were sorted and indexed using samtools. The sorted and filtered BAM files were then converted to bedgraph format using bedtools and peaks were called using SEACR calling enriched regions in target data by selecting the top 1% of regions by AUC (stringent threshold). Transcription factor motif enrichment analysis was performed using Basepair software with HOMER following MACS2 peak calling (basepairtech.com). EOMES CHIP-seq data from Teo et al.⁴⁶ was aligned to GRChg38 using bowtie2 with Partek software (partek.com). CUT&RUN and CHIP-seq data were visualized using Integrative Genomics Viewer (IGV) 2.16.2, distal enhancer-like signatures (dELS) defined by the ENCODE project⁴⁵ were identified using the UCSC Genome Browser for hg38 (genome.ucsc.edu).

RIME analysis

About ∼50 million H1-GATA6^+/+ cells at day 2 or day 4 of cardiac differentiation were fixed with 1% formaldehyde (Electron Microscopy Sciences 15710) for 8 min at RT with gentle agitation. Fixation was quenched with 125 mM glycine (Bio-Rad 1610718) for 5 min at RT before cells were manually scraped, transferred to conical tubes, and pelleted at 800g at 4°C for 10 min. Pellets were resuspended in PBS with 0.5% IGEPAL CA-630 (Sigma-Aldrich I8896) and centrifuged again before resuspending in PBS-IGEPAL with 1 mM PMSF (Thermo Fisher Scientific 36978) and centrifuged one final time. Supernatant was then removed from pellets which were then snap frozen on dry ice and transferred to -80°C. RIME analysis was performed by Active Motif (Carlsbad, CA) according to published methods⁸⁴ using 150 μg of chromatin per sample immunoprecipitating with anti-GATA6 (Cell Signaling Technology 5851) or anti-IgG control (Cell Signaling Technology 2729) antibodies (15 μg antibody per IP) and performing mass spectrometry with enriched GATA6 bound proteins defined as not present in the IgG controls and having a spectral count (the total number of all peptides detected) of at least 5 and >1 unique peptides.

Statistical analysis

Data is represented as the mean ± SEM from independent biological replicates (as indicated in figure legends). Graphing and statistical analysis was performed using GraphPad Prism v.10.2.3 software with two-tailed Student’s t test when comparing two groups. One-way or two-way ANOVA was performed when comparing three or more groups, unless otherwise indicated.

Acknowledgements

JAB, DH, and TE conceived and planned the study. JAB generated most of the data, generated figures, and wrote the first manuscript draft. MG and KMB provided significant assistance and consult on experimental strategies and bioinformatics. RK and ZS generated and validated the iPSC lines. NdS and EY carried out experiments. ZS, KL, and DY generated and provided the hESC mutant lines prior to publication, and provided significant consult. WKC provided patient-derived fibroblasts used for iPSC generation. All authors edited and agreed on the final manuscript version. This work was supported by grants from the National Institutes of Health to TE (R35HL135778), DH (R01DK096239), and WKC (P01HD068250), an NIH postdoctoral fellowship to JAB (F32HL152575), the MSK Cancer Center Support Grant/Core Grant (P30CA008748), the Tri-Insitutional Medical Scientist Training Program to KMB (T32GM007739) and additional grants (TE and DH) from the Starr Tri-I Stem Cell Initiative (2016-004) and the New York State Department of Health (NYSTEM, C029567).

Competing interests

T.E. is a co-founder of OncoBeat, LLC.

Data availability

All RNA-seq and CUT&RUN data is available at the GEO repository and are publicly available as of the date of publication. The series accession number is GSExxxxxx. Requests for primary data, resources and reagents should be directed to and will be fulfilled by the lead contact, Todd Evans (tre2003@med.cornell.edu).