Research ethics for donor tissues
All heart tissue samples were obtained from transplant donors after Research Ethics Committee approval and written informed consent from donor families as previously described2. The following ethics approvals for donors of additional heart tissue were obtained: D8 and A61 (REC reference 15/EE/0152, East of England Cambridge South Research Ethics Committee); AH1 (DN_A17), AH2 (DN_A18), AH5 (DN_A19) and AH6 (DN_A20) (REC reference 16/LO/1568, London, London Bridge Research Ethics Committee); AV1 (HOPA03), AV3 (HOPB01), AV10 (HOPC05), AV13 (HOPA05) and AV14 (HOPA06) (REC reference 16/NE/0230, North East, Newcastle & North Tyneside Research Ethics Committee). Samples of failing hearts used for validation were obtained under the Research Ethic Committee approval given to the Royal Brompton & Harefield Hospital Cardiovascular Research Centre Tissue Bank (REC reference 19/SC/0257).
Tissue acquisition and processing
Cardiovascular history was unremarkable for all donors (Supplementary Table 1). Hearts contributing to the SAN and AVN regions were from donors confirmed to be in sinus rhythm with normal conduction parameters by echocardiogram before donation (Supplementary Table 1). Hearts were acquired after circulatory death (D8, A61, AH1, AH2, AH5 and AV3) and after brain death (AV1, AV10, AV13, AV14 and AH6). For donation after circulatory death donors, after confirmation of death there was a mandatory 5-min stand-off before sternotomy. In all cases, the aorta was cross-clamped and cold cardioplegia solution was administered to the aortic root before cardiotomy. Samples AV1, AV3, AV10, AV13 and AV14 were procured in the standard fashion and then immediately preserved with and transported on a hypothermic perfusion machine. Sample AH2 was similarly preserved, but with immediate normothermic perfusion. It then underwent 4 h of normothermic perfusion before samples were taken. For single nuclei sequencing, all donor samples were full-thickness myocardial tissues from the SAN, the AVN, the LA, the RA, the LV, the RV, the SP and the AX. For spatial transcriptomics, ventricle regions, the thicknesses of which were larger than one side of the Visium frame (6.5 mm), were separated into epicardial and endocardial portions. Samples used for single nuclei isolation and spatial transcriptomics were flash-frozen or frozen in OCT and stored at −80 °C, or formalin-fixed and subsequently embedded in paraffin blocks. All tissue samples were stored and transported on ice at all times until freezing or tissue dissociation to minimize any transcriptional degradation.
Preparation of node samples
For the SAN, a 6 × 3 cm portion of the posterolateral RA with its long axis parallel to and centred on the crista terminalis was dissected. This was then divided into 5-mm-thick strips cut perpendicular to the crista terminalis (Supplementary Fig. 1a). For the AVN region, a tissue sample including the triangle of Koch (bordered by Todaro’s tendon, the coronary sinus ostium and the tricuspid valve annulus) as well as the basal septum (spanning from the interatrial to interventricular septum and including the membranous septum) was dissected (Supplementary Fig. 1b). As before, the sample was then divided into strips that were cut perpendicular to the tricuspid annular plane. Each strip was then embedded in OCT medium and frozen, retaining information of its position in the septum–lateral wall axis. As confirmed by H&E staining, the lateral portions captured AV nodal tissue, whereas the septal portions included the AVB and its branches.
Single nuclei isolation
Single nuclei were obtained from flash-frozen tissues using sectioning and mechanical homogenization as previously described2,54. Slices of 5–10 mm thickness from frozen tissue were first sectioned with a cryostat in a 50-μm thickness section. All sections from each sample were homogenized using a 7 ml glass Dounce tissue grinder set (Merck) with 8–10 strokes of a loose pestle (A) and 8–10 strokes of a tight pestle (B) in homogenization buffer (250 mM sucrose, 25 mM KCl, 5 mM MgCl2, 10 mM Tris-HCl, 1 mM dithiothreitol (DTT), 1× protease inhibitor, 0.4 U μl−1 RNaseIn, 0.2 U μl−1 SUPERaseIn and 0.1% Triton X-100 in nuclease-free water). Homogenate was filtered through a 40 μm cell strainer (Corning). After centrifugation (500g, 5 min, 4 °C), the supernatant was removed and the pellet was resuspended in storage buffer (1× PBS, 4% BSA and 0.2 U μl−1 Protector RNaseIn). Nuclei were stained with 7-AAD viability staining solution (BioLegend), and positive single nuclei were purified by FACS using a MA900 Multi-Application Cell Sorter (Sony) and its proprietary software (Cell Sorter v.3.1.1) (Supplementary Fig. 7). Nuclei purification and integrity were verified by microscopy, and nuclei were further processed for multiome paired RNA and ATAC-seq using Chromium Controller (10x Genomics) according to the manufacturer’s protocol.
Chromium 10x library preparation
Single nuclei were manually counted by Trypan blue exclusion. Nuclei suspension was adjusted to 1,000–3,000 nuclei per microlitre and loaded on a Chromium Controller (10x Genomics) with a targeted nuclei recovery of 5,000–10,000 per reaction. Next, 3′ gene expression libraries and ATAC libraries were prepared according to the manufacturer’s instructions from Chromium Single Cell ATAC and multiome ATAC+Gene Expression kits (10x Genomics). Quality control of cDNA and final libraries was done using Bioanalyzer High Sensitivity DNA Analysis (Agilent) or a 4200 TapeStation System (Agilent). Libraries were sequenced using a NovaSeq 6000 (Illumina) at the Wellcome Sanger Institute with a minimum depth of 20,000–30,000 read pairs per nucleus.
Visium slides and library preparation
For fresh-frozen samples, samples were frozen and embedded in OCT medium using a dry ice-cooled bath of isopentane at −45 °C. OCT-embedded samples were sectioned using a cryostat (Leica CX3050S) and were cut at 10 μm.
For formalin-fixed paraffin-embedded (FFPE) samples, fresh samples were fixed in >5 times their volume of 4% v/v formalin at ambient temperature for 24 h before processing to paraffin on a Tissue-Tek Vacuum Infiltration Processor 5 (Sakura Finetek). FFPE blocks were sectioned at 5 μm using a microtome (Leica RM2125RT).
Samples of the microanatomical regions of interest (ROIs) were selected on the basis of morphology with expert review (S.Y.H.), orientation (based on H&E staining) and either RNA integrity number (fresh-frozen samples) or DV200 (formalin-fixed) that was obtained using an Agilent 2100 Bioanalyzer. In addition, FFPE tissues were checked for possible detachment issues using 10x Genomics Adhesion test slides. FFPE Visium Spatial Gene Expression (10x Genomics) was performed following the manufacturer’s protocol. For fresh-frozen samples, the Tissue Optimization protocol from 10x Genomics was performed to obtain a permeabilization time of 45 min, and the FF Visium Spatial Gene Expression experiment was performed as per the manufacturer’s protocol (10x Genomics). H&E-stained Visium Gene Expression slides were imaged at ×40 magnification on a Hamamatsu NanoZoomer S60. After transcript capture, the Visium Library Preparation protocol from 10x Genomics was performed. Eight cDNA libraries were diluted and pooled to a final concentration of 2.25 nM (200 μl volume) and sequenced on 2× SP flow cells of an Illumina NovaSeq 6000.
Read mapping
After sequencing, samples were demultiplexed and stored as CRAM files. Each sample of sc/snRNA-seq was mapped to the human reference genome (GRCh38-3.0.0) provided by 10x Genomics and using the CellRanger software (v.3.0.2) or STARsolo (v.2.7.3a) with default parameters. For single nuclei samples, the reference for pre-mRNA was created using the method provided by 10x Genomics (https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/advanced/references). Each sample of multiome, or Visium, were mapped to the human reference genome (multiome: GRCh38-2020-A-2.0.0; Visium: GRCh38-3.0.0) provided by 10x Genomics using CellRanger ARC (v.2.0.0) or SpaceRanger (v.1.1.0) with default parameters. For Visium samples, SpaceRanger was also used to align paired histology images with mRNA capture spot positions in the Visium slides. Part of the SAN samples were mixed with different donors after the nuclei isolation for cost-efficient experimental design (Supplementary Table 8) and computationally demultiplexed (Soupercell, v.2.0)55 on the basis of genetic variation between the donors.
Quality control and processing of data
For sc/snRNA-seq and multiome gene expression data, the CellBender algorithm56 (remove-background) was applied to remove ambient and background RNA from each count matrix produced using the CellRanger pipeline. Python (v.3), Pandas (v.1.3.5), NumPy (v.1.21.5), Matplotlib (v.3.5.2) and Scanpy (v.1.8.2 and v.1.9.1) were used for quality control and downstream processing. Cells or nuclei for each sample were filtered for more than 200 genes and less than 20% (cells) or 5% (nuclei) mitochondrial and ribosomal reads. A Scrublet57 (v.0.2.3) score cutoff value of 0.3 of was applied to remove doublets. The Scanpy toolkit was used to perform downstream processing.
For multiome ATAC data (10x Genomics), the data processed using CellRanger ARC were further analysed using ArchR (v.1.0.2)58. Quality control was performed, considering, among other factors, transcription start site enrichment, nucleosomal banding patterns, the number and fraction of fragments in peaks, reads falling into ENCODE blacklist regions as well as doublet scores computed by ArchR. For high-quality cells, reads were mapped to 500-bp bins across the reference genome (hg38) (TileMatrix). Gene scores based on chromatin accessibility around genes were computed from TileMatrix using the createArrowFiles function to check their consistency with measured expression values. Before peak calling, pseudo-bulk replicates were generated (addGroupCoverages) for each fine-grained cell state annotated using the paired gene expression data. Peak calling (501 bp fixed-width peaks) was performed for each cell state, and the peak sets were merged to obtain a unified peak set (addReproduciblePeakSet). A cell-by-peak count matrix was obtained using the addPeakMatrix function.
For Visium data, the Scanpy toolkit was used for quality control and downstream processing. Visium spots of each sample were filtered for more than 500 UMI counts and 300 genes.
Data integration and cell-type annotation
All transcriptome data were integrated using scVI59 (v.0.14.5, n_hidden=128, n_latent=50, n_layers=3, dispersion=‘gene-batch’) and scArches60 (v.0.5.5, n_hidden=128, n_latent=50, n_layers=3, dispersion=‘gene-batch’) with correcting batch effects (donor, cells or nuclei, and 10x Genomics library generation kits) and removing unwanted source of variations (total counts, per cent mitochondrial genes and per cent ribosomal genes for ‘continuous_covariate_keys’). Scanpy functions were used to compute a neighbourhood graph of observations based on the scVI latent space (scanpy.pp.neighbors) and to perform dimensionality reduction (scanpy.tl.umap) and Leiden clustering (scanpy.tl.leiden, resolution = 1.0). Clusters showing hybrid transcriptional signatures that also had a high scrublet score were removed. After re-clustering, cell lineages were annotated on the basis of the expression of major marker genes and statistically identified marker gene expression for each cluster (scanpy.tl.rank_genes_groups).
To identify fine-grained cell states of the CCS, aCMs of the SAN were subclustered; thus we identified a cluster of P cells (SAN_P_cell) that expressed canonical channel genes (HCN1 and HCN4)61 and a TF (TBX3) (Fig. 1c and Extended Data Fig. 2a). Subclustering of AVN aCMs and vCMs showed two CCS cell state clusters (AVN pacemaker cell; AVN_P_cell and AVB cell; AV_bundle) (Fig. 1d). AVB cells formed a distinct cluster defined by an enrichment in GJA5 (which encodes Cx40), CRNDE and CNTN5, which were previously identified as a marker of His bundle cells in the mouse heart9 (Extended Data Fig. 2c). To identify Purkinje cells, CM populations from AVN samples were integrated and clustered with those from the AX. This analysis showed one cluster (cluster 8) that contained not only the CCS cells from the AVN but also a population derived from the AX expressing Purkinje cell marker genes (GJA5, IRX3, KCNJ3 and MYL4)62 (Fig. 1e and Extended Data Fig. 2d,e).
Cell states of other cell types and other regions were defined by label transfer (scNym)63 using a published dataset2 as a reference with revised annotations. The new annotations include neural cell populations, which express pan-glial markers and lack core neuronal markers (Extended Data Fig. 2i); therefore this compartment is best described as glia and will be described below with the ‘_glial’ suffix. FB4 was renamed as FB4_activated based on FB-activation signature genes (POSTN and TNC) and genes encoding ECM proteins (COL1A1, COL1A2, COL3A1 and FN1) (Extended Data Fig. 2h,j). vCM3 was renamed as vCM3_stressed based on the specific expression of NPPB, which encodes B-type natriuretic peptide (BNP), which is a diagnostic marker for HF and a valuable prognostic predictor for the entire spectrum of disease severity and expressed in stressed CMs64 (Extended Data Fig. 2h,k). EC7_atria was renamed as EC7_endocardial based on a recently published study45. For myeloid cells, dimensionality reduction and batch correction (scVI) with Leiden clustering were repeated to identify and annotate fine-grained cell states, such as tissue-resident LYVE1+ MPs65, based on the markers (Supplementary Fig. 2c,d). The transferred cell state labels that were not consistent with the coarse-grained cell-type labels based on the global clusters were replaced with ‘unclassified’ and excluded from downstream analyses.
snATAC-seq data were integrated using cell-by-peak count matrix and peakVI66 (v.0.19.0) with correction for batch effects (donor). Scanpy functions were used to compute a neighbourhood graph of observations based on the peakVI latent space (scanpy.pp.neighbors) and to perform dimensionality reduction (scanpy.tl.umap).
Spatial mapping of cell states with cell2location
To spatially map heart cell states defined by single-cell transcriptomics data analysis in the Visium data, we used our cell2location (v.0.1) method11,67. In brief, we first estimated reference signatures of cell states using sc/snRNA-seq data of each region and a negative binomial regression model provided in the cell2ocation package. For the cell types that had fewer than 100 cells or nuclei per region, cells or nuclei from all the regions were used for the reference signature inference. The inferred reference cell state signatures were used for cell2location cell-type mapping for corresponding regions that estimate the abundance of each cell state in each Visium spot by decomposing spot mRNA counts. The H&E images of the Visium slides were used to determine the average number of nuclei per Visium spot (n = 7) in the tissue and used as a hyperparameter in the cell2location pipeline. For Visium–FFPE sections (Extended Data Fig. 6c,d), cell state proportions in each Visium spot were calculated based on the estimated cell state abundances.
Cell state spatial enrichment analysis
Anatomical microstructures of spatial transcriptomics data were manually annotated using the paired histology images as follows: epicardium, subepicardium, endocardium, myocardium, vessel, nerve, adipose tissue, cardiac_skeleton, fibrosis, node, AVB and Purkinje cell. Cell state proportions per spot were calculated based on the estimated abundance of cell states (cell2location). Cell state enrichments (odds ratio) in each structure were calculated by dividing the odds of target cell state proportions by the odds of the other cell state proportions. Odds of cell proportions were calculated as the ratio of cell proportion in the spots of a structure of interest to that in the other spots. Significance was obtained by chi-square analysis (scipy.stats.chi2_contingency) and the P value was corrected using the Benjamini–Hochberg method.
The mapping of expected cell types to histologically defined structures such as EC7_endocardial in the endocardium, glial cells (NC1_glial) in the nerve, arterial smooth muscle cells (SMC2_art) in the vessel, and FBs and MPs in fibrosis structures (Supplementary Fig. 3a–d) provided further validation of the spatial mapping method.
Cellular microenvironment discovery
The NMF analysis implemented in the cell2location pipeline was performed on the spatial mapping results of each anatomical region of the heart. The NMF model was trained for a range of cell state combinations (number of factors: n_fact) N={5,…,14}, and the effect size of the cell state group abundance between the spots within a given structure against the spots in the other areas was calculated for each factor (95 factors in total) (Supplementary Fig. 4, step 1). To test the significance, we permuted the annotation labels of all spots and generated a null distribution of the effect size. The P values were calculated on the basis of the proportion of the value that is as high as or higher than the actual effect size. For each given structure, we first selected the factor that has the highest significant effect size (best-factor) (Supplementary Fig. 4, step 2). Next, we selected the n_fact that had multiple numbers of factors (fine-factors) with an effect size more than an arbitrary proportion (0.5) of the best-factor (Supplementary Fig. 4, step 3). We considered the fine-factors as refined microenvironments, which were identified using the NMF method (and not with the knowledge-based structural annotations).
For the myocardial-stress niche, first the estimated abundance values (cell2location) of FB4_activated and vCM3_stressed cells were multiplied. Based on the multiplied abundance, clusters of neighbouring spots (n > 5) that had higher than a value of (0.03) were selected.
CellPhoneDB neural–GPCR expansion module
Using the HUGO Gene Nomenclature Committee (HGNC)68 library of GPCRs as a master list (HGNC group 139), we used publicly available databases (UniProt, Reactome, IUPHAR and GPCRdb (https://gpcrdb.org/)69) to generate a set of GPCRs with known ligands. To generate ligand–receptor interactions for these GPCRs, we used ligand genes (for gene-encoded ligands). For non-gene-encoded ligands (such as small-molecule ligands), we used ligand proxies in the form of specific biosynthetic enzymes or transporter genes. Additionally, we added new trans-synaptic adhesion molecule interactions29,70,71. Together, this formed more than 800 new interactions (Supplementary Tables 3 and 4), which we used with the user-defined ‘database generate’ function in CellPhoneDB7.
Spatially resolved cell–cell interaction analysis
CellPhoneDB7,72 analyses with a custom neural–GPCR expansion module were performed on the identified niches and the cell state components. Overall, the cell–cell interaction inferences were performed using single-cell transcriptomics data of each anatomical region and by restricting to the cell states that were colocalized in the identified cellular niches.
For CCS and myocardial-stress niches, we selected the cell states that were either in the node of the SAN or the AVN and retrieved the interacting pairs of ligands and receptors that satisfied the following criteria: (1) all the members were expressed in at least 10% of the cells in the cell states; and (2) ligand–receptor complexes specific to two cell states were inferred by the statistical method framework in CellPhoneDB (‘statistical_analysis’, P value threshold = 0.05).
For epicardial–subepicardial niches, the ligand–receptor interactions of the colocalized cell states were retrieved on the basis of the following criteria: (1) all the members were expressed in at least 10% of the cell states; and (2) at least one of the members in the ligand or the receptor was a differentially expressed gene (DEG) compared with other cell states (scanpy.tl.rank_genes_groups, P value threshold = 0.05, log2(fold change) threshold = 0.1). The ligand–receptor interactions were further selected on the basis of mean expression levels and the biological questions as indicated in the Results and the figure legends.
The following HGNC annotations were used for selecting some of the ligand–receptor interactions: chemokines (HGNC GID:189), cytokines (HGNC GID:599, 602, 1932, 781 and 1264), and LGIC (HGNC GID:GID161).
The following cell states were used in each analysis: SAN (SAN_P_cell, FB2, FB4_activated, FB5, FB6, NC2_glial_NGF+ and LYVE1+IGF1+ MP); AVN (AVN_P_cell, aCM2, FB1, FB2, FB5, SMC1_basic, SMC2_art, NC2_glial_NGF+, LYVE1+IGF1+ MP and mast); epicardial–subepicardial niche (meso, EC6_ven, EC8_ln, PC2_atria, LYVE1+IGF1+ MP, B_plasma, T/NK_cycling, FB2, FB3, FB5 and NC1_glial); and myocardial-stress niche (FB3, FB4_activated, FB5, vCM3_stressed, EC2_cap, EC3_cap, EC4_immune, EC6_ven, PC3_str, SMC2_art, LYVE1+IGF1+ MP, MoMP and NK_CD16hi).
Ion channel and GPCR profile
Differential gene expression analysis with t-test method was performed using the Scanpy function scanpy.tl.rank_genes_groups. Only multiome gene expression data were used to avoid the technical batch effects due to kit differences. P value correction was performed using the Benjamini–Hochberg method. Each of the CCS cells (SAN_P_cell, AVN_P_cell, AV_bundle_cell and Purkinje cells) was compared with non-CCS aCMs as a reference (Supplementary Table 2). Genes were deemed differentially expressed with an adjusted P value of < 0.05. DEGs encoding for ion channels and GPCRs were selected based on HGNC groups 177 and 139, respectively. Upregulated (log2(fold change) > 0) DEGs were depicted in the GPCR overview schematic (Fig. 3b). To compare the transcriptional similarity of working (aCM, vCM) and CCS cell states, a dendrogram based on principal-component-analysis-reduced gene expression was computed using the Scanpy function scanpy.tl.dendrogram (Extended Data Fig. 8c).
Mouse DEG analysis
A list of upregulated ion channel genes (log2(fold change) > 1, adjusted P value < 0.01) for SAN_P_cells was made. Differentially expressed testing summary statistics from two mouse single-cell studies were obtained52,73. In both studies, differentially expressed testing was conducted by comparing sinoatrial CMs (P cells) against all other cells. Genes orthologous to the list of upregulated (human) P cell genes were identified in the mouse differentially expressed summary statistics (using the NCBI Orthologs database as a reference). The differentially expressed statistics (log2(fold change), –log(P value)) from the mouse are plotted. Human and mouse data were not integrated.
Identification of ligands in Visium data
Four spatial transcriptomics, sections were identified as containing the RAGP by an expert anatomist (S.Y.H.). Spot counts from the Visium–FFPE sections that contained RAGP were normalized, and spots were scored for the expression of four generic pan-neuronal cytoskeletal markers (PRPH, NEFL, NEFM and NEFH) using the SCANPY sc.tl.gene_score() function. Correlation of individual gene expression profiles with this score was calculated (Pearson r and P value for each gene). The ligand/ligand-proxy list created as part of the CellPhoneDB neural–GPCR module was used to identify ‘ligand’ genes among the set of correlated genes.
Gene regulatory network
The Scenic pipeline74,75 was used (pySCENIC, v.0.11.2) to predict TFs and putative target genes regulated in P cells. First, gene regulatory interactions were calculated based on co-expression (either positively or negatively correlated) across the single-cell transcriptomics datasets of aCMs (using only multiome data) with GRNBoost2 (ref. 76). This was followed by pruning interactions using known TF-binding motifs and the construction of dataset-specific regulatory modules (regulons)77. Regulons were then scored in each individual cell using AUCell. P-cell-relevant TFs and target genes were retrieved based on the following criteria: (1) regulator TFs that are DEGs in P cells compared with other aCMs (scanpy.tl.rank_genes_groups, with only multiome data, P value threshold = 0.05, log2(fold change) threshold = 0.5); (2a) for activating regulons, target genes that were expressed in at least 10% and differentially expressed in P cells compared with other aCMs (same criteria as TF selection); and (2b) for repressor regulons, target genes that were expressed at low levels in P cells compared with other aCMs (P value threshold = 0.05, log2(fold change) threshold = −0.5). A network of regulatory TFs and target genes was then constructed by linking individual regulons to create a graph (NetworkX, v.2.6.3) (Fig. 3c and Extended Data Fig. 8i). The node colour of the target genes is based on the class GPCR (HGNC and GID139), ion channel (HGNC and GID177) or TFs78.
The interactions of regulatory TFs and target genes were also inferred using the snATAC-seq data and ArchR (v.1.0.2)58. TF-binding motifs in the identified peaks were searched (addMotifAnnotations, motifSet=“cisbp”), and correlations between peak accessibility and gene expression were analysed (addPeak2GeneLinks and getPeak2GeneLinks, correlation > 0.2 or < −0.2, FDR < 1 × 10–4) using the multiome data of aCMs. TFs and potential target gene interactions were obtained by combining the two results (Supplementary Fig. 5c). The inferred interactions are highlighted in red in the activator network graph (Fig. 3c) and blue in the repressor network graph (Extended Data Fig. 8i).
GWAS SNP enrichment analysis
To create a list of SNPs associated with various physiological and pathological cardiovascular traits, index SNPs (meeting genome-wide significance) were first extracted from the NHGRI-EBI GWAS Catalog79 using the R package gwasrapidd80. Index SNPs were added to by SNPs in tight linkage disequilibrium (r2 > 0.8) based on 1000 Genomes (phase 3) European samples, obtained using the Ensembl API80,81, with the window size set to the default 500 kb.
For each nucleus barcode in the snATAC-seq dataset, counts of each identified peak was binarized (1 if the read count was >0). The binarized barcode by peak matrix was then aggregated by cell state (defined using the paired gene expression data) to form a cell state-by-peak matrix, in which a peak was defined as open for a given cell state if the binarized count was 1 in at least 5% of that population.
To calculate enrichment, a permutation test was performed using a previously described method82. For each cell state, a random background was created by shuffling the open/closed labels of the peaks of that cell state such that a random set of peaks were annotated as open, the number of which equalled the actual number of open peaks for that cell state. This was repeated to create 1,000 random permutations. For each trait and cell state, the proportion of trait-associated SNPs falling within the open peaks of that cell state was calculated (the SNP proportion). This SNP proportion was also calculated for each of the 1,000 random permutations. A P value could then be calculated as the fraction of times the random SNP proportion exceeded or was equal to the real SNP proportion. Finally, these P values were corrected for multiple testing using the Benjamini–Hochberg method.
Drug2cell
Drug and target gene information was obtained from ChEMBL83,84 (v.30). Drugs were filtered based on targeting organisms (Homo sapiens), achieved phase in a clinical trial (max_phase=4, clinically approved), and functionally active or not. The activity (pChEMBL) threshold was specifically set for each family of target molecules according to the IDG project (https://druggablegenome.net/ProteinFam) (kinases: ≦30 nM; GPCRs: ≦100 nM; nuclear hormone receptors: ≦100 nM; ion channels: ≦10 μM; others: ≦30 nM) (Supplementary Table 9). Clinically approved drugs were categorized based on the WHO ATC classification (https://www.who.int/tools/atc-ddd-toolkit/atc-classification). Drug scores in each single cell were calculated based on the target gene expression levels. Score were obtained by taking the mean of a set of target genes either without (method=’mean’) or with (method=’seurat’) subtracting with the mean expression of a reference set of genes85. The reference set was randomly sampled from the gene pool for each binned expression value. For the drug repurposing analysis, all the drugs tested and selected were ones that had the statistically highest score in a cell type of interest by testing with Wilcoxon sum test, and P values were adjusted using the Benjamini–Hochberg method. For the drugs targeting GPCRs or ion channels, we searched the clinically approved (maximum phase: 4) and preclinical bioactive molecules with drug-like properties (maximum phase: 1–3) that target genes encoding GPCRs (HGNC GID:139) or ion channels (HGNC GID:177). The drug2cell Python package is available at GitHub (https://github.com/Teichlab/drug2cell).
In vitro validation of chronotropic effects of GLP-1
Cell culture
hiPSC-derived CMs (hiPSC-CM) were differentiated and maintained as previously described86,87. In brief, IMR90 hiPSCs from WiCell88 were seeded onto plates coated with Matrigel (Corning, 356231) in TeSR-E8 medium (StemCell Technologies, 05990) with 10 μM Y-27632 dihydrochloride (Sigma-Aldrich, Y0503). The next day, the medium was switched to TeSR-E8 without Y-27632. Media were subsequently changed daily. Before starting the differentiation into CMs, the cells were replated twice using 0.5 mM EDTA (Thermo Fisher, 15575020, diluted 1:1,000 in PBS (Gibco, 20012027)) at room temperature for 6 min to detach them before plating. For the differentiation process, cells that reached 90% confluency were treated for 2 days with 5 μM CHIR-99021 (Tocris, 4423) in RPMI 1640 (Gibco, 11875119) supplemented with B27, minus insulin (Gibco, A1895602). On day 2, the cell culture medium was replaced with RPMI/ B27, minus insulin. For the next 2 days, 2 μM Wnt-C59 (Biorbyt, orb181132) was added in RPMI/B27, minus insulin. Cells were then cultured in RPMI/B27, minus insulin, with media changes every 2 days. hiPSC-CM contraction was observed between days 8 and 10. On day 11, the cells were placed in starvation medium (RPMI without glucose (Gibco, 11879020) supplemented with B27 (Gibco, 17504044)) for 2 days to improve purity. On day 15, the cells were detached using TrypLE select enzyme (10×) (Life Technologies, A1217702) and replated at a density of 2 × 106 per well in RPMI/B27, 10% KOSR (Thermo Fisher, 10828028) and 10 μM Y-27632 dihydrochloride. The cells were then cultured in RPMI/B27, and media were changed every 2 days.
Gene expression analysis of hiPSC-CMs
Bulk RNA-seq data from IMR90-derived CMs, deposited into the Sequence Read Archive public repository with accession number PRJNA629893, were analysed as previously described89. Transcripts per million (TPM), indicating the normalized amplitude of gene expression, were used to examine gene expression. The TPM was calculated by dividing the read counts by gene length and the total number of exon reads, and then multiplied by the scaling factor of 1,000,000 (ref. 90).
Calcium imaging and quantification
hiPSC-CMs were stained for calcium imaging 30–35 days after the differentiation protocol was started. Cells were seeded in 96-well plates at 100,000 cells cm–2 density 1 week before imaging and staining. On the day of imaging, cells were gently washed with phenol-red-free RPMI (Thermo Fisher, 11835063), then stained with 1 μM Fluo-4 (Thermo Fisher, F14201) solution and incubated at 37 °C for 40 min. After incubation, the Fluo-4 solution was replaced with fresh, pre-warmed phenol-red-free RPMI. Then, hiPSC-CMs were transferred to a microscopy humidified chamber (pre-set at 37 °C with 5% CO2) and allowed to acclimatize for 10 min. Cells were imaged using a Zeiss Axio Observer inverted widefield microscope with a ×20/0.8 objective. The time series experiment was performed with 10 ms exposure time on the Fluo-4 channel (excitation 494 nm, emission 516 nm), and recorded for 10 s at 100 f.p.s. Stage positions for each well were stored to allow recurrence after drug treatment. All wells were scanned for baseline calcium transients. Then, hiPSC-CMs were treated with ivabradine (Sigma, SML0281), GLP-1 (Tocris, 2082) or the corresponding vehicles (DMSO for ivabradine and water for GLP-1). After 20 min, cells were scanned again using the same configuration.
Quantification of calcium imaging videos was performed using Fiji (v.2.1.0). To correct the intensity decay due to photobleaching, the Fiji bleach correction plugin (exponential fit method) was applied to the raw image series. Then, corrected image series were divided by their minimum intensity to remove background. A confluent region (containing at least 20 cells) was selected for quantification. Spiky plugin (https://github.com/PCCV/Spiky)91 was used to automate peak detection and quantification. Amplitude (peak value minus baseline value), Pk2Pk (time between two consecutive peaks), Time2Pk (time between threshold to the peak) and RW50 (time at 50% of amplitude from peak to next baseline at right) were averaged from multiple peaks detected. All results were normalized to the corresponding baseline values.
Statistical tests were performed in Prism 9 using D’Agostino and Pearson test to test for normal distribution. Unpaired t-test was used to compare vehicle and treated groups. P < 0.05 was regarded as significant.
Immunofluorescence staining of cells
hiPSC-CMs were fixed in 3.7% formalin in PBS then either blocked and permeabilized in 4% (w/v) BSA (Sigma Aldrich, A3059) and 0.2% (v/v) Triton X-100 (Thermo Fisher, 85111) in PBS (Gibco, 20012027) (for HCN1 and HCN4 staining) or blocked in 1% (w/v) BSA and 5% (v/v) normal goat serum (EDM Millipore, S26-100ML) in PBS (for GLP1R staining) for 30 min at room temperature. Incubation with primary antibodies (Supplementary Table 10) diluted in BSA, Triton X-100 and PBS (HCN1 and HCN4) or BSA, goat serum and PBS (GLP1R) was done overnight at 4 °C. Isotype controls and secondary antibody only (Supplementary Table 10) stainings were performed as negative controls. Following three washes with PBS, cells were stained with secondary antibodies diluted in BSA, Triton X-100 and PBS or BSA, goat serum and PBS was done for 1 h at room temperature (Supplementary Table 10). After three washes in PBS, cell nuclei were stained with DAPI (Invitrogen, D1306) for 15 min at room temperature. DAPI was rinsed and PBS was added to the cells. Confocal imaging acquisition was performed using a Zeiss LSM-780 inverted microscope with a EC Plan Neofluar ×40/1.3 oil objective at the Imperial College London Hammersmith FILM facility using 405 nm, 488 nm and 633 nm lasers for excitation. Image processing was performed in Fiji (v.2.1.0).
smFISH
The FFPE-embedded heart tissue sections (with a thickness of 5 μm) were placed onto SuperFrost Plus slides. Staining with a RNAscope Multiplex Fluorescent Reagent kit v2 assay (Advanced Cell Diagnostics, Bio-Techne) was automated using a Leica BOND RX, according to the manufacturer’s instructions. The tissues were baked and dewaxed on the Leica Bond RX, followed by the application of a heat-induced epitope retrieval step with epitope retrieval 2 for 15 min at 95 °C and protease digestion with protease III for 15 min. Subsequent processing included RNAscope probe hybridization and channel development with Opal 520, Opal 570 and Opal 650 dyes (Akoya Biosciences) at a concentration of 1:1,000, and streptavidin-conjugated Atto-425 (Bio Trend) at a concentration of 1:400 using TSA-biotin (TSA Plus Biotin Kit, Perkin Elmer). All nuclei were DAPI stained. Stained sections were imaged using a Perkin Elmer Opera Phenix High-Content Screening System with a ×20 water-immersion objective (NA of 0.16, 0.299 μm per pixel). The following channels were used: DAPI (excitation, 375 nm; emission, 435–480 nm); Opal 520 (excitation, 488 nm; emission, 500–550 nm); Opal 570 (excitation, 561 nm; emission, 570–630 nm); Opal 650 (excitation, 640 nm; emission, 650–760 nm); and Atto 425 (excitation, 425 nm; emission, 463–501 nm). Stained sections were also imaged on a Hamamatsu S60 with a ×40 objective (0.23 μm per pixel).
RNAscope quantification
Quantification of RNAscope images was performed using ImageJ. To remove background, each channel was initially normalized using the following steps: (1) subtracting the raw image with a Gaussian blur transformation with σ = 50; (2) performing a background subtraction with rolling ball radius of 50 pixels; and (3) setting every pixel with an intensity lower than 40 (of an 8-bit image) to 0. Following normalization, each section was quantified with sequential ROIs of 200 × 200 μm. To avoid any bias due to the placement of the initial ROI, quantification was performed over three rounds, with the initial ROI displaced by 50 μm in the x and y axis in each round. The following parameters were recorded per ROI: number of nuclei; COL1A1 area; and NPPB area.
For data analysis, to avoid variation due to cell density, for each ROI, COL1A1 and NPPB areas were normalized by dividing the area value by the number of nuclei in the ROI. Only ROIs with more than ten nuclei were considered for analysis. An ROI was only considered to contain a COL1A1 or NPPB niche if the staining area was equal to or higher than 1 μm2. To avoid any confounding effect due to different section sizes, we quantified the increase in the number of COL1A1 and NPPB niches through the average ROI percentage, consisting of the number of COL1A1 or NPPB niche ROI divided by the total number of ROI, averaged over the three rounds of quantitation. To quantify the expression of COL1A1 and NPPB, we averaged the normalized area of each ROI per niche over the three rounds of quantitation.
Immunofluorescence staining of tissues
The FFPE-embedded heart tissue samples were sectioned at 6 µm thickness and placed on VWR Superfrost Plus Microscope slides. Deparaffinization was performed in xylene (twice for 10 min), followed by graded washes in 100% ethanol (twice for 10 min), 95% ethanol for 5 min, 70% ethanol for 5 min, 50% ethanol for 5 min, and incubated in deionized water for rehydration. Antigen retrieval was performed using a proteinase K kit (Abcam ab64220) for 5 min at room temperature. Following antigen retrieval, sections were permeabilized and blocked in 0.1 M Tris containing 0.1% Triton X-100 (Sigma), 1% normal mouse serum, 1% normal goat serum and 1% BSA (R&D Systems). Samples were stained for 2 h at room temperature in a wet chamber with the appropriate antibodies, washed three times in PBS and mounted in Fluoromount-G (Southern Biotech). Images were acquired using a TCS SP8 (Leica) inverted confocal microscope with a ×40/1.1 NA water objective. Raw imaging data were processed using Imaris (Bitplane). The antibody information is provided in Supplementary Table 10.
For node and glial cell staining of fresh-frozen heart tissue, 10 µm sections were cut on to Superfrost plus slides and stored at −80 °C. These were thawed for 10 min at room temperature, briefly rehydrated in TBS and fixed with room temperature 4% paraformaldehyde for 5 min. Slides were then immersed for 10 min in TBS, and a hydrophobic pen was used to delimit the area of staining before starting the permeabilization process (0.25% saponin in TBS for 10 min). Blocking buffer (0.3 M glycine in antibody dilution buffer) was applied for 1 h at room temperature, followed by an overnight incubation with primary antibody at 4 °C in a humidified chamber. The primary antibodies used were mouse anti-HCN1 (Abcam, ab84816) at 1:100 dilution and rabbit anti-PLP1 (Abcam ab254363) at 1:500 in 10% normal goat serum in 0.2% Tween-20 and TBS. The isotype controls were Abcam ab37355 and ab172730. Slides were washed three times (5 min each) in 0.2% Tween-20 and TBS with gentle shaking, then secondary antibody solution was applied for 1 h at room temperature in the dark (1:1,000 goat anti-rabbit IgG AF555, LifeTech A21428 and 1:1,000 goat anti-mouse IgG AF647Plus, Fisher 15627898). Slides were washed three times (5 min each) in 0.2% Tween-20 and TBS, incubated with DAPI for 15 min at room temperature in the dark (Invitrogen D1306, 5 mg ml–1 stock then diluted 1:50,000), washed again briefly and mounted with ProLong gold antifade mountant (Thermo Fisher). Slides were scanned at ×40 magnification on a Hamamatsu NanoZoomer S60, and ROIs were imaged on a Leica SP8 confocal at ×20 magnification. The antibody information is provided in Supplementary Table 10.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.