Epistasis regulates genetic control of cardiac hypertrophy

Figures

• 

  Fig. 1. Schematic of the study workflow. The study workflow includes four major stages: (1) derivation of LV mass from cardiac MRIs (green boxes), (2) computational prioritization of epistatic drivers (orange), (3) functional interpretation of the hypothesized epistatic genetic loci (purple) and (4) experimental confirmation of epistasis in cardiac tissues and cells (blue). LVM, left ventricular mass; LVMi, LVM indexed by body surface area; PLINK26 and BOLT-LMM27, two different GWAS software packages; ANNOVAR74, a software for functional annotation of genetic variants; CADD36 scores the deleteriousness of variants; RegulomeDB35, a database that scores functional regulatory variants; ChromHMM34, a multivariate Hidden Markov Model for chromatin state annotation; PheWAS, phenome-wide association study; D40, cultured for 40 days post-differentiation from iPSCs.

• 

  Fig. 2. Lo-siRF workflow. a, Lo-siRF took as input SNV data and cardiac MRI-derived LV mass indexed by body surface area (LVMi) from 29,661 UKBB participants. b, Two GWAS studies were applied using different software to reduce the dimensionality of SNVs. c, LVMi was binarized into high and low categories using three thresholds (shown as stacked boxes). d, For each threshold, an siRF was fitted using the GWAS-filtered SNVs to predict the binarized LVMi, achieving prediction accuracy better than or comparable to other methods, before interpreting the model fit. PRS, polygenic risk score; ML, machine learning; DL, deep learning. e, siRF aggregates SNVs into genetic loci using ANNOVAR74 and finally ranks loci and their interactions based on a stability-driven importance score averaged across the three binarization thresholds.

• 

  Fig. 3. Lo-siRF finds epistasis between genetic risk loci for LV hypertrophy. a, Circos plot illustrating genetic risk loci identified by lo-siRF (green, circle 2) and regions after clumping FUMA-extracted SNVs in LD (r2 > 0.6) with lo-siRF-prioritized SNVs (black, circle 3). Circle 1 shows the top 300 epistatic SNV–SNV pairs ranked by occurrence frequency in lo-siRF (green), SNV–gene linkages (FDR < 0.5) from GTEx37 v8 cis-eQTL data of heart and skeletal muscle tissues (purple), and 3D chromatin interactions33 from Hi-C data of LV tissue (GSE87112). Circle 5 and 6 show occurrence frequencies and the number of partner SNVs in epistasis (normalized by the maximum value per locus) identified by lo-siRF, respectively. Circle 7 displays the ChromHMM34 core-15 chromatin state for the LV, right ventricle (RV), right atrium (RA) and fetal heart (Fetal), using the standard color scheme (bottom right)34. Circle 8 presents a GWAS Manhattan plot (points, two-sided P < 0.05 from PLINK26), in which SNVs are color coded by their maximum r2 value relative to the 283 lo-siRF-prioritized SNVs. The outer heatmap layer represents LD-linked SNVs (r2 > 0.6) extracted from the FUMA reference panel, which lack…

• 

  Fig. 4. Genes mapped from epistatic loci show strong functional co-associations. a, UpSet plot showing the number of lo-siRF-prioritized SNVs (dark blue) and their LD-linked (r2 > 0.6) SNVs (light blue, circle 8 in Fig. 3a) functionally mapped to 20 protein-coding genes via positional, eQTL and/or chromatin interaction (CI) mapping in FUMA33. CCDC141 and IGF1R (highlighted red) are prioritized by all three mapping strategies (details in Supplementary Data 4). b, Heatmap of averaged expression (GTEx, 50% winsorization, log2(TPM + 1)) across different tissues for these functionally mapped genes. c, Co-association network built from an enrichment analysis integrating multiple annotated gene set libraries for GO and pathway terms from Enrichr43. Nodes represent genes (green) functionally linked from lo-siRF-prioritized epistatic and hypostatic loci (Table 1) and top enriched GO and pathway terms. Edge width indicates the co-association strength between enriched terms and genes ranked by two-sided P values from an exhaustive permutation of the co-association score for all possible gene–gene combinations in the network (details in Methods and Supplementary Data 6). d, Comparison between the…

• 

  Fig. 5. Network analysis using transcriptomic data from 313 human hearts reveals strong correlations between suggestive statistical epistasis contributors. a, Gene co-expression networks for control (blue) and heart failure (red) conditions were established using weighted gene co-expression network analysis (WGCNA) on transcriptomic data from 313 non-failing and failing human heart tissues47. b,c, The connectivity between lo-siRF-prioritized genes was assessed against the full connectivity distributions of all possible gene–gene combinations in the control (b) and heart failure (c) networks. CCDC141 showed a significant connectivity to SYNM and LYSMD4 (IGF1R lo-siRF locus) and TTN and FKBP7 (TTN lo-siRF locus) in the control network. The color scales indicate two-sided empirical P values. d, Comparing between the control and heart failure networks indicated a significant connectivity decrease of these gene pairs during heart failure progression. e, A Sanky plot showing the rewired gene modular assignments for the lo-siRF loci-associated genes (middle column) in the control versus heart failure networks. Module names47 in the control (left column) and heart failure (right column) networks…

• 

  Fig. 6. Single-cell image analysis of epistatic gene-silencing effects on cardiomyocyte hypertrophy. a, hiPSC-CMs, with and without HCM (MYH7-R403Q mutation), were transfected with scramble siRNA or siRNAs specifically targeting single (CCDC141, IGF1R, TTN) or combined (CCDC141–IGF1R, CCDC141–TTN) loci prioritized by lo-siRF. RISC, RNA-induced silencing complex. b, Gene-silenced cardiomyocytes were bifurcated into large and small cells using a spiral microfluidic device (Extended Data Fig. 7) for high-resolution single-cell imaging. c, Image analysis workflow. Single-cell images were analyzed using a customized MATLAB program to extract cell size and shape features through a multistep process.

• 

  Fig. 7. CCDC141 nonadditively interacts with TTN and IGF1R to modify cardiomyocyte morphology. a, Violin plots showing total cell diameter distributions for unaffected (blue, n = 36,327) and MYH7-R403Q variant (red, n = 46,542) cardiomyocytes from three independent cell batches (P = 2.1 × 10−36, two-sided Mann–Whitney U test). In all panels, box plots show the interquartile range (IQR) with the median line, and whiskers extend to 1.5× IQR. Overlaid points represent median diameters from individual batches. b, Gene-silencing efficiency measured by RT-qPCR in unaffected (blue) and MYH7-R403Q variant (red) cells from at least n = 3 cell batches (exact n numbers provided in Source Data). The bars, overlaid points and error bars show mean values, individual biological replicates and standard deviations, respectively. c–h, Gene-silencing effects on single-cell morphology analyzed for n = 3 independent cell batches (n = 4 for unaffected cells silencing CCDC141), with over 33,000 cells measured per condition per batch. c, Percentage change in median cell diameter for each gene-silencing condition relative to scramble controls (that is, relative size difference). The squares represent mean…

• 

  Extended Data Fig. 1. Distribution of LVM and LVMi measurements for 29,661 UK Biobank participants. Left ventricular mass (LVM, a) and LVM indexed to body surface area (LVMi, b) measurements were extracted from cardiac magnetic resonance imaging for 29,661 unrelated White British individuals via deep learning20. c, A high Pearson correlation of 0.92 was observed between these LVM and LVMi measurements. Source data

• 

  Extended Data Fig. 2. LVMi GWAS using BOLT-LMM and PLINK. Two-sided GWAS p-values from BOLT-LMM (a) and PLINK (b) identified LVMi-associated SNVs, of which the top lead SNV rs3045696 showed the highest significance in both analyses. Other labeled lead SNVs were significant in only one of the two methods. The red dashed line indicates the genome-wide significance threshold (p < 5E-8, two-sided). Notably, rs3045696 and rs67172995 were also stably prioritized by lo-siRF as epistasis interactor variants. Further details are provided in Supplementary Data 1. Source data

• 

  Extended Data Fig. 3. Differences in local stability scores between high and low LV mass highlight the importance of the lo-siRF-prioritized interactions between genetic loci. a, Schematic of local stability importance score computation. Given a locus (light blue transcript), the local stability importance score for an individual is defined as the proportion of trees for which at least one SNV (shaded black region) in the locus is used in the individual’s decision path. This computation (top) was performed for each individual (denoted by the stacked boxes). Then, a two-sided permutation test was conducted to assess the difference in these local stability importance scores between the low and high LVMi individuals (bottom). b, Differences in the distribution of local stability importance scores suggest that the identified interactions between genetic loci are important for differentiating individuals with high (dark gray) and low (orange) LVMi in the siRF fit. This result, evaluated on the validation data, is stable across the three binarization thresholds and is quantified by a nominal two-sided p-value from a permutation test given in the top right corner of each subplot. Given the use of 10,000…

• 

  Extended Data Fig. 4. Top SNVs from lo-siRF-prioritized loci and interactions between loci. The most important SNVs and SNV-SNV pairs, as measured by their proportion of occurrence in the siRF fit, are annotated for the top lo-siRF-prioritized interactions between loci in a and top genetic loci in b. The y-axis shows the proportion of decision paths in siRF, for which the SNV or SNV-SNV pair occurs, averaged across all three binarization thresholds. In each of the interactions between genetic loci, rs7591091 in the CCDC141 locus appears most frequently, suggesting a key role in cardiac hypertrophy. Source data

• 

  Extended Data Fig. 5. Genes mapped from epistatic loci share transcription factors and splicing regulators. Gene pairs exhibiting strong co-association with transcription factors (TFs) and RNA-binding regulators (p < 0.05, two-sided Fisher’s exact test, Fig. 4e) are displayed. Horizontal bars indicate the number of shared TFs or RNA-binding regulators (bottom axis), with the top-ranked TF or regulator, determined by the lowest two-sided enrichment p-value (dots, top axis), labeled on each bar. Bars are color-coded by the corresponding gene set library. Detailed information can be found in Supplementary Data 7.