Dynamics of genome reorganization during human cardiogenesis reveal an RBM20-dependent splicing factory

Abstract

Functional changes in spatial genome organization during human development are poorly understood. Here we report a comprehensive profile of nuclear dynamics during human cardiogenesis from pluripotent stem cells by integrating Hi-C, RNA-seq and ATAC-seq. While chromatin accessibility and gene expression show complex on/off dynamics, large-scale genome architecture changes are mostly unidirectional. Many large cardiac genes transition from a repressive to an active compartment during differentiation, coincident with upregulation. We identify a network of such gene loci that increase their association inter-chromosomally, and are targets of the muscle-specific splicing factor RBM20. Genome editing studies show that TTN pre-mRNA, the main RBM20-regulated transcript in the heart, nucleates RBM20 foci that drive spatial proximity between the TTN locus and other inter-chromosomal RBM20 targets such as CACNA1C and CAMK2D. This mechanism promotes RBM20-dependent alternative splicing of the resulting transcripts, indicating the existence of a cardiac-specific trans-interacting chromatin domain (TID) functioning as a splicing factory.

Fig. 1

Hi-C across cardiac differentiation. a Schematic of the cardiomyocyte differentiation. b Log transformed contact maps of chromosome 1. c t-SNE plot of PC1 scores on the contact matrices. d Fraction of genome in A and B compartment by sample. e PC1 scores for a region of chromosome 2, gray boxes highlight regions transitioning from A to B and B to A. f Genomic regions divided by stable (81%) and switching (19%) A/B compartment (PC1 scores significantly different by one-way ANOVA, p-value < 0.05; n = 2 independent differentiations). g Regions switching A/B compartment divided by types of transitions. A–B (33%), B–A (49%), A–B–A (8%), B–A–B (10%). Two percent of switching regions were A–B–A–B or B–A–B–A and were combined with A–B and B–A. h Heatmap of the PC1 scores of the compartment switching regions. Clustering of rows based on the four time points of differentiation. Dendrogram of columns was ordered to match the temporal status of differentiation. i Delta compartmentalization saddle plot in cis contacts CM vs. hESC. Bins were assigned to ten deciles based on PC1 score, average observed/expected distance-normalized scores for each pair of deciles were calculated. j Distance plot of A–A, B–B, and A–B interactions for hESC and CM, values are normalized to all contacts at a given distance. Data was smoothed using R, raw maps in Supplementary Fig. 3g. Source data are provided as a Source Data file

Fig. 2

Dynamic gene expression correlations with changes in genome architecture. a PCA plot on the replicates of the RNA-seq samples. b Heatmap of the differentially expressed genes during differentiation, clustering of rows based on differentiation time points. c Enrichment of differentially regulated genes by time point of peak expression against A/B compartment dynamics. Log₂ values are observed/gene density; p-values by chi-squared test for the indicated gene set overlaps, * < 0.05, ** < 0.01, *** < 0.001. d GO term enrichment in CM peak expression genes in B to A compartments, p-value plotted on log scale. e Gene track of Hi-C PC1 and RNA-seq reads of ACTN2. f Gene track of Hi-C PC1 and RNA-seq reads of BMPER. Source data are provided as a Source Data file

Fig. 3

TADs are dynamically regulated independent of A/B compartment changes. a TAD boundaries shared between time points as calculated on the union set of all four time points using both DI and insulation method. b TAD size and number within A and B compartments across differentiation for DI method. Boxplots present the median and 25th and 75th percentile, with the whiskers extending to 1.5 times the inter-quartile range. n = TAD number; p-values by Wilcoxon test, *** < 0.001. c Enrichment of TAD boundaries between hESC and CM state within A/B compartment dynamics for DI method. Log₂ values are observed/TAD union set; p-values by chi-squared test for the indicated TAD boundary set overlaps, * < 0.05, ** < 0.01, *** < 0.001. d Expression of nearest gene to TAD boundaries that are either stage specific or shared between hESC and CM for DI method. Box and whisker plots as in panel (b). n = 308 genes for hESC (lost) boundaries, 207 genes for CM (gained) boundaries, and 1326 genes for shared boundaries; p-values by one-sample, two-sided t-test relative to a hypothesized mean value of 0 (no expression change vs hESC), ** < 0.01. e Gene track of DI score, DI-determined TADs and insulation score-determined TADs and RNA, along with a Hi-C heatmap of the COMMD6 and LMO7 locus in hESC and CM. Source data are provided as a Source Data file

Fig. 4

Changes in local chromatin accessibility occur coincident with global architecture changes. a ATAC peaks across differentiation. b Fraction of ATAC peaks divided by A and B compartment for stage specific and total peaks. c Enrichment of stage-specific ATAC peaks within A/B compartments. Log₂ values are observed/ATAC peak density; p-values by chi-squared test for the indicated ATAC peak set overlaps, * < 0.05, ** < 0.01, *** < 0.001. d Enrichment of peaks within the B to A compartment regions by stage specificity, normalized to ATAC peak density. Plot and statistical analysis as in panel (c). e Motifs from DREME and TOMTOM within CP and CM-specific peaks in B to A compartment bins. f Enrichment of GATA4, NKX2-5, TBX5 (CM), and GATA4 (ESC) ChIP-seq peaks within constitutive A and B regions and B–A and A–B regions. g Overlap of GATA4 peaks with NKX2-5 and TBX5 in B–A regions or the rest of the genome. p-values by chi-squared test for the indicated ChIP-seq peak set overlaps, *** < 0.001. h Metaplot of GATA4, NKX2-5, TBX5 (CM), and GATA4 (ESC) binding across stage specific and constitutive ATAC peaks. Peaks are centered at the mid-point of ATAC peaks and extended ± 1000 bp. i Metaplot of GATA4, NKX2-5, TBX (CM) across CM or CP&CM-specific peaks in B–A regions or the rest of the genome. j Gene track of RNA-seq, ATAC-seq and ChIP-seq reads for NEBL gene and promoter region. Source data are provided as a Source Data file

Fig. 5

TTN locus is spatially regulated during differentiation. a Gene size of upregulated genes peaking in CM stage subdividing either B to A compartment or heart development genes (GO term). Boxplots present the median and 25th and 75th percentile, with the whiskers extending to 1.5 times the inter-quartile range. n = gene number; p-values by Wilcoxon test, ** < 0.01, *** < 0.001. b PC1 value of the TTN locus and the FPKM expression levels across differentiation. c Log transformed contact maps of the long arm of chromosome 2 containing the TTN locus at 500 kb resolution and a zoomed in view at 40 kb resolution, along with the PC1 values. d Gene track of the ATAC-seq signal, CTCF ChIP-seq signal (ESC and CM) and motif orientation, and GATA4, NKX2-5, and TBX5 (CM) ChIP-seq signal across the TTN locus. Dynamic ATAC peaks are highlighted, including increase (in red) at the promoters and decreases (in blue). e Virtual 4C of the TTN promoter at 40 kb resolution. f Fraction of stage-specific ATAC peaks (ESC and CM) overlapping CTCF peaks in the corresponding time point in B–A regions and the rest of the genome. Source data are provided as a Source Data file

Fig. 6

TTN locus becomes associated with RBM20 target genes. a Trans contacts of TTN by compartment, Z-scored by time point at 500 kb resolution. Boxplots within violin plots present the median and 25th and 75th percentile, with the whiskers extending to whole data range. The number of genomic bins is indicated; p-values by Wilcoxon test, *** < 0.001. b FPKM values for RBM20 across differentiation. c, d Association of upregulated RBM20 target genes in hESCs (c) and CM (d). Dashed red line indicates cumulative sum of ICE normalized Hi-C reads between target genes in trans. Histogram represents the background model of 1000 random permutations of selected genes from similar chromosomal distribution. The resulting random shuffling p-value is indicated. e Network of upregulated RBM20 genes, line thickness proportional to contact score, only scores greater than median trans contact displayed. f Representative 3D DNA FISH images in hESC and CM, and 3D reconstructions after spot calling processing (nuclei counterstained with DAPI); scale bars: 5 µm. g On the left, normalized minimum distance per diploid cell between the indicated loci (the number of cells is indicated). On the right, normalized distance of each locus from the nuclear periphery (the number of loci is indicated); A/B compartment transitions based on Hi-C data are reported. Box and whiskers plots present aggregated data from two independent cultures, and indicate median, 25th and 75th percentile, and the 10–90 percentile range. p-values by Kruskal–Wallis test followed by Dunn’s multiple comparisons vs hESC (unless otherwise indicated); ns ≥ 0.05; ** < 0.01; *** < 0.001. Source data are provided as a Source Data file

Fig. 7

TTN trans interactions are transcription-dependent. a Representative results of immunofluorescence for RBM20 and α-actinin combined with DNA FISH for the TTN locus in CM (immunoFISH; nuclei counterstained with DAPI); scale bars: 10 µm. Insets show magnified views. Cells were maintained in standard culture conditions or treated with 5 µM Actinomycin D. b Quantification of distance relationships between RBM20 foci and TTN loci in control conditions. Center-to-center distance below 1 µM (twice the radius of 3D-resconstructed spots) was used to determine overlap (see Methods). On the left, box and whiskers plots showing median, 25th and 75th percentile, and the 10–90 percentile range plus outliers; n = 40 diploid cells. On the right mean ± s.e.m.; n = 5 field of views. c Quantification of RBM20 foci in control and Actinomycin D-treated CM. Graphs are as described for panel (b) except that on the left n = 64 and 55 cells for CTR and ActD, respectively. Note that all CM were considered, including polyploid cells with more than two RBM20 foci. p values by Mann–Whitney test (left) or Welch’s t-test (right); **** < 0.0001. d Representative 3D FISH images of control and Actinomycin D-treated CM; scale bars: 5 µm. e Normalized minimum distance per diploid CM between TTN and the indicated loci (the number of cells is reported). Box and whiskers plots show median, 25th and 75th percentile, and the 10–90 percentile range. p-values by Kruskal–Wallis test followed by Dunn’s multiple comparisons vs Control; *** < 0.001. Source data are provided as a Source Data file

Fig. 8

TTN transcription is required for nucleation of RBM20 into foci. a Schematic of the hESC gene editing strategies used to test the role of transcription at the TTN locus. ∆Prom: promoter (P) deletion; KO: functional knockout by frameshift mutation (indicated by the X). Aspects of TTN biogenesis predicted to be impaired by each modification are indicated by increased transparency (loss of transcription following promoter deletion; loss of translation following knockout). Alternatively spliced exons regulated by RBM20 are shown in blue. b Representative immunofluorescence results in CM derived from the indicated hESC lines (nuclei counterstained with DAPI); scale bars: 10 µm. c RT-qPCR in CM derived from the indicated hESC lines. Expression is relative to the housekeeping gene RPLP0 and mean ± s.e.m. is shown. n = 3, 5, or 7 independent differentiations for TTN KO, TTN ΔProm, and WT, respectively. d As in panel (b); insets show magnified views. e Quantification of RBM20 foci in CM derived from the indicated hESC lines. On the left, box and whiskers plots showing median, 25th and 75th percentile, and the 10–90 percentile range plus outliers; the number of cells is indicated (non-diploid CM were included in the analysis). On the right mean ± s.e.m.; n = 5 field of views. All p-values in this figure are calculated vs WT by one-way ANOVA followed by Holm-Sidak’s multiple comparisons (for (c) and the right graph in (e)) or by Kruskal–Wallis test followed by Dunn’s multiple comparisons (for the left graph in (e)); ns ≥ 0.05; *** < 0.01. Source data are provided as a Source Data file

Fig. 9

Knockout of RBM20 dysregulates alternative splicing in hESC-derived cardiomyocytes. a Predicted effects of the knockout of RBM20 on TTN biogenesis (loss of RBM20-dependent alternative exon exclusion). Refer to the legend of Fig. 8a for the abbreviations. b Representative immunofluorescence results in CM derived from wild-type or RBM20 knockout hESCs generated with two independent CRISPR/Cas9 single guide RNAs (sgRNAs). Nuclei are counterstained with DAPI, and insets show magnified views; scale bars: 10 µm. c Western blot validation of RBM20 knockout (uncropped images shown in Supplementary Fig. 11). g: sgRNA. d RT-qPCR in CM derived from the indicated hESC lines. Expression is relative to the housekeeping gene RPLP0 and mean ± s.e.m. is shown. n = 5 (WT and RBM20 KO g1) or 2 (RBM20 KO g2) independent differentiations. e Expected major TTN isoforms in CM. Exons predicted to be excluded by RBM20 are in blue, and arrows indicate the location of RT-qPCR primers. f As in panel (d), but expression of the indicated TTN isoform is relative to the total amount of TTN. g Expected CACNA1C isoforms due to RBM20-dependent alternative exon selection. RBM20 is predicted to favor inclusion of exon 9 over exon 9*. h As in panel (d), but expression of the CACNA1C exon 9* isoform is relative to the total amount of CACNA1C. i Expected CAMK2D isoforms due to RBM20-dependent alternative exon selection. RBM20 is predicted to favor inclusion of exon 14 over exons 15–16. j Semiquantiative RT-PCR in CM for the indicated CAMK2D isoforms and the housekeeping gene HPRT1; n = 2 independent differentiations. All p-values in this figure are vs WT and calculated by one-way ANOVA followed by Holm-Sidak’s multiple comparisons (ns ≥ 0.05; * < 0.05; ** < 0.01; *** < 0.001). Source data are provided as a Source Data file

Fig. 10

RBM20 foci promote alternative splicing of genes interacting with TTN in trans. a Representative 3D DNA FISH images in CM derived from the indicated hESC lines (nuclei counterstained with DAPI); scale bars: 5 µm. b Normalized minimum distance per diploid CM between TTN and the indicated loci (the number of cells is reported). Box and whiskers plots show median, 25th and 75th percentile, and the 10–90 percentile range. p-values by Kruskal–Wallis test followed by Dunn’s multiple comparisons vs WT; ns ≥ 0.05; ** < 0.01; *** < 0.001. c RT-qPCR in CM derived from the indicated hESC lines. Expression of the CACNA1C exon 9* isoform is relative to the total amount of CACNA1C and mean ± s.e.m. is shown. n = 3, 5, or 7 independent differentiations for TTN KO, TTN ΔProm, and WT, respectively. p-values calculated by one-way ANOVA followed by Holm-Sidak’s multiple corrections vs WT; ns ≥ 0.05; * < 0.05. d Semiquantiative RT-PCR in CM for the indicated CAMK2D isoforms and the housekeeping gene HPRT1; n = 2 independent differentiations. e Proposed model for the regulation of global and local chromatin organization during human cardiogenesis. Upon differentiation the heterochromatin compacts while large cardiac genes such as TTN transition from the inactive to the active compartment. Transcription of TTN nucleates foci of its splicing regulator RBM20 leading to a trans-interacting chromatin domain (TID) involving other RBM20 targets. This mechanism promotes alternative splicing of the resulting transcripts, and can be disrupted by preventing TTN transcription. Source data are provided as a Source Data file