Patient-Specific TBX5-G125R Variant Induces Profound Transcriptional Deregulation and Atrial Dysfunction

Supplemental Digital Content is available in the text.

Tissue harvest. Adult mice were asphyxiated using CO2 followed by cervical dislocation. Pregnant mice used for embryo collection were anaesthetized using 5% isoflurane for 5 minutes followed by cervical dislocation and embryo harvest was performed in 1x PBS at 4°C.
ECGs and Trans-esophageal burst pacing. All procedures as described elsewhere. In this experiment male (control n=22, HET n=10) as well as female mice (control n=12, HET n=8) were used of ages between 12-16 weeks, and the scientist was blinded for genotype during ECG measurements, pacing and data analysis. In short, mice were anaesthetized using 5% isoflurane, and after disappearance of reflexes the animals were placed on a heating pad of 37°C with a steady flow of isoflurane 1.5%.
Electrodes were inserted subcutaneously in the limbs and connected to an electrocardiogram (ECG) amplifier (PowerLab 26T, ADInstruments, Oxford, United Kingdom). ECG was measured for 5 min before further pacing experiments. Electrophysiology studies were performed as described before. 54,55 In short: For atrial stimulation, an octapolar catheter was advanced through the esophagus to the site where capture of the atria could be established. The transesophageal stimulation electrodes were connected to an isolation unit through which triggered stimuli with variable current output were initiated. For sinus node recovery time (SNRT) measurements, a 4 s pacing train with a cycle length of 100 ms was applied. SNRT definition: interval between the last pacing stimulus and onset of first sinus return. To control for differences in sinus rate, SNRT was normalized to resting heart rate (cSNRT = SNRT -RR interval) before pacing. Wenckebach period was determined by incremental increases in atrial pacing frequency performed for 2 s, beginning at 100 ms. Atrial and atrioventricular (AV) nodal effective refractory periods were determined by 1-s stimulus trains at 100 ms cycle length with a shorter last stimulus cycle of 90 ms decreasing successively by 2 ms (s1-s2 loop) until AV block was achieved. Atrial arrhythmia inducibility was determined by 1 or 2 s bursts with cycle length of 60 ms, decreasing successively with a 2 ms decrement -measured down to a cycle length of 10 ms. AA definition: period of rapid irregular atrial rhythm lasting at least 1 s.
Optical mapping. The mice (8 control and 8 Tbx5 +/G125R ) were killed by cervical dislocation after which the hearts were excised and the aorta cannulated. The hearts were Langendorff-perfused with Tyrode's solution (37°C) containing: NaCl 128 mM, KCl 4.7 mM, CaCl2 1.45 mM, MgCl2 0.6 mM, NaHCO3 27 mM, NaH2PO4 0.4 mM, and glucose 11 mM. The solution was maintained at 7.4 pH by equilibration with a mixture of 95% O2 and 5% CO2 and contained 10 μM Blebbistatin (Bio-Techne LtD) to reduce motion artefacts in the optical signals. The hearts were incubated -not perfused -in 10 ml Tyrode's containing 15 μM Di-4 ANEPPS (Molecular Probes, Eugene, OR) after which the hearts were placed in the optical mapping setup and submerged in HEPES buffered Tyrode's solution. Pseudoelectrograms were recorded by placing 3 electrodes at ±0.5 cm distance of the heart in the Einthoven configuration. Electrode R and L were placed alongside the right and left atrium, respectively, whereas electrode F was placed alongside the apex. Electrode R was used as negative input for both lead I and lead II. Recordings were made using Labchart amplifier (AD Instruments, Model 15T, sample frequency 1 kHz) and analyzed in LabChart Pro v8.1.13. Excitation light was provided by a 5 Watt power LED (filtered 510 ± 35 nm). Fluorescence (filtered > 610nm) was transmitted through a lens system on CMOS sensor (100 x 100 elements, sampling rate 5 kHz, MICAM Ultima). Optical action potentials were analyzed with custom made software. Atrial conduction velocity was determined during central stimulation using a custom made current controlled stimulator (Antronics).
Patch clamp. Mouse hearts were harvested from male mice of 8 to 16 weeks of age (n=7 wildtype and n=4 Tbx5 +/G125R ). Single cells were isolated from left atria by enzymatic dissociation. Therefore, excised hearts were perfused for 5 minutes in a Langendorff system with a modified Tyrode's solution containing (in mmol/L): NaCl 140, KCl 5.4, CaCl2 1.8, MgCl2 1.0, glucose 5.5, HEPES 5.0; pH 7.4 (set with NaOH). Subsequently, the hearts were perfused with Tyrode's solution containing a low Ca 2+concentration (10 μmol/L) for 10 minutes, after which Liberase TM research grade (Roche Diagnostics, GmbH, Mannheim, Germany) and Elastase from porcine pancreas (Bio-Connect B.V., Huissen, Netherlands) were added for 12 minutes at a concentration of 0.038 mg/mL and 0.01 mg/mL, respectively. All solutions were saturated with 100% O2 and the temperature was maintained at 37 °C.
To obtain single cells, the digested left atria was cut into small pieces which were triturated for 4 minutes through a pipette (tip diameter: 0.8 mm) in the low Ca 2+ Tyrode's solution, supplemented with 4 stored at room temperature for at least 45 min before they were used. Quiescent single rod-shaped cells with smooth surfaces were selected for electrophysiological measurements.
Action potentials (APs) were recorded at 36±0.2 °C using the amphotericin-perforated patch-clamp technique with an Axopatch 200B amplifier (Molecular Devices, Sunnyvale, CA, USA). Data acquisition and analysis were performed with custom-made software. Patch pipettes (resistance 2-3 MΩ) were pulled from borosilicate glass capillaries (Harvard Apparatus, UK) using a custom-made microelectrode puller. Action potentials were low-pass-filtered with a cut-off of 5 kHz, and digitized at 40 kHz.
Potentials were corrected for the calculated liquid junction. 56 AP measurements. AP were measured using the modified Tyrode's solution as bath solution.. Pipettes were filled with solution containing (in mmol/l): K-gluc 125, KCl 20, NaCl 5, amphotericin-B 0.44, HEPES 10, pH 7.2 (KOH). APs were elicited at 2 to 8 Hz by 3-ms, ~1.2× threshold current pulses through the patch pipette. We analyzed resting membrane potential (RMP), maximal AP amplitude (APA), maximum AP upstroke velocity, and AP duration at 20, 50 and 90% repolarization (APD20, APD50 and APD90, respectively). Parameters from 10 consecutive APs were averaged. albumin, pH 7.3) using the fluorescent probe Indo-1 as described previously. 57 In brief, myocytes were exposed to the acetoxymethyl esters of either 5 µmol/l indo-1 during 30 min at 37°C. They were washed twice with fresh extracellular solution (without albumin) and kept for 15 min to ensure complete de-esterification. Myocytes were attached to a poly-D-lysine (0.1 g/l) treated cover slip placed on a temperature controlled (37°C) microscope stage of an inverted fluorescence microscope (Nikon Diaphot) with quartz optics. A temperature controlled perfusion chamber (height 0.4 mm, diameter 10 mm, volume 30 µL, temperature 37°C), with two needles at opposite sides for perfusion purposes, was tightly positioned over the cover slip. The contents of the chamber could be replaced within 100 ms. Bipolar square pulses for field stimulation (40 V/cm) were applied through two thin parallel platinum electrodes at a distance of 8 mm. One quiescent single myocyte was selected (myocyte with more than one spontaneous oscillation per 10 s were excluded) and the measuring area was adjusted to the rod-shaped surface with a rectangular diaphragm. The wavelength of excitation of Indo-1 was 340 nm, applied with a stabilized xenon-arc lamp (100 W) with flashes of 100 at duration at a frequency of 2 Hz. Fluorescence was measured in dual emission mode at 410 and 516 nm. Emitted light passed a barrier filter of 400 nm, a dichroic mirror (450 nm) and respective narrow band interference filters in front of two photomultipliers (Hamamatsu R-2949). The sample rate of signals recorded during each flash was 10 kHz and averaged. Averaged signals were corrected for background as well as for separately measured signals from non-loaded myocytes. Apparent [Ca 2+ ]i was calculated according to the ratio equation. 58 RNA isolation. RNA isolation of tissue for both qPCR and RNA-seq was done with ReliaPrep RNA Tissue Miniprep System (Promega), using the kit DNase step. RNA quantity and quality was tested using respectively Qubit RNA HS Assay (ThermoFisher) and Bioanalyzer High Sensitivity RNA analysis chip (Agilent).
Nuclei isolation. Nuclei were isolated as described before. 59 All steps were performed at 4 °C. In short, tissue was cut into pieces <1mm on ice using razor blades. The tissue pieces were turraxed for 45 s at highest setting using Turrax in 2 mL Lysis buffer (10 mM Tris-HCl pH 8, 5 mM CaCl2, 2 mM EDTA, 0.5 mM EGTA, 1 mM DTT, 3 mM MgAc) containing RNAse inhibitor (NEB) in a round bottom stoppertube.
Nuclei suspension is then diluted with 3 mL Lysis buffer supplemented with 0.4% Triton-X, and transferred to a Douncer. Next, 10 strokes were performed with a loose pestle, 10 with a tight pestle.
The lysate is passed through a 100 μm strainer (Sysmex), filter washed with 2 mL Lysis buffer supplemented with 0.2% Triton-X, passed through a 30 μm strainer (Sysmex), and the filter washed again with 1 mL Lysis buffer supplemented with 0.2% Triton-X. The resulting lysate is centrifuged for 5 minutes at 1000xg, the supernatant is poured off carefully, and the pellet of nuclei resuspended in Staining buffer (1x PBS with 1% BSA supplemented with RNAse inhibitor).
Droplet-based Sn-RNA-seq. Female mice of 9 to 14 weeks of age were used. Nuclei isolation was performed as described above on 3 right atria of each genotype pooled together loaded (resulting in one wildtype sample (containing n=3 RA) and one Tbx5 +/G125R sample (containing n=3 RA)). Then nuclei were FACS sorted using a gate for DAPI into 0.4% BSA-PBS solution. Nuclei were counted manually using cell counter and 7000 nuclei were loaded. As per 10x Genomics protocol, sorted single nuclei were processed through the GemCode Single Cell Platform using GemCode Gel Bead, Chip, and Library Kits (10x Genomics). The nuclei were then partitioned into gel beads in emulsion in the GemCode instrument followed by processing. Libraries were sequenced on an Illumina Hiseq 4000 PE150. We filtered for at least 1000 reads per nucleus using R2 (R2; Genomics analysis and visualization platform (http://r2.amc.nl)). t-Distributed Stochastic Neighbor Embedding (t-SNE) is a dimensionality reduction algorithm that reduces high dimensional datasets to 2 or 3 dimensions. Differential gene expression in snRNA-seq was determined by Kruskal-wallis test and Bonferroni correction in R2.
Cleavage under targets and release under nuclease (CUT&RUN). 60 Nuclei were isolated as described above from 8 week old wildtype (n=5) and Tbx5 +/G125R (n=6) pooled left and right atrium. 100ul was taken from the isolated nuclei as a negative control for FACS gating. Nuclei were incubated with rabbit PCM-1-antibody (AtlasAB) (1:1000) for 30 minutes rotating at 4 C. Secondary antibody Donkey-antirabbit Alexa 647 (ThermoFisher) (1:500) was added and incubated for a further 30 minutes rotating at 4 C. Then, the sample was centrifuged at 1000xg for 5 minutes and resuspended in 500 μl staining buffer. DAPI was added (1:1000) and the sample was FACS sorted for PCM-1+ nuclei. 3000-5000 nuclei or 6000-30000 nuclei were sorted. Samples were collected in 300 μl staining buffer, centrifuged at 1000xg for 5 minutes, and resuspended in 600 μl NE buffer as stated in step 7 of the C&R procedure. 61 Protocol was followed as published, with 1:00 H3K27Ac incubation overnight. NEBNext Ultra II DNA library Prep Kit for Illumina was used with 15 cycli PCR. Fragment size was determined on a Tapestation high sensitivety chip, DNA concentration was determined with Qubit, and samples were pooled equimolilarly and sequenced PE150 on Hi-seq4000.
Chromatin accessibility analysis using ATAC-seq. Nucleus isolation was performed on control (n=4) and Tbx5 G125R/+ (n=4) atria as described above. We used two male and two female mice of 8-12 weeks for each group. ATAC-seq was performed on FACS-sorted PCM-1 positive nuclei (as described under CUT&RUN) using Illumina Tagment DNA enzyme and buffer (20034197) according to protocol. Libraries were sequenced paired end (PE150) on an Illumina Hiseq4000.
Differential expression analysis (RNA-seq). Reads were mapped to mm10 build of the mouse transcriptome using STAR. Differential expression analysis was performed using the DESeq2 package based on a model using the negative binomial distribution. 50 P-values were corrected for multiple testing using the false discovery rate (FDR) method of Benjamini-Hochberg. We have used 0.05 as FDR control level. Peak calling and differential acetylation analysis (CUT&RUN). Reads from H3K27ac CUT&RUN data were mapped to mm10 build of the mouse genome using BWA, 62 the default settings were used. The BEDTools suite was used to distribute the genome wide signal into bins of 500 bp. 63 The background read level was determined per sample per 500 bp by taking the average read count in bins with less than 20 tags. The background level of tags was subtracted per bin per sample. Bins with less than 101 cumulative tags across all 11 samples were discarded as noise. Differential acetylation was assessed using the DESeq2 package 50 based on a model using the negative binomial distribution P-values were corrected for multiple testing using the FDR. We have used 0.05 as FDR control level. Continuous bins with differential signal were subsequently merged together using the BEDTools suite.
Expression of variant alleles. STAR was used to map reads to mm10 build of the mouse transcriptome. 64 We used VarScan 2 to detect variants and their allele specific occurrence within the RNA-seq. 65 We used a minimum read depth of 100 and a minimum variant allele frequency threshold of 0.20 in VarScan 2 for variant detection. To determine whether the G125R variant was more often present in the RNA-seq reads compared to the wildtype allele, we performed Fisher's exact test on the observed allele frequencies.
ncRNA analysis. The RA RNA-seq dataset was used to calculate de novo transcripts using StringTie, 66 with a minimal transcript length 50bp. To distill the ncRNA from all detected transcripts, genes encoding for proteins (gencode.vM3.annotation on mouse genome build mm10 obtained from Biomart) subtracted from the StringTie result.
The Limma package was used to determine differential expression with settings from referenced paper, as it has shown superior performance in detecting ncRNAs. 67 P-values were corrected for multiple testing using FDR. We have used 0.05 as FDR control level. The resulting list of 159 differentially expressed ncRNAs contains a mix of MGI, RIKEN and not previously annotated transcripts; not all of them will be eRNAs. We considered a differential eRNA detection valid when a clear TSS could be determined (given the signature distribution of forward and/or reverse reads) and when we could establish colocalization with enhancer signatures at the TSS using EMERGE heart enhancer prediction 68 and the CUT&RUN and ATAC-seq data described in this manuscript. We detected 128 eRNAs to be differentially expressed in Tbx5 G125R/+ .
To determine which ncRNAs are potentially relevant for nearby gene expression, promoter coordinates of genes within 2Mb of the differentially expressed ncRNAs were obtained by intersecting these coordinates using the BEDTools suite in Galaxy, 63 and expression results of these genes were extracted from the RNA-seq in this study. Additionally, myoblast TAD coordinates were obtained 69 to determine whether eRNA and potential (differentially expressed) target gene were within the same TAD (using BEDTools suite in Galaxy) (Supplemental Table 17).
Chromatin accessibility analysis using ATAC-seq. CM nucleus isolation was performed on control (n=4) and Tbx5 G125R/+ (n=4) atria as described above. ATAC-seq was performed on FACS-sorted PCM-1 positive nuclei using Illumina Tagment DNA enzyme and buffer (20034197) according to protocol.
Libraries were sequenced paired end on an Illumina Hiseq4000. Peak calling and differential accessibility analysis (ATAC-seq). Reads from ATAC-seq data were mapped to mm10 build of the mouse genome using BWA, 62 the default settings were used. The BEDTools suite was used to distribute the genome wide signal into bins of 500 bp. 63 Bins with less than 89 cumulative tags across all 10 samples were discarded as noise. Differential accessibility was assessed using the DESeq2 package based on a model using the negative binomial distribution. 50 P-values were corrected for multiple testing using FDR. We have used 0.05 as FDR control level. Continuous bins with differential signal were subsequently merged together using the BEDTools suite.
TF motif analysis using HOMER. 200bp summits were determined of ATAC-seq performed in mouse SAN cells. 70 These regions were intersected with H3K27Ac CUT&RUN peaks. In total, HOMER 53 was performed on sequences of 1,714 randomly sampled neutral peaks, 1,393 peaks up in Tbx5 +/G125R , 324 peaks down in Tbx5 +/G125R , all with random genome as background in HOMER.
Nkx2.5 immunostaining and cardiomyocyte counting. Mouse hearts were harvested from female mice of 30 to 36 weeks of age (n=4 wildtype and n=6 Tbx5 +/G125R ), and fixed in 4% PFA for 24 hours. After an EtOH 70% wash, the hearts were embedded in paraffin and sectioned (5um). 1:10 sections with visible atrial tissue were used from immunohistochemistry of Nkx2.5. Sections were also stained for DAPI. Images were taken of 5-15 sections throughout the atria on a DSM5000 with 2.5 objective. In ImageJ, the Nkx2.5 and DAPI positive cells were quantified using automated counting of thresholded (threshold 200) particles. By dividing the number of nkx2.5 positive cells by the number of DAPI positive cells in the left atrium, we calculate the ratio of CMs. The average ratios of CMs per individual were compared using an unpaired t-test.
Fibrosis measurements (picrosirius red stainings). The same processed mouse hearts (n=4 wildtype and n=5 Tbx5 +/G125R ) were used as described in the previous section ("cardiomyocyte counting"). 1:50 sections with visible atrial tissue were stained with picrosirius red. Stained sections were imaged (on a DSM5000 with 2.5x objective) with white balancing, corresponding to an image every 250um. Using ImageJ, (protocol "Quantifying Stained Liver Tissue" https://imagej.nih.gov/ij/docs/examples/stainedsections/index.html) the image is converted to grayscale, and all of the non-atrial tissue (and blood) is blacked out of the picture. The red-stained collagen isolated using thresholding, and measured (for whole tissue measurement threshold 175, for fibrosis 110). Ratios are calculated by dividing percentage area of fibrosis by percentage area of whole tissue. On average 10 images were analyzed per atrium, of which the average was taken per individual.

Legends Supplemental Tables
Supplemental Table 1. Pedigree showing updated information from Postma et al., Circulation Research 2008, with additional electrical anomalies from re-analyzed ECG data.
Supplemental Table 3. Heart size as a function of heart weight divided by tibia length of neonatal day 20 littermates.
Supplemental Table 4. Conduction velocity in isolated right atria.
Supplemental Table 5. Average calcium measurements of isolated atrial cardiomyocytes.
Supplemental Table 6. Single nucleus RNA-seq differential expression in different clusters. Differential expression determined by Kruskal-Wallis test and Bonferroni correction in R2. Table 7. GO term analysis of snRNA-seq differentially expressed genes in cardiomyocyte cluster.

Supplemental
Supplemental Table 8. Whole tissue RNA-seq of the right and left atrium.
Supplemental Table 9. Gene expression of 300 ion channel genes extracted from whole tissue RNA-seq. The 33 significantly different genes are depicted in Figure 5.
Supplemental Table 10. GO-term analysis using PANTHER of upregulated genes in whole-tissue RNA-seq of the right atrium.
Supplemental Table 11. GO-term analysis using PANTHER of downregulated genes in whole-tissue RNAseq of the right atrium.  Table 17. ncRNA expression data. Curated list of ncRNAs differential expression, overlap with epigenetic data (ATAC-seq, H3K27ac and Tbx5 occupancy) and presence of differentially expressed genes in a 1Mb and within the same TAD.

Supplemental Figure 1
Abbreviations: E, embryonic developmental day. (differential) accessibility (ATAC-seq), (differential) H3K27ac, Tbx5 binding in the fetal mouse heart 24 and expression in RA whole adult tissue as well was tag direction of the transcription in Tbx5 G125R/+ and control mice. Tag direction of the transcription in RA is depicted in blue and red as well as arrows in the last two tracks.
Tracks shown ate predicted regulatory regions (EMERGE), (differential) accessibility (ATAC-seq), (differential) H3K27ac, Tbx5 binding in the fetal mouse heart 24 and expression in RA whole adult tissue as well was tag direction of the transcription in Tbx5 G125R/+ and control mice. Tag direction of the transcription in RA is depicted in the last track in blue and red.