Genetic architecture of natural variation of cardiac performance from flies to humans

  1. Saswati Saha
  2. Lionel Spinelli
  3. Jaime A Castro Mondragon
  4. Anaïs Kervadec
  5. Michaela Lynott
  6. Laurent Kremmer
  7. Laurence Roder
  8. Sallouha Krifa
  9. Magali Torres
  10. Christine Brun
  11. Georg Vogler
  12. Rolf Bodmer
  13. Alexandre R Colas  Is a corresponding author
  14. Karen Ocorr  Is a corresponding author
  15. Laurent Perrin  Is a corresponding author
  1. Aix-Marseille University, INSERM, TAGC, Turing Center for Living systems, France
  2. Centre for Molecular Medicine Norway (NCMM), Norway
  3. Development, Aging and Regeneration Program, Sanford Burnham Prebys Medical Discovery Institute, United States
  4. CNRS, France

Abstract

Deciphering the genetic architecture of human cardiac disorders is of fundamental importance but their underlying complexity is a major hurdle. We investigated the natural variation of cardiac performance in the sequenced inbred lines of the Drosophila Genetic Reference Panel (DGRP). Genome-wide associations studies (GWAS) identified genetic networks associated with natural variation of cardiac traits which were used to gain insights as to the molecular and cellular processes affected. Non-coding variants that we identified were used to map potential regulatory non-coding regions, which in turn were employed to predict transcription factors (TFs) binding sites. Cognate TFs, many of which themselves bear polymorphisms associated with variations of cardiac performance, were also validated by heart-specific knockdown. Additionally, we showed that the natural variations associated with variability in cardiac performance affect a set of genes overlapping those associated with average traits but through different variants in the same genes. Furthermore, we showed that phenotypic variability was also associated with natural variation of gene regulatory networks. More importantly, we documented correlations between genes associated with cardiac phenotypes in both flies and humans, which supports a conserved genetic architecture regulating adult cardiac function from arthropods to mammals. Specifically, roles for PAX9 and EGR2 in the regulation of the cardiac rhythm were established in both models, illustrating that the characteristics of natural variations in cardiac function identified in Drosophila can accelerate discovery in humans.

Editor's evaluation

The authors investigated natural variation and new genetic mechanisms underlying cardiac performance using sequenced inbred lines of the Drosophila Genetic Reference Panel. The study provides insights into the genetic architecture of complex cardiac performance traits and represents an important resource for researchers studying cardiac performance.

https://doi.org/10.7554/eLife.82459.sa0

Introduction

Heart diseases is a major cause of mortality (Bezzina et al., 2015). Although a large number of genome-wide association studies (GWAS) have identified hundreds of genetic variants related to cardiovascular traits (Roselli et al., 2018; van Setten et al., 2018; Shah et al., 2020; Verweij et al., 2020), we are very far from a comprehensive understanding of the genetic architecture of these complex traits. Deciphering the impact of genetic variations on quantitative traits is however critical for the prediction of disease risk. But disentangling the relative genetic and environmental contributions to pathologies is challenging due to the difficulty in accounting for environmental influences and disease comorbidities. Underlying epistatic interactions may also contribute to problems with replication in human GWAS performed in distinct populations which rarely take epistatic effects into account. In addition, linking a trait associated locus to a candidate gene or a set of genes for prioritization is not straightforward (Mackay, 2014, Boyle et al., 2017). Furthermore, the analysis of genetic factors related to cardiac traits is complicated by their interactions with several risk factors, such as increasing age, hypertension, diabetes mellitus, ischemic, and structural heart disease (Paludan-Müller et al., 2016).

These pitfalls can be overcome using animal models. Model organisms allow precise controlling of the genetic background and environmental rearing conditions. They can provide generally applicable insights into the genetic underpinnings of complex traits and human diseases, due to the evolutionary conservation of biological pathways. Numerous studies have highlighted the conservation of cardiac development and function from flies to mammals. Indeed, orthologous genes control the early development as well as the essential functional elements of the heart. The fly is the simplest genetic model with a heart muscle and is increasingly used to identify the genes involved in heart disease and aging (Ocorr et al., 2007b; Diop and Bodmer, 2015; Rosenthal et al., 2010). Although a large number of genes are implicated in establishing and maintaining cardiac function in Drosophila (Neely et al., 2010), the extent to which genes identified from mutant analysis reflect naturally occurring variants is neither known, nor do we know how allelic variants at several segregating loci combine to affect cardiac performance.

We previously showed that wild populations of flies harbor rare polymorphisms of major effects that predispose them to cardiac dysfunction (Ocorr et al., 2007a). Here, we analyzed the genetic architecture of the natural variation of cardiac performance in Drosophila. Our aims were to (i) identify the variants associated with cardiac traits found in a natural population, (ii) decipher how these variants interact with each other and with the environment to impact cardiac performance, and (iii) gain insights into the molecular and cellular processes affected. For this, we used the Drosophila Genetic Reference Panel (DGRP) (Mackay et al., 2012; Huang et al., 2014), a community resource of sequenced inbred lines. Previous GWAS performed in the DGRP indicate that inheritance of most quantitative traits in Drosophila is complex, involving many genes with small additive effects as well as epistatic interactions (Mackay and Huang, 2018). The use of inbred lines allows us to assess the effects of genetic variations in distinct but constant genetic backgrounds and discriminate genetic and environmental effects.

We demonstrated substantial among-lines variations of cardiac performance and identified genetic variants associated with the cardiac traits together with epistatic interactions among polymorphisms. Candidate loci were enriched for genes encoding transcription factors (TFs) and signaling pathways, which we validated in vivo. We used non-coding variants - which represented the vast majority of identified polymorphisms – for predicting transcriptional regulators of associated genes. Corresponding TFs were further validated in vivo by heart-specific RNAi-mediated knockdown (KD). This illustrates that natural variations of gene regulatory networks have widespread impact on cardiac function. In addition, we analyzed the phenotypic variability of cardiac traits between individuals within each of the DGRP lines (i.e., with the same genotype), and we documented significant diversity in phenotypic variability among the DGRP lines, suggesting genetic variations influenced phenotypic variability of cardiac performance. Genetic variants associated with this phenotypic variability were identified and shown to affect a set of genes that overlapped with those associated with trait means, although through different genetic variants in the same genes.

Comparison of human GWAS of cardiac disorders with results in flies identified a set of orthologous genes associated with cardiac traits both in Drosophila and in humans, supporting the conservation of the genetic architecture of cardiac performance traits, from arthropods to mammals. siRNA-mediated gene KD were performed in human induced pluripotent stem cells derived cardiomyocytes (hiPSC-CMs) to indeed show that dmPox-meso/hPAX9 and dmStripe/hEGR2 have conserved functions in cardiomyocytes from both flies and humans. These new insights into the fly’s genetic architecture and the connections between natural variations and cardiac performance permit the accelerated identification of essential cardiac genes and pathways in humans.

Results

Quantitative genetics of heart performances in the DGRP

In this study, we aimed to evaluate how naturally occurring genetic variations affect cardiac performance in young Drosophila adults and identify variants and genes involved in the genetic architecture of cardiac traits. To assess the magnitude of naturally occurring variations of the traits, we measured heart parameters in 1-week-old females for 167 lines from the DGRP, a publicly available population of sequenced inbred lines derived from wild caught females (Figure 1A). Briefly, semi-intact preparations of individual flies (Ocorr et al., 2007c) were used for high-speed video recording combined with Semi-automated Heartbeat Analysis (SOHA) software (http://www.sohasoftware.com/) which allows precise quantification of a number of cardiac function parameters (Fink et al., 2009; Cammarato et al., 2015). Fly cardiac function parameters are highly influenced by sex (Wessells et al., 2004). Due to the experimental burden of analyzing individual cardiac phenotypes, we focused on female flies only and designed our experiment in the following way: we randomly selected 14 lines out of 167 that were replicated twice. The remaining 153 lines were replicated once. Each replicate was composed of 12 individuals. No block effect was observed due to the replicates in the 14 selected lines (see Supplementary file 1a). This allowed us to perform our final analysis on one replicate of each of the 167 lines. A total sample of 1956 individuals was analyzed. Seven cardiac traits were analyzed across the whole population (Figure 1—source data 1 and Table 1). As illustrated in Figure 1A, we analyzed phenotypes related to the rhythmicity of cardiac function: the systolic interval (SI) is the time elapsed between the beginning and the end of one contraction, the diastolic interval (DI) is the time elapsed between two contractions and the heart period (HP) is the duration of a total cycle (contraction+relaxation (DI+SI)). The arrhythmia index (AI, std-dev(HP)/mean (HP)) is used to evaluate the variability of the cardiac rhythm. In addition, three traits related to contractility were measured. The diameters of the heart in diastole (end diastolic diameter [EDD]), in systole (end systolic diameter [ESD]), and the fractional shortening (FS), which measures the contraction efficacy (EDD −ESD/EDD). We found significant genetic variation for all traits (Figure 1B and Figure 1—figure supplement 1) with broad sense heritability ranging from 0.30 (AI) to 0.56 (EDD) (Table 1). Except for EDD/ESD and HP/DI, quantitative traits were poorly correlated with each other (Figure 1—figure supplement 1).

Figure 1 with 2 supplements see all
Quantitative genetics and genome-wide associations studies (GWAS) for cardiac traits in the Drosophila Genetic Reference Panel (DGRP).

(A) Left: Cardiac performance traits were analyzed in 167 sequenced inbred lines from the DGRP population. Approximately 12 females per line were analyzed. Right panels: Schematic of the Drosophila adult heart assay and example of M-mode generated from video recording of a beating fly heart. Semi-intact preparations of 1-week-old adult females were used for high-speed video recording followed by automated and quantitative assessment of heart size and function. The representative M-mode trace illustrate the cardiac traits analyzed. DI: diastolic interval; SI: systolic interval; HP: heart period (duration of one heartbeat); EDD: end diastolic diameter (fully relaxed cardiac tube diameter); ESD: end systolic diameter (fully contracted cardiac tube diameter). Fractional shortening (FS=EDD − ESD/EDD) and arrythmia index (AI=Std Dev (HP)/HP) were additionally calculated and analyzed. (B) Distribution of line means and within lines variations (box plots) from 167 measured DGRP lines for HP and EDD. DGRP lines are ranked by their increasing mean phenotypic values. For both phenotypes, representative M-modes from extreme lines are shown below (other traits are displayed in Figure 1—figure supplement 1). (C) Pearson residuals of chi-square test from the comparison of indicated single nucleotide polymorphism (SNP) categories in the DGRP and among variants associated with cardiac traits. According to DGRP annotations, SNPs are attributed to genes if they are within the gene transcription unit (5’ and 3’ UTR, synonymous and non-synonymous coding, introns) or within 1 kb from transcription start and end sites (1 kb upstream, 1 kb downstream). NA: SNPs not attributed to genes (>1 kb from transcription start site [TSS] and transcription end sites [TES]). (D) Comparison of gene sets identified by single marker using Fast-LMM (LMM) and in interaction using FastEpistasis (Epistasis). The Venn diagram illustrates the size of the two populations and their overlap. (E) Overlap coefficient of gene sets associated with the different cardiac traits analyzed.

Figure 1—source data 1

Individual values for cardiac traits analyzed across the 167 Drosophila Genetic Reference Panel (DGRP) lines.

Individual and DGRP line number are indicated. Phenotypic values were determined from high-speed video recording on dissected flies and movie analysis using Semi-automated Heartbeat Analysis (SOHA) (Mackay et al., 2012).

https://cdn.elifesciences.org/articles/82459/elife-82459-fig1-data1-v1.xlsx
Figure 1—source data 2

Variants identified by FastLMM as associated to indicated phenotypes.

Among the 100 best ranked associations, only variants with MAF >4% were retained. Tables for variants mapped to genes and for variants that are not within gene mapping criteria (>1 kb from transcription start site [TSS] and transcription end sites [TES]) are indicated.

https://cdn.elifesciences.org/articles/82459/elife-82459-fig1-data2-v1.xlsx
Figure 1—source data 3

All FastEpistasis data on mean phenotypes, per quantitative trait.

Single nucleotide polymorphism (SNP) ID, position, associated genes, and statistics are indicated for both focal SNPs (left) and their interacting SNPs (right). Each sheet displays the results for indicated quantitative traits, except for the first one which is a merge of all quantitative traits association analyses.

https://cdn.elifesciences.org/articles/82459/elife-82459-fig1-data3-v1.xlsx
Table 1
Quantitative genetics of cardiac traits in the Drosophila Genetic Reference Panel (DGRP).

Summary statistics over all DGRP genotypes assayed. Number of lines and individuals (after outlier removal, see Materials and methods) analyzed for each cardiac trait is indicated. Mean, standard deviation (Std dev.), and coefficient of variation (Coef. Var) among the whole population are indicated. Genetic, environment, and phenotypic variance (respectively Genet. var, Env. var, and Phen. var) were calculated for each trait. Broad sense heritability of traits means (H2) suggested heritability of corresponding traits. Levene test indicated significant heterogeneity of the variance among the lines.

DiastolicintervalsSystolicintervalsHeartperiodDiastolic DiameterSystolic diameterFractional shorteningArrhythmia Index
total.nb.lines167167167167167167167
mean0.46380.21660.688379.420051.05000.35380.2475
Std dev.0.263300.032160.2769014.090009.493000.068370.29230
Coef. var0.56770.14850.40220.17740.18600.19331.1810
lines (mean)165166165159157158166
Indiv. (mean)1914191119201779175317671832
lines (Cve)165166165159157158166
Indiv. (Cve)1914191119201779175317671832
Genet. var2.59e-025.03e-042.87e-021.13e+024.39e+011.57e-032.21e-02
Env. var4.36e-025.35e-044.82e-028.64e+014.65e+013.11e-036.35e-02
Phen. var6.95e-021.04e-037.68e-021.99e+029.04e+014.68e-038.56e-02
H20.3730.4850.3730.5660.4850.3350.258
F value76,86411,68674,71546,95015,04111,16465,308
Pr(F)8.8e-1202.3e-1875.8e-1167.1e-621.9e-2318.8e-1751.8e-96
Levene test1.9e-101.9e-101.7e-081.6e-052.1e-131.6e-052.1e-13

GWAS analyses of heart performance

To identify candidate variants associated with cardiac performance variation, we performed GWAS analyses and evaluated single marker associations of line means with common variants using a linear mixed model (Lippert et al., 2011) and after accounting for effects of Wolbachia infection and common polymorphism inversions (see Materials and methods). Genotype-phenotype associations were performed separately for all seven quantitative traits and variants were ranked based on their p-values. For most of the phenotypes analyzed, quantile-quantile (QQ) plots were uniform (Figure 1—figure supplement 2) and none of the variants reached the strict Bonferroni correction threshold for multiple tests (2 · 10–8), which is usual in the DGRP given the size of the population. However, the decisive advantage of the Drosophila system is that we can use GWA analyses as primary screens for candidate genes and mechanisms that can be subsequently validated by different means. We therefore chose to analyze the 100 top ranked variants for each quantitative trait. This choice is based on our strategy to test the selected single nucleotide polymorphisms (SNPs) and associated genes by a variety of approaches – data mining and experimental validation (see below) – in order to provide a global validation of association results and to gain insights into the characteristics of the genetic architecture of the cardiac traits. This cut-off was chosen in order to be able to test a significant number of variants while being globally similar to the nominal cut-off (10–5) generally used in DGRP analyses. A large proportion of the variants retained have indeed a p-value below 10–5. Selected variants were further filtered on the basis of minor allele frequency (MAF >4%) (Figure 1—source data 2, Supplementary file 1b). Among the seven quantitative traits analyzed, we identified 530 unique variants. These variants were associated to genes if they were within 1 kb of transcription start site (TSS) or transcription end sites (TES). Using these criteria, 417 variants were mapped to 332 genes (Supplementary file 1c). We performed a chi-squared test to determine if the genomic location of variants associated with cardiac traits is biased toward any particular genomic region when compared with the whole set of variants with MAF >4% in the DGRP population and obtained a p-value of 2.778e-13. Genomic locations of the variants were biased toward regions within 1 kb upstream of genes TSS, and, to a lesser extent, to genes 5’ UTR (Figure 1C). Variants not mapped to genes (located at >1 kb from TSS or TES) were slightly depleted in our set.

In GWAS analyses, loci associated with a complex trait collectively account for only a small proportion of the observed genetic variation (Manolio et al., 2009) and part of this ‘missing heritability’ is thought to come from interactions between variants (Flint and Mackay, 2009; Manolio et al., 2009; Huang et al., 2012; Shorter et al., 2015). As a first step toward identifying such interactions, we used FastEpistasis (Schüpbach et al., 2010). SNP identified by GWAS were used as focal SNPs and were tested for interactions with all other SNPs in the DGRP. FastEpistasis reports best ranked interacting SNP for each starting focal SNP, thus extending the network of variants and genes associated to natural variation of cardiac performance, which were used for hypothesis generation and functional validations; 288 unique SNPs were identified, which were mapped to 261 genes (Figure 1—source data 3, Supplementary file 1e). While none of the focal SNPs interacted with each other, there is a significant overlap between the 332 genes associated with single marker GWAS and the 261 genes identified by epistasis (n=31, Figure 1D and Supplementary file 1e, fold change (FC)=6; hypergeometric pval=6.8 × 10–16). This illustrates that the genes that contribute to quantitative variations in cardiac performance have a tendency to interact with each other, although through distinct alleles.

Taken together, single marker GWAS and epistatic interactions performed on the seven cardiac phenotypes identified a compendium of 562 genes associated with natural variations of heart performance (Supplementary file 1f). In line with the correlation noted between their phenotypes (Figure 1—figure supplement 1B), the GWAS for HP and DI identified partially overlapping gene sets (overlap index 0.23, Figure 1E). The same was true, to a lesser extent, for ESD and EDD (0.15). Otherwise, the sets of genes associated with each of the cardiac traits are poorly correlated with each other.

Functional annotations and network analyses of association results

Our next objective was to identify the biological processes potentially affected by natural variation in cardiac performance. Gene Ontology (GO) enrichment analysis of the combined single marker GWAS and epistatic interactions analyses indicated that genes encoding signaling receptors, TFs, and cell adhesion molecules were over-represented among these gene sets (pval=1.4 × 10–9 [FC=2.9], 5×10–4 [FC=2], and 3×10–3 [FC=4.6], respectively). There was also a bias for genes encoding proteins located at the plasma membrane, at ion channel complexes as well as components of contractile fibers (pval=3.4 × 10–10 [FC=3], 7×10–5 [FC=4.2], and 4×10–2 [FC=3.6]; Figure 2A; Supplementary file 2a). Of note, although a number of genes have previously been identified as being required during heart development or for the establishment and maintenance of cardiac function by single gene approaches, we found no enrichment for these gene categories in our analysis. In addition, genes identified in a global screen for stress-induced lethality following heart-specific RNAi KD (Neely et al., 2010) were also not enriched in GWAS detected genes (FC=1; Supplementary file 2b). This indicates that genes associated to natural variations of cardiac traits are typically missed by traditional forward or reverse genetics approaches, which highlights the value of our approach.

Figure 2 with 2 supplements see all
Functional annotations and validations of genes associated with genome-wide associations studies (GWAS) for cardiac performance.

(A) Gene Ontology (GO) enrichment analyses. Selected molecular functions (MF, left) and cellular components (CC, right) associated with cardiac performances at FDR < 0.05 are shown. Enrichment analysis was performed using G:profiler with a correction for multitesting (see Materials and methods). (B) Interaction network of genes associated with natural variations of cardiac performance. Direct genetic and physical interactions between cardiac fly GWAS genes are displayed. Nodes represent genes and/or proteins, edges represent interactions (blue: genetic; black: physical). Node shapes refer to single marker and/or epistasis GWAS, node color to the cardiac performances phenotype(s) for which associations were established. Genes and proteins highlighted in pink point to transcription factors, in green and blue to signaling pathways (FGF and TGFb, respectively), and in yellow to ion channels. (C) Heatmap representing the effects on indicated cardiac traits of heart-specific RNAi-mediated knockdown (KD) of 42 genes identified in GWAS for cardiac performance. Results of Wilcoxon rank sum test of the effects of indicated heart-specific RNAi-mediated gene KD (rows) for cardiac performance traits (columns), analyzed on semi-intact 1-week females flies. Detailed data are presented in Figure 2—source data 2. Thirty-eight (out of 42) genes tested lead to significant effects on cardiac performance traits upon KD. Black dots indicate the trait(s) for which the corresponding gene was associated in GWAS. ns: not significant; *: pval <0.05; **: pval <0.01; ***: pval <0.005; ****pval <0.0001 (p-values were adjusted for multiple testing using Bonferroni correction). Comparison with heart-specific effect of random selected genes is displayed in Figure 2—figure supplement 1. (D) Schematic drawing of BMP and activin pathways in Drosophila. (E) Genetic interactions between BMP and activin pathway genes. Genetic interactions tested between snooBSC234 and sogU2 for SI and between dppd14 and babo32 for FS (other phenotypes are shown in Figure 2—figure supplement 2). Cardiac traits were measured on each single heterozygotes and on double heterozygotes flies. Two-way ANOVA reveals that the interaction between snooBSC234 and sogU2 for SI and between dppd14 and babo32 for FS are significant. Detailed data for interaction effect corresponding to all phenotypes are displayed in Figure 2—figure supplement 2.

Figure 2—source data 1

Collection of physical (IP, Y2H) and genetic interactions identified in Drosophila.

https://cdn.elifesciences.org/articles/82459/elife-82459-fig2-data1-v1.xlsx
Figure 2—source data 2

Data from validation experiments.

(Heart-specific RNAi validations and tests for genetic interactions among BMP members.)

https://cdn.elifesciences.org/articles/82459/elife-82459-fig2-data2-v1.xlsx

In order to gain additional knowledge about the cellular and molecular pathways affected by natural variations of cardiac traits, we have mapped the associated genes and gene products onto characterized interaction networks. Of the 562 identified genes, 419 were mapped to the fly interactome that includes both physical (protein-protein) and genetic interactions from both DroID (Murali et al., 2011) and flybi databases (see Materials and methods and Figure 2—source data 1). Remarkably, a high proportion (148) of the GWAS identified genes were directly connected within the fly interactome and formed a large network of interacting genes/proteins (Supplementary file 2c and d, Figure 2B), suggesting that they participate in common biological processes. This network encompasses several TFs and ion channel complex genes, consistent with their potential role in the genetic architecture of natural variation of heart performance. Several components of signaling pathways are also present in the network, including members of the FGF and TGFβ pathways (see below).

Functional validations of candidate genes

To assess in an extensive way whether mutations in genes harboring SNPs associated with variation in cardiac traits contributed to these phenotypes, we selected 42 GWAS associated genes for cardiac-specific RNAi KD and tested the effects on cardiac performance. We selected genes that were identified in at least two independent GWAS for two traits or that were known to be dynamically expressed in the adult heart (Monnier et al., 2012) and for which inducible RNAi lines were available. Genes were tested in 1-week-old adult female flies, using the heart-specific Hand>Gal4 driver (Popichenko et al., 2007) and the same semi-intact heart preparations and SOHA analyses as for DGRP lines screening. Notably, 38 of the 42 selected genes led to various levels of cardiac performance defects following heart-specific KD (Figure 2C). In parallel, we tested the effect on cardiac performance of knocking down 18 genes randomly selected in the genome – the GWAS associated genes being excluded from the selection (see Materials and methods and Figure 2—figure supplement 1). Although a number of these genes lead to cardiac phenotypes when inactivated – which is consistent with published observations that quantitative traits can be influenced by a large number of genes (Zhang et al., 2021) – when inactivated in the heart, the genes selected from GWAS lead to significantly more frequent phenotypes compared to the randomly selected genes (Figure 2—figure supplement 1). These results therefore supported our association results. It is important to emphasize that our approach is limited to testing the effects of tissue-specific gene KD. Since some of the variants may lead to increased gene function and/or expression, this can lead to a false negative rate that is difficult to estimate. In addition, some of the associated variants may influence heart function by non-cell-autonomous mechanisms, which would not be replicated by cardiac-specific RNAi KD.

We further focused on the TGFβ pathway, since members of both BMP and activin pathways were identified in our analyses. We tested different members of the TGFβ pathway for cardiac phenotypes using cardiac-specific RNAi KD (Figure 2C), and confirmed the involvement of the activin agonist snoo (Ski orthologue) and the BMP antagonist sog (chordin orthologue). Notably, activin and BMP pathways are usually antagonistic (Figure 2D). Their joint identification in our GWAS suggest that they act in a coordinated fashion to regulate heart function. Alternatively, it may simply reflect their involvement in different aspects of cardiac development and/or functional maturation. In order to discriminate between these two hypotheses, we tested if different components of these pathways interacted genetically. Single heterozygotes for loss of function alleles show dosage-dependent effects of snoo and sog on several phenotypes, providing an independent confirmation of their involvement in several cardiac traits (Figure 2—figure supplement 2). Importantly, compared to each single heterozygotes, snooBSC234/sogUU2 double heterozygotes flies showed non-additive SI phenotypes (two-way ANOVA pval: 2.1 · 10–7), suggesting a genetic interaction (Figure 2E and Figure 2—figure supplement 2). It is worth noting however that snoo is also a transcriptional repressor of the BMP pathway (Takaesu et al., 2006). The effect observed in snooBSC234/sogUU2 double heterozygotes can therefore alternatively arise as a consequence of an increased BMP signaling without affecting the activin pathway. We thus tested other allelic combinations for loss of function alleles of BMP and activin pathways. babo/tkv heterozygotes (respectively activin and BMP type 1 receptors) displayed non-additive ESD and EDD phenotypes (Figure 2—figure supplement 2). Synergistic interaction of BMP and activin pathways was also suggested by the analysis of FS in loss of function mutants for babo and dpp, the BMP ligand (Figure 2—figure supplement 2). Of note, babo/tkv double heterozygotes also displayed a tendency to non-additive effects on FS albeit non-significant (two-way ANOVA p=0.054). In addition, mad/smox heterozygotes (specifc downstream TFs of BMP and activin pathways) displayed non-additive effects on several traits, including phenotypes related to rhythmicity (HP, SI, DI) and contractility (ESD and EDD) (Figure 2—figure supplement 2). Altogether, cardiac performance in response to allelic combinations of activin and BMP supported a coordinated action of both pathways in the establishment and/or maintenance of cardiac activity. This was further supported by the observation that simple heterozygotes for the tested loss of function alleles displayed similar trends with respect to cardiac performance, irrespective of the pathway considered (dpp, tkv, and mad for BMP; babo and smox for activin). Indeed, they displayed either no effect or increased FS and rhythmicity phenotypes (HP, DI, SI, AI), and decreased cardiac diameters (ESD and EDD). This suggests coordinated activity of both pathways. Importantly, the genetic interactions were tested using amorphic alleles that lead to systemic loss of function. The observed phenotypes may thus not unravel cardiac-specific effects of the pathways. In support of this, mad cardiac-specific RNAi KD was tested (see below, Figure 3) and lead to a decreased HP, DI, SI, and FS whereas heterozygotes for mad Neely et al., 2010 have either no (FS) or opposite (HP, DI, SI) effect on these phenotypes (Figure 2—figure supplement 2). Inversely, mad RNAi caused a significant increase in AI whereas mad (Neely et al., 2010) had no effect. However, heart-specific dpp RNAi KD (Figure 2—figure supplement 2) lead to similar phenotypic trends compared to dppd14 (increased HP, DI, SI, decreased EDD and ESD) with the notable exception of FS which was reduced following cardiac-specific KD (Figure 2—figure supplement 2), but increased in dppd14 heterozygotes (Figure 2—figure supplement 2). Taken together, these data point to a complex picture of TGFβ pathway activity in regulating cardiac performance, involving both the activin and the BMP pathways as well as gene-specific effects with both systemic and tissue-specific contributions.

Transcription factor binding site (TFBS) predictions from sequences surrounding non-coding variants and in vivo validation of cognate transcription factors (TFs).

(A) Schematic representation of the workflow used for motifs analysis in the vicinity of the non-coding variants and comparison with TFBSs databases. (B) Logo representation of position-specific scoring matrices (PSSMs) identified in corresponding sets of non-coding variants. Potential cognate TFs are indicated and were identified by comparing PSSMs to TFBSs databases. (C) Cardiac-specific RNAi knockdown phenotypes of TFs tested in vivo. Hand>Gal4; UAS RNAi 1-week adult females (n>30 per cross) were analyzed for cardiac performance in intact flies. Scatter plots of individual data, mean, and SD are represented. Wilcoxon test results are indicated (*: pval <0.05; **: pval <0.01; ***: pval <0.005; ****pval <0.0001). Data were normalized to genotype matched controls. Raw data are presented in Figure 2—source data 2. Detailed statistical analyses are in Supplementary file 2e.

Figure 3—source data 1

Variants in non-coding regions used for motif discovery, associated with end diastolic diameter (EDD), end systolic diameter (ESD), and fractional shortening (FS) (‘structure’) and to AI, DI, SI, and HP (‘rhythm’) for both mean and CVe of cardiac traits.

https://cdn.elifesciences.org/articles/82459/elife-82459-fig3-data1-v1.xlsx

As a whole, the cardiac phenotypes observed following RNAi-mediated gene KD of GWAS selected genes and the complex interplay characterized between different activin and BMP pathways members are strong validations of the GWAS results that are further supported by the functional analysis of non-coding variants and by the correlations we observed between identified fly genes and human genes associated to cardiac dysfunction (discussed below).

Natural variations of heart performance primarily affect gene regulatory networks

The over-representation of genes coding for TFs found in GWAS for cardiac performance traits (pval=5559 · 10–4, Figure 2A, Supplementary file 2a) suggested that perturbations of gene regulatory networks (GRNs) play a central role in the genetic architecture of natural variations of cardiac performance. In addition, most variants (95%) were localized in non-coding regions of the genome. Remarkably, their distribution was enriched in gene proximal (0–1 kb) upstream regions (Figure 1C) suggesting that a high proportion of variants affects cis-regulatory sequences.

Thus, we sought to determine whether potential TF binding motifs are enriched in the vicinity of variants identified by GWAS. Analysis of the sites of the variants themselves did not reveal any particular pattern. However, oligo analysis (van Helden et al., 1998) of the 150 bp sequences around these non-coding variants detected significant over-represented k-mers (Evalue ≤ 1e-4). These were further assembled and converted into positions-specific scoring matrices (PSSMs) and compared with PSSM databases of known fly TFs binding motifs were used for annotations (see Materials and methods, Figure 3—source data 1 and Figure 3A). Variants associated with quantitative traits linked to heart size and contractility (ESD, EDD, and FS; 367 variants) and those related to heart rhythmicity (AI, HP, DI, SI; 435 variants) were analyzed separately. We identified enrichment of potential DNA binding motifs for both sets of traits (Figure 3B). Although this approach predicted motifs that were specific for each group of variants, a significant proportion of potential DNA binding motifs were predicted by both groups, suggesting commonalities in the genetic architecture of the corresponding traits.

Remarkably, we identified DNA binding motifs for TFs that were themselves associated with natural variations of the cardiac traits. Some examples are klumpfuss (klu/ZBTB7A), luna (KLF5-7), longitudinals lacking (lola/BTB18), serpent (srp/GATA1-6), mirror (mirr/IRX3-6), Helix Loop Helix protein 3B (HLH3B/TAL2), and defective proventiculus (dve/SATB1) which were also identified by GWAS (see Supplementary file 1f). Potential binding sites for mothers against dpp, Mad (SMAD) were enriched in the vicinity of the non-coding variants associated with heart size and contractility. SMAD is the downstream effector of the BMP signaling pathway that we showed was also involved in regulating cardiac contractility (Figure 2D–E). Motifs predicted to be bound by the NK TFs Tinman (Tin/NKX2.5) and Bagpipe (Bap/NKX3.2) are over-represented in both sets of variants. tin has a well-characterized function in the cardiac GRN, both during embryogenesis (Bodmer, 1993) and at adulthood (Qian et al., 2011) and bap is also a member of the cardiac network (Seyres et al., 2016). Lastly, we found Jun-related antigen (Jra/JUN) TF binding sites (TFBS) enriched in sequences associated with natural variations of cardiac rhythm. We previously reported that djun KD leads to altered HP and AI in young flies (Monnier et al., 2012), supporting a function for this stress-activated TF in the cardiac GRN.

We did not test individually the effects on cardiac performance of mutations in predicted TFBSs located near the SNPs because any individual effect would probably be too small to be detectable by the available methods. Rather, we tested the potential involvement of their cognate TFs by cardiac-specific RNAi-mediated KD. luna, lola, and srp KD had already been analyzed in semi-intact preparations and shown to induce strong cardiac phenotypes (Figure 2C). We repeated these heart-specific manipulations employing additional, independent RNAi lines targeting these genes and an in vivo approach that uses the red fluorophore TdTomato driven by a cardiac enhancer (Klassen et al., 2017) to monitor heart wall movements. To validate roles for the predicted TFBSs, we also performed heart-specific KD of mirr, mad, klu, and HLH3B using this same in vivo approach. As illustrated in Figure 3C, all tested TFs are associated with changes in multiple cardiac traits. mad KD affected FS but also impacted several traits related to rhythmicity. Consistent with the predictions of its binding sites, klu KD altered HP and SI but also EDD and ESD. A PSSM for one isoform of lola (Lola-RC) was predicted from variants associated with rhythmicity, and another one (for Lola-RL) from variants related to contractility. Of note, lola KD strongly impacted all analyzed cardiac phenotypes. luna cardiac-specific RNAi KD had an impact on EDD, FS, HP, and SI, in line with its binding sites predictions. Similarly, srp KD impacted HP, DI, and SI and also modified cardiac diameters. Lastly, mirr KD affected both SI and ESD and HLH3B impacted EDD, HP, and SI (Figure 3C and Supplementary file 2e).

Taken together, the high rate of non-coding SNPs among variants identified in GWAS, the TFBS predictions with in vivo validation of potential cognate TFs, and enrichment for TFs among genes associated with natural variations of heart performance support a central role of GRN deviations in the natural variation of cardiac performance.

Identifying natural variations of phenotypic variability of cardiac performance

The use of inbred lines offers the possibility of studying the heritability of the within-line variance of the phenotypes. Several studies have indeed shown that the within-line environmental variance of inbred lines is heritable (Morgante et al., 2015; Harbison et al., 2013; Ayroles et al., 2015). These genetic variations in phenotypic variability suggest that different genotypes respond differently in response to micro-environmental changes (Geiler-Samerotte et al., 2013; Anholt and Mackay, 2018). We used the Levene test Levene, 1960 to examine heterogeneity of the ‘micro-environmental’ variance among lines for the seven analyzed cardiac traits. We found evidence of significant heterogeneity among the lines in the variance for all analyzed cardiac traits (Table 1). This suggests that the lines differ in their ability to adapt to micro-environmental variations, and supports the heritability of this differing adaptability. The coefficient of within-line variation (CVe) was used as the metric for within-line environmental variance. We found significant genetic variation for all CVe traits (Figure 4A and B and Figure 4—figure supplement 1A) indicating that within some lines, individuals display relatively constant cardiac performance traits, whereas within others, individuals display widely varying cardiac traits. Similar to trait means, CVe traits were poorly correlated with each other (Figure 4—figure supplement 1B), except for HP/DI and EDD/ESD. Correlations of CVe traits with their respective trait means ranged from moderate (FS: –0.53) to none (DI: 0.064, SI: 0.09, ESD: 0.08) (Figure 4—figure supplement 1C). This suggested that distinct loci may affect natural variations of CVe and means of cardiac performance parameters.

Figure 4 with 2 supplements see all
Quantitative genetics and genome-wide associations studies (GWAS) for cardiac traits CVe in the Drosophila Genetic Reference Panel (DGRP).

(A) Distribution of end systolic diameter (ESD) CVe values among DGRP lines. Lines are arranged by increasing CVe values (other phenotypes are displayed in Figure 4—figure supplement 1). (B) Distribution of ESD values across individuals for two representative DGRP lines with low and high intra-genotypic variability. Each dot represents ESD value of a single fly. (Individual values scaled by line mean). (C) Venn diagram illustrating the overlap between single nucleotide polymorphisms (SNPs) (left) and genes (right) associated with mean (yellow) and CVe traits (pink). While only three SNPs were retrieved from GWAS of both mean and CVe traits, genes associated with CVe traits are largely overlapping with those associated with mean traits. (D) Miami plots showing the results of GWAS performed on mean of ESD (top, red) and ESD CVe (bottom, blue). No overlap between sets of variants is observed. Miami plots for other traits are displayed in Figure 4—figure supplement 2. (E) Pearson residuals of chi-square test from the comparison of indicated SNPs categories in the DGRP and among variants associated with cardiac traits CVe (see legend to Figure 1C for detailed description of SNPs categories). (F) Gene Ontology (GO) enrichment analyses of genes associated with CVe traits. Selected molecular functions and cellular components associated with variance of cardiac performances at FDR < 0.05 are shown. (G) Logo representation of position-specific scoring matrices (PSSMs) identified in corresponding sets of non-coding variants associated with cardiac traits CVe. Potential cognate transcription factors (TFs) indicated below were identified by comparing PSSMs to TF binding sites databases.

Figure 4—source data 1

Variants identified associated to phenotypic variance (Cve).

Sheets ‘phenotypes_CVe’: Variants among the 100 best ranked associations for each cardiac trait, with information about the associated genes and the number of Drosophila Genetic Reference Panel (DGRP) lines in which they are present. Variants among the 100 best ranked associations for each cardiac trait that are not within gene mapping criteria (>1 kb from transcription start site [TSS] and transcription end sites [TES]) are shown below. sheets ‘FastEpistasis_CVe_...’: All FastEpistasis data (on traits CVe), per phenotypes. SNPs ID, position, associated genes, and statistics are indicated for both focal SNPs (left) and their interacting SNPs (right). Each sheet displays the results for indicated quantitative traits, except for the first one which is a merge of all quantitative traits association analyses.

https://cdn.elifesciences.org/articles/82459/elife-82459-fig4-data1-v1.xlsx

We performed GWAS analysis to identify candidate variants and genes associated with variation of CVe of cardiac performance using both single marker and epistasis GWAS analyses. Similar to trait means we retained the 100 top ranked variants for reporting single marker associations. These variants were further used as focal variants for epistasis detection and the best ranked interacting SNP for each starting focal SNP was identified. Overall, this led to the identification of 887 variants and 566 genes associated with natural variations of CVe of cardiac performance (Figure 4—source data 1, Supplementary file 3a), which were used for hypothesis generation and functional annotations. Although the number of individuals we were able to analyze – due to the experimental burden of analyzing individual cardiac phenotypes – reduced the statistical power of the analysis of micro-environmental variance, the functional enrichment analyses of variants and genes described below provide support for our association results and highlight important features of the genetic architecture of cardiac traits.

Consistent with the weak correlation between cardiac trait CVe and trait means, we found only three variants associated with both (Figure 4C and D and Figure 4—figure supplement 2; Supplementary file 3c). Two variants in total linkage disequilibrium (LD) (5 bp apart) are localized in an intron of Ca-β, encoding a voltage gated calcium channel subunit, and are associated for both mean and CVe of FS (Figure 4—figure supplement 2, Supplementary file 3c). The other variant, highlighted by both GWAS, was situated at more than 25 kb away from any gene. Despite unraveling largely distinct set of variants, GWAS analyses on CVe and trait means are associated with two sets of genes that overlap significantly (Figure 4C, Supplementary file 3b). Among 566 detected respectively from trait CVe and 562 genes detected by trait mean GWAS, 92 genes are found associated with both (FC=5, hypergeometric pval=2.2 × 10–39). Therefore, a significant number of genes contribute to natural variation of both mean and CVe traits, although through distinct variants.

Within the CVe-associated genes we observed an over-representation of cell adhesion processes, transmembrane signaling, and ion channel activity (pval=1.6 × 10–8 (FC=7), 1.1×10–3 (FC=2,2), and 4.5×10–2 (FC=2,5) respectively). Enriched cellular component annotations included cell-cell junction (pval=5.7 × 10–3 (FC=3.2), receptor complexes (pval=3.3 × 10–4 (FC=4.6), sarcolemma (pval=1.3 × 10–2 (FC=25), as well as several annotations related to neuronal projections suggesting both cell autonomous and non-autonomous activity of related genes (Figure 4F, Supplementary file 3d). In addition, genes leading to stress-induced lethality upon heart-specific KD (Neely et al., 2010) were significantly under-represented in this gene set (FC=0,6; p=0,05; Supplementary file 3e). This may suggest that a significant proportion of CVe GWAS-associated genes affect cardiac trait variance non-cell autonomously. As for trait means, variants localized in the proximal region of the genes (0–1 kb upstream of TSS) are over-represented among SNPs associated with traits CVe (Figure 4E). This suggested that natural variations of cardiac performance variance are also mainly affecting cis-regulatory regions. We also identified over-represented PSSMs in sequences within 150 pb of non-coding CVe variants (Figure 4G). Some of these identified PSSMs were specifically enriched in the sequences around CVe variants. This is the case for instance for Ladybird-early, Lbe, the fly orthologue of the mammalian LBX TF. However, several identified CVe PSSMs are similar to those identified in the vicinity of variants associated with trait means. This is the case for motifs predicted to be bound by Luna and Lola – variants of which are associated with CVe traits (Supplementary file 3a) – and it is also the case for Klu and Mad. Therefore, similar to the situation described for trait means, natural variations of phenotypic plasticity of cardiac performance appears, for a large part, to be driven by modifications of gene regulatory networks.

Conserved features of natural variations of heart performance between flies and humans

Next, we compared our GWAS results for cardiac performance in flies with similar data in humans. A literature survey identified a comprehensive set of human genes associated with cardiac disorders and coronary artery diseases (CAD) in human populations (Supplementary file 3f and h). Using DIOPT (Hu et al., 2011), we identified human orthologues of the genes associated with cardiac trait means and CVe in flies (Supplementary file 4a and d). Compared to the whole fly genome gene set (17,500 genes), both sets of GWAS-associated genes were significantly enriched for genes that have a human orthologue (FC = 1.4, pval = 4.3 · 10–21 for genes associated with trait means and FC = 1.25; pval = 7 · 10–11 for genes associated with traits CVe, Table 2), indicating that, overall, the loci associated with natural variations in cardiac traits have evolutionary conserved functions. More importantly, GWAS loci were enriched for genes whose human orthologues were associated with cardiac disorders (FC = 1.66 [trait means] and 1.7 [traits CVe], pval = 7 · 10–6, Table 2 and Supplementary file 4b and e). Lower enrichment was found for orthologues of human genes associated with CAD (FC = 1.35, pval = 0.04 [trait means] and 1.25, pval = 0.07 [traits CVe]). This analysis strongly supports the view that the cellular and molecular processes of heart development and physiology – which is already well characterized between fly and humans (Ocorr et al., 2007b; Diop and Bodmer, 2015; Rosenthal et al., 2010; Bier and Bodmer, 2004) – are not only conserved but also have a common genetic basis underlying their natural variation.

Table 2
Conserved genes associated with natural variations of cardiac traits from flies to humans.

Enrichment analyses for genes conserved in human and for genes whose human orthologue is associated with either coronary artery diseases (CAD) or cardiac disorders. Numbers in parenthesis indicate the total number of genes in each analyzed gene set (whole fly genome/mean and CVe GWAS-associated genes). First row (human orthologue): Number of fly genes that display a human orthologue according to DIOPT (high and moderate rank). Second and third rows: Number of genes whose human orthologue (high and moderate rank) has been associated with CAD or cardiac disorders by genome-wide associations studies (GWAS) in human populations (Supplementary file 3f-i). Fold change (FC): Ratio between expected (based on successes observed on fly genome) and observed number of successes in respective GWAS gene sets. pval: hypergeometric p-value.

Fly genome (17500)GWAS mean (562)FCpvalGWAS Cve (566)FCpval
Human orthologue94634101.44.3 · 10–213791.257 ·10–11
CAD321191.350.04161.250.07
Cardiac disorders944681.667 · 10–6641.77 · 10–6

Building on these observations, we next wanted to determine whether our data in flies could guide the identification of novel players in the human cardiac system. We focused on the paired box TF gene pox-meso (poxm/PAX9) and on the zinc finger containing TF stripe (sr/EGR2) – both were associated with natural variations of cardiac trait means (Supplementary file 1f). Notably, neither the fly genes nor their human orthologues have a known function in the cardiac system. In flies, poxm contributes to the development and specification of somatic muscles (Duan et al., 2007), while sr is required for induction of tendon cell fate (Gunage et al., 2017). We first tested the effects of heart-specific KD on cardiac performance in flies using the in vivo assay. Both induced an increased heart rate (reduced HP and SI; Figure 5A) but did not affect other cardiac traits analyzed (not shown). The function of their human orthologues PAX9 and EGR2 was then tested in hiPSC-CMs. We asked whether siRNA-mediated KD affected action potential kinetics using a recently developed in vitro assay (Figure 5B; McKeithan et al., 2017). As shown in Figure 5C–E, both gene KDs induced shortening of the AP duration, thus suggesting conserved function in setting cardiac rhythm for both TFs. The effect was stronger for PAX9 KD and seems specific to the repolarization phase, as suggested by APD50 and 75 shortening (Figure 5F–I). APD shortening for PAX9 KD was coincident with increased expression of Na+ and K+ ion channels (SCN5A, KCNH2, and KNCQ1) (Figure 5J), supporting the APD shortening phenotype. In this context, the AP kinetics also correlated with shorter calcium transient duration (Figure 5—figure supplement 1A-D and H-K), including faster upstroke and downstroke calcium kinetics and increased beat rate (peak frequency) (Figure 5—figure supplement 1E-G and L, M), consistent with decreased expression of Calsequestrin 2 isoform (CASQ2) associated with PAX9 KD (Figure 5J). Finally, assessment of the PAX9 KD effect on sarcomeric content revealed an increase in sarcomeric gene expression (Figure 5K), and an upregulation of genes associated with an hypertrophic response (NPPA, NPPB, and NPR1) (Battistoni et al., 2012), which was coincident with increased CM size as quantified by the area of TNNT2 staining per cardiac nuclei (Figure 5L and M).

Figure 5 with 1 supplement see all
In vivo and in vitro assays for poxm/Pax9 and sr/Egr2 in flies and human induced pluripotent stem cells derived cardiomyocyte (hiPSC-CM).

(A) Heart-specific knockdown (KD) of poxm/Pax9 and sr/Egr2 lead to increased heart rate in flies. Cardiac-specific RNAi KD phenotypes of tested transcription factors (TFs). Hand>Gal4; UAS RNAi 1-week adult females (n>30 per cross) were analyzed for cardiac performance using the in vivo assay on intact flies. Scatter plots of individual data, mean, and SD are represented. Wilcoxon test results are indicated (***: pval <0.005). Data were normalized to genotype matched controls. Right panel: representative M-modes (5 s) of RNAi KD and their respective control illustrating the increased heart rate observed upon poxm and sr inactivation. Number of heartbeats counted from M-modes are indicated. (B–I) Effects of Pax9 and Egr2 siRNA on hiPSC-CM action potentials. (B) Schematic overview of single cell and high-throughput voltage assay. (C) Two-dimensional graph for APD75 and Kolgomorov-Smirnov distance (KS-D) representing screen results for PAX9 and EGR2 KD. (D–E) Population distribution of APD75 measurements for siPAX9 and siEGR2 vs. siControl-transfected hiPSC-CMs, respectively. (F) Median action potential traces for siPAX9, siEGR2, and siControl-transfected hiPSC-CMs. (G–I) Histograms showing median APD50, APD75, and APD90 for siPAX9, siEGR2, and siControl-transfected hiPSC-CMs. n > 4 in all experiments. SD are represented. (J,K) Heatmaps of differentially expressed genes in siControl- and siPAX9-transfected hiPSC-CMs. (L) Histogram showing the effect of PAX9 KD on the average TNNT2 staining area per CM as compared to siControl. n > 4 in all experiments. SD are represented. (M) Representative images for siPAX9 and siControl. TNNT2 is shown in red and nuclei are stained with DAPI in blue. Scale bar = 25 μm. t-Test was used to calculate p-values. *p<0.05, ***p<0.001, ****p<0.0001.

Collectively, these data illustrate conserved functions for poxm/PAX9 and sr/EGR2 in setting the cardiac rhythm and identify PAX9 as a novel and key regulator of cardiac performance at the cellular level, via the integrated regulation of expression of genes controlling electrophysiology, calcium handling, and sarcomeric functions in hiPSC-CMs.

Discussion

Disentangling the relative contribution of genetics and environmental factors to cardiac pathologies and the influence of comorbidities is of fundamental importance from a medical perspective. From a fundamental perspective, natural selection is thought to operate on complex traits rather than on single genes, further illustrating the need to characterizing and understanding the genetic architecture of natural variations in cardiac performance. Although mammalian studies provide important information about the molecular and cellular basis of heart development and physiology, these systems are less suited to dissecting the relationship between naturally occurring genetic variations and variations in individual cardiac performance.

Here, we report on the quantitative genetics of cardiac performance in Drosophila. We performed a standard GWAS of variants associated with cardiac trait looking at both within-line and intra-line variance for seven weakly correlated traits that also included an epistasis analysis. The cornerstone of our approach is a series of downstream analyses, which strongly implied that top hits are indeed genetic polymorphisms influencing quantitative variation in cardiac performance. These findings include: (i) the enrichment of pathways known to be related to cardiac development and physiology; (ii) the generation of an extensive network of protein-protein/genetic interactions that plausibly describes the architecture of cardiac performance; (iii) identification of significant overlap between gene occurrences for trait mean and variance; (iv) enrichment of likely causal variants near binding motifs for a set of cardiac-related TFs; (v) observation that several of these TFs are themselves GWAS candidate genes; and (vi) enrichment of fly genes for overlap with known human cardiac-related genes. In addition, we performed validation experiments that showed enrichment of cardiac effects for GWAS identified genes compared to randomly selected control genes. Importantly, we confirmed synergistic genetic interactions between the BMP and activin pathways; and defined two novel regulators of cardiac function that may well be relevant to human disease.

Our results therefore validate Drosophila as a unique model that permits the analysis of the etiology of the genetic foundations of cardiac function. Importantly, we identified many genes not previously associated with cardiac traits at the genetic, cellular, or molecular levels. This suggested that the space for natural variation associated with cardiac performance extends far beyond the genes and gene networks already identified as participating in the establishment and maintenance of cardiac function and illustrates the value of complementing classical genetic studies with an in-depth analysis of the genetic architecture of corresponding complex traits. A large number of candidate genes were validated in vivo using mutations or RNAi KD generated in unrelated genetic backgrounds, suggesting that the underlying biological processes are common to most genotypes.

Most of the variants identified affected non-coding regions of the genome. Recent studies in humans have highlighted that non-coding SNPs may affect the function of enhancers or promoters, thereby influencing the expression levels of surrounding or distant genes and which would explain how some non-coding SNPs may contribute to pathogenesis (Tak and Farnham, 2015; Kikuchi et al., 2019). Although we failed to find enriched potential TFBS at SNPs positions, extending the search space to sequences of ±75 bp around the variants revealed enriched DNA binding motifs, and pointed to potential roles for several TFBS. Our results indicate that variants likely affect sequences outside the core TF binding motif (e.g., flanking sequences) that in turn affect TF binding (Grossman et al., 2017). Supporting our findings, different reports on human GWAS analyses have shown that many risk-associated SNPs are not precisely located in conserved TF binding motifs but are in nearby regions (Heinz et al., 2013; Farh et al., 2015). Several observations indicate that our predictions are accurate: (i) a large number of similar motifs were independently predicted from distinct sets of variants, (ii) some predicted motifs correspond to TFs whose role during cardiogenesis or for the maintenance of cardiac function have already been characterized, (iii) several discovered motifs correspond to TFs for which SNPs were themselves associated with the variation of some of the cardiac traits studied, and (iv) a high percentage of predicted motifs were validated by testing their cognate TFs in vivo for cardiac performance. Of note, while the cardiac GRN has been extensively studied in flies (Seyres et al., 2016; Seyres et al., 2012), a large proportion of the TFBS and cognate factors identified in this study have not been previously characterized as cardiac TFs, again illustrating the complementary contribution represented by the analysis of the genetic architecture of these phenotypes compared to more classical genetic, molecular, and cellular studies.

Although rarely studied, the genetic control of phenotypic variability is of primary importance. Indeed, variability offers adaptive solutions to environmental changes and may impact developmental robustness to micro-environmental perturbations (Geiler-Samerotte et al., 2013). From a medical genetics point of view, highly variable genotypes may produce individuals beyond a phenotypic threshold, leading to disease state (Morgante et al., 2015; Gibson, 2009). Having assessed the cardiac traits for multiple individuals from inbred lines, we were able to examine phenotypic variability among individuals with the same genotypes and therefore to analyze how natural variations impinge on such complex traits. Our results indeed highlight the important contribution that genetic control of variability can play in our understanding of the cause of phenotypic variation. An important feature concerns the nature of the variants and associated genes compared to trait means. In particular, genes identified by a global screen for genes cell autonomously affecting cardiac function (Neely et al., 2010) are under-represented among the genes associated with CVe traits. This suggests that this subset of genes influences the variance of cardiac traits in a cell non-autonomous manner. The separation between variants associated with trait means and those associated with the variance of the same traits was another remarkable feature of the genetic architecture of cardiac phenotypic plasticity. Although based on non-overlapping sets of non-coding regions, predictions of TFBS in the vicinity of non-coding variants identified overlapping sets of motifs. This suggests that natural variation of trait means and trait variance affect partially overlapping GRNs, albeit by distinct polymorphisms. Our analysis allowed us to identify common variants associated with phenotypic variability, and to identify essential characteristics of the genetic architecture of these traits. However, we cannot exclude that part of the observed variability is due also to rare alleles in partial LD with the lead associations (Ek et al., 2018) or to within-line accumulation of mutations. Nevertheless, the former is likely to be marginal, given the rapid decay of LD with local physical distance in Drosophila (Mackay and Huang, 2018).

Several previous studies have established that the principles identified on the DGRP population are universal and can be applied across phyla (Anholt and Mackay, 2018). In particular, the identification of human genes orthologous to those identified by GWAS in the DGRP may allow identification of a translational blueprint and prioritize candidate genes for focused studies in human populations (Zhou et al., 2016; Carbone et al., 2016). Our results provide additional support that the genetic architecture of quantitative traits is conserved. Indeed, a remarkable characteristic of the gene sets associated either with trait means or with trait variance is that they are both enriched for orthologues of human genes that are also associated with disorders affecting the cardiovascular system. This is a strong argument for considering that the genetic architecture of cardiac performance traits is conserved, at least in part, from flies to humans. Hence, the features learned from our study in flies should help reveal the characteristics of the natural variations of cardiac performance in human populations. One intriguing observation is that genes associated with intra-genotypic variance are highly significantly enriched for orthologues of human genes associated with cardiac disorders. Yet, GWAS on human populations are not designed to test the variants involved in intra-genotypic variations. One possible explanation may be that GWAS in humans can pick up variants that do not necessarily affect the mean of cardiac traits, but that cause extreme phenotypes in some individuals by affecting genes that function as buffers.

The high proportion of orthologous genes identified in GWAS in both flies and humans – for cardiac traits in Drosophila and for cardiac disorders in humans – suggested that our results in flies may accelerate discovery of causal genes in humans. We suggest that our study may be used as an aid in prioritizing genes at loci identified by GWAS for cardiac pathologies in humans. Indeed, we show that our data allow identification of new actors that have a conserved role in cardiac function from flies to humans. Although not previously revealed by conventional approaches, we show that poxm and sr – two genes identified by GWAS as being associated with natural variations of cardiac traits – have a conserved function in the regulation of cardiac activity.

In summary, the analysis of natural variations of cardiac performance in Drosophila, combined with functional validations in vivo and in hiPSC-CM, represents a major advance in understanding the mechanisms underlying the genetic architecture of these complex traits, and holds promise for the identification of the genes and mechanisms underlying cardiac disorders in humans.

Materials and methods

Flies’ strains and husbandry

Request a detailed protocol

Fly strains were obtained from the Bloomington Drosophila Stock Center or the Vienna Drosophila Resource Center. Flies were reared at a controlled density on corn flour/yeast/sucrose-agar medium at 25°C, 50% relative humidity and 12 hr light/dark. For movies acquisitions, newly emerged adults (0–24 hr) were collected and aged for 7 days. For screening of the DGRP population, 12 flies were analyzed for each DGRP strain. For validations (RNAi-mediated gene KD or mutant analyses), 30–40 flies per genotype were analyzed. All fly strains used in this study are listed in Supplementary file 4h and i.

Fly heart-specific RNAi KD and mutant analyses

Request a detailed protocol

Driver-line (Hand>Gal4 4.2) (Popichenko et al., 2007) virgins were crossed to UAS-RNAi males or corresponding isogenic control males. BMP mutants were back-crossed into the same isogenic control background. Flies were raised on standard fly food and kept at 25°C. Female progeny were collected and aged to 1 week at 25°C, at which point they were imaged and analyzed.

Semi-intact preparations and SOHA analyses

Request a detailed protocol

To measure cardiac function parameters, denervated, semi-intact Drosophila females adult fly hearts were prepared according to previously described protocols (Ocorr et al., 2007c). High-speed ~3000 frames movies were taken at a rate of 140 frames per second using a Hamamatsu ORCAFlash4.0 digital CMOS camera (Hamamatsu Photonics) on a Zeiss Axioplan microscope with a ×10 water immersion lens. The live images of the heart tube within abdominal A3-4 segments were captured using HCI imaging software (Hamamatsu Photonics). M-modes and cardiac parameters were generated using SOHA (http://www.sohasoftare.com/), as described previously (Ocorr et al., 2007c). M-mode records provide a snapshot of the movement of heart wall over time. Cardiac parameters were calculated as below: DI is the period when the heart is completely relaxed (diastole). SI is the period when the heart is actively contracted. HP is the time between the two consecutive diastoles (i.e., DI+SI) and heart rate is calculated from HP (1/HP). AI is the standard deviation of the HP mean for each fly normalized to the median HP. FS was calculated as (diastolic diameter – systolic diameter/diastolic diameter).

In vivo analysis of heart function (tdtk)

Request a detailed protocol

R94C02::tdTomato (attP2) (Klassen et al., 2017) was combined with Hand-Gal4 as stable stock. Flies were then crossed to RNAi lines and reared to 3 weeks of age. To assay fluorescent hearts in vivo flies were mounted as described (Klassen et al., 2017) and were illuminated with green light (3 mW power). Five s movies of the beating heart were acquired at 270 fps and high-speed recordings were then processed using a custom R script (Vogler, 2022a).

Source code for analyses

Request a detailed protocol

We have developed a GitHub project containing the instructions and material to reproduce the analyses (Saha et al., 2021). Source code are available in the GitHub repository. Required data and built singularity images are available on download. Instructions to reproduce the analyses are provided.

Quantitative genetic analyses

Request a detailed protocol

Phenotype data were analyzed separately to perform quality control. In each line, data outside the range [Q1 − 1.5*IQR; Q3+1.5*IQR] were identified (where Q1/3 are first/third quartile and IQR = Q3 −Q1). Most of the films corresponding to these outlier data were re-analyzed. It was identified that these extreme results corresponded to technical problems (poor quality film, bad placement of the measurement points, etc.) and were removed from the analysis to avoid technical bias.

Mean and coefficient of variation (CV) were computed for each line and each phenotype. Only lines with at least seven remaining observations after QC were kept for the computation.

Correlation between phenotypes was computed using Spearman correlation between mean/CV of lines. Genetic indicators, like broad sense heritability, were computed from the QC data.

Genotype-phenotype associations: GWAS and epistasis study

Request a detailed protocol

We used mean and CV summary statistics from our GWAS analysis to study the leverage effect of extreme values. We identified strong leverage effects in some phenotypes and treated them for outliers. We decided to remove lines with large leverage effect from GWAS analysis of the phenotype, since they imply a strong bias in the linear model. Hence, mean/CV values outside the range of [Q1 − 1.5*IQR; Q3+1.5*IQR] were discarded from GWAS analysis of the corresponding phenotype, thus removing all leverage effects.

GWAS analysis was performed using FastLMM (Lippert et al., 2011). Linear mixed model was built using Wolbachia infection and common polymorphic inversions (data provided by the DGRP Consortium) as covariates. A GitHub project containing the instructions and material to reproduce the analyses reported is available at: https://doi.org/10.5281/zenodo.5582846.

A QQ plot was generated to display the quantile distribution of observed marker-phenotype association p-values vs. the distribution of expected p-values by chance. The QQ plot was plotted to check if some polymorphisms are more associated with the phenotype under a null hypothesis, which would be represented by a uniform distribution of p-values. For each phenotype, we selected the 100 variants with the lowest GWAS p-value. These variants were mapped to genes using the following rules: (i) if the variant was located in a gene (INTRON, EXON, 3’or 5’ UTR), the SNP was then associated with this gene, (ii) if not, the variant was associated with a gene if its distance from the boundary of the gene (TSS and TES was less than 1 kb); (iii) otherwise the variant was not mapped to a gene. Of the mapped variants we then examined variants with an MAF greater than 4% and discarded SNPs with synonymous effects. Epistasis analysis was performed using FastEpistasis (Schüpbach et al., 2010). For each phenotype, we used the SNP selected from the GWAS analysis as focal SNP and the interaction was tested with all the SNPs with MAF greater than 5% in the DGRP population. The best matching test SNP was selected for each focal SNP. Each resulting SNP was mapped to genes using the same rules as for GWAS variants.

Miami plots were designed according to Vogler, 2022a.

Interaction network analysis

Request a detailed protocol

Drosophila interactome: A network generated from genetic and physical interaction databases has been built. Briefly, known protein-protein interactions, identified using Y2H and AP-MS techniques, and genetic interactions, are curated in the DroID database (http://www.droidb.org/) (Murali et al., 2011), to which the Y2H FlyBi dataset (https://flybi.hms.harvard.edu/results.php) has been added. Overall, this Drosophila interaction network contains 79,698 interactions between 11,022 genes/proteins. Network representations were made using Cytoscape (Shannon et al., 2003).

Statistics

Enrichment p-values were based on a test following the hypergeometric distribution. Enrichment (FC) were calculated as the ratio between observed and expected successes in samples (based on the observations in the population). GO enrichment analyses were conducted using the online version of G:profiler tool (Raudvere et al., 2019) with default parameters, including g:SCS method for computing multiple testing corrections for p-values gained from GO and pathway enrichment analysis.

(https://biit.cs.ut.ee/gprofiler/page/docs#significance_threhshold).

Statistical analysis of cardiac parameters in RNAi and mutant experiments were performed using R. The non-parametric Wilcoxon-Mann-Whitney test was used to compare the median of the difference between two groups (tested genes/respective control).

Comparing Z-scores (Figure 2—figure supplement 1B): For each phenotype and each gene, Z-scores summarizes the difference in the phenotype distribution from tested genes to their respective control. Z-scores were calculated for all phenotypes (7 phenotypes) and all tested genes (GWAS-associated genes: 42, randomly selected genes: 18). The Z-scores for all the phenotypes of the tested genes were pooled together (GWAS-associated genes: 294 Z-scores; randomly selected genes: 126 Z-scores) and the distribution of absolute Z-scores between the two groups was analyzed using Mann-Whitney-Wilcoxon test.

Tests for genetic interactions (Figure 2E and Figure 2—figure supplement 2): Two-way ANOVA with interactions was conducted to test if the interaction of two genes is significantly different from the two-way ANOVA model with the additive effects of the two genes (Kim, 2014).

Phenotypic correlation between each trait pair was computed using Spearman correlation. Pearson chi-squared test was applied to test if the genomic location of variants associated with cardiac traits is biased toward any particular genomic region when compared against the whole set of variants in the entire genome of the DGRP population. The overlap coefficient/overlap index is a similarity measure that measures the overlap between two finite sets. This is used to quantify the overlap between the gene sets associated with the different cardiac traits.

Motif’s analysis

Request a detailed protocol

For every set of SNPs (rhythmicity CVe and mean, contractility CVe, and mean), the coordinates (dm3 genome assembly) were extended 75 nt to both sides and their corresponding fasta sequences were retrieved. TF binding motifs were identified using RSAT peak motifs (Thomas-Chollier et al., 2012) (default parameter) using two distinct background models: intrinsic (frequencies of short k-mers in the input sequences to calculate the enrichment of longer k-mers) and random (set of random genomic regions matching length and nucleotide composition using biasaway; Worsley Hunt et al., 2014). To reduce the redundancy of the found motifs (PSSMs), we applied the RSAT matrix-clustering (Castro-Mondragon et al., 2017) algorithm (default parameters) to find clusters of motifs. The clustered motifs were manually annotated by comparison with PSSM databases (Cis-BP Weirauch et al., 2014, JASPAR Fornes et al., 2020, and FlyFactorSurvey Zhu et al., 2011) using the tool SAT compare matrices (Nguyen et al., 2018) (default parameters).

Generation of hiPSC-CMs

Request a detailed protocol

Id1 overexpressing hiPSCs (Cunningham et al., 2017) were dissociated with 0.5 mM EDTA (Thermo Fisher Scientific) in PBS without CaCl2 and MgCl2 (Corning) for 7 min at room temperature. hiPSC were resuspended in mTeSR-1 media (STEMCELL Technologies) supplemented with 2 µM Thiazovivin (STEMCELL Technologies) and plated in a Matrigel-coated 12-well plate at a density of 3×105 cells per well. After 24 hr after passage, cells were fed daily with mTeSR-1 media (without Thiazovivin) for 3–5 days until they reached ≥90% confluence to begin differentiation. hiPSC-ACMs were differentiated as previously described (Burridge et al., 2015). At day 0, cells were treated with 6 µM CHIR99021 (Selleck Chemicals) in S12 media (Pei et al., 2017) for 48 hr. At day 2, cells were treated with 2 µM Wnt-C59 (Selleck Chemicals), an inhibitor of WNT pathway, in S12 media. Forty-eight hr later (at day 4), S12 media is fully changed. At day 5, cells were dissociated with TrypLE Express (Gibco) for 2 min and blocked with RPMI (Gibco) supplemented with 10% FBS (Omega Scientific). Cells were resuspended in S12 media supplemented with 4 mg/L Recombinant Human Insulin (Gibco) (S12+media), 2 µM Thiazovivin and plated onto a Matrigel-coated 12-well plate at a density of 9×105 cells per well. S12+media was changed at day 8 and replaced at day 10 with RPMI (Gibco) media supplemented with 213 µg/µL L-ascorbic acid (Sigma), 500 mg/L BSA-FV (Gibco), 0.5 mM L-carnitine (Sigma), and 8 g/L AlbuMAX Lipid-Rich BSA (Gibco) (CM media). Typically, hiPSC-ACMs start to beat around day 10. At day 15, cells were purified with lactate media (RPMI without glucose, 213 µg/µL L-ascorbic acid, 500 mg/L BSA-FV, and 8 mM Sodium-DL-Lactate, Sigma), for 4 days. At day 19, media was replaced with CM media.

Voltage assay in hiPSC-CMs

Request a detailed protocol

Voltage assay is performed using labeling protocol described in McKeithan et al., 2017. Briefly, hiPSC-CMs at day 25 of differentiation were dissociated with TrypLE Select ×10 for up to 10 min and action of TrypLE was neutralized with RPMI supplemented with 10% FBS. Cells were resuspended in RPMI with 2% KOSR (Gibco) and 2% B27 ×50 with vitamin A (Life Technologies) supplemented with 2 µM Thiazovivin and plated at a density of 6000 cells per well in a Matrigel-coated 384-well plate. hiPSC-CMs were then transfected with siRNAs directed against PAX9 and EGR2 (ON-TARGETplus Human) using Lipofectamine RNAi Max (Thermo Fisher Scientific). Each siRNA was tested individually in 8-plicates. Three days post-transfection, cells were first washed with pre-warmed Tyrode’s solution (Sigma) by removing 50 µL of media and adding 50 µL of Tyrode’s solution five times using a 16-channel pipette. After the fifth wash, 50 µL of ×2 dye solution consisting in voltage-sensitive dye Vf2.1 Cl (Fluovolt, 1:2000, Thermo Fisher Scientific) diluted in Tyrode’s solution supplemented with 1 µL of 10% Pluronic F127 (diluted in water, Thermo Fisher Scientific) and 20 µg/mL Hoescht 33258 (diluted in water, Thermo Fisher Scientific) was added to each well. The plate was placed back in the 37°C 5% CO2 incubator for 45 min. After incubation time, cells were washed four times with fresh pre-warmed Tyrode’s solution using the same method described above. hiPSC-CMs were then automatically imaged with ImageXpress Micro XLS microscope at an acquisition frequency of 100 Hz for a duration of 5 s with excitation wavelength of 485/20 nm and emission filter 525/30 nm. A single image of Hoescht was acquired before the time series. Fluorescence over time quantification and trace analysis were automatically quantified using custom software packages developed by Molecular Devices and Colas lab. Two independent experiments were performed.

Calcium assay in hiPSC-CMs

Request a detailed protocol

Calcium assay is performed using labeling protocol as described (Cerignoli et al., 2012). Cells were prepared and transfected with siRNAs against PAX9 and EGR2 as described above in the voltage assay. Three days post-transfection, 50 µL of media was removed and replaced in each well by a ×2 calcium dye solution consisting in Fluo-4 NW dye (Invitrogen), 1.25 mM Probenecid F-127 (Invitrogen), and 100 µg/mL Hoescht 33258 (diluted in water, Thermo Fisher Scientific) diluted in warm Tyrode’s solution (Sigma). The plate was placed back in the 37°C 5% CO2 incubator for 45 min. After incubation time, cells were washed four times with fresh pre-warmed Tyrode’s solution and imaged as described above.

RNA-seq and data analysis

Request a detailed protocol

Day 25 hiPSC-CMs were transfected with 25nM final concentrations of siRNA against PAX9 and with scrambled control siRNAs. After 3days, cells were pelleted and resuspended in 500µL TRIzol reagent followed by total RNA extraction. Library preparation and sequencing of the samples was done at La Jolla Institute of Immunology (La Jolla, CA). FASTQ files were processed using nf-core/rnaseq (Ewels et al., 2021) with Nextflow Version: 21.03.0.edge. Differential gene expression was determined using R/DESeq2 (Love et al., 2014) and GO term enrichment was done using gprofiler2 (Kolberg and Raudvere, 2021). All analysis scripts can be downloaded at https://github.com/gvogler/eLife-2022-Saha-et-al (copy archived at swh:1:rev:b102c48b85976216c37d2e5aed060670535192f6; Vogler, 2022b). All sequencing data were deposited on the Gene Expression Omnibus (GEO) repository (accession number GSE217655).

Immunostaining of hiPSC-CMs

Request a detailed protocol

Three days post-transfection cells were fixed with 4% paraformaldehyde for 30 min and blocked in blocking buffer (10% Horse Serum, 10% Gelatin, and 0.5% Triton X-100) for 20 min. hPSC-CMs were stained with sarcomeric protein TNNT2 (Catalog #: HPA017888 Lot #: A91786, Sigma, dilution 1/200), secondary antibody Goat anti-Rabbit IgG Alexa Fluor 568 (Invitrogen, 1/500) and DAPI (1/1000) in Blocking Buffer. Cells were imaged with ImageXpress Micro XLS microscope (Molecular Devices) and custom algorithms were used to quantify the area of TNNT2 mask per cardiac nuclei. All experiments were performed at least in biological quadruplicates.

Data availability

All data generated or analysed during this study are included in the manuscript and supporting file and available on Zenodo (Saha et al., 2021).

References

  1. Book
    1. Levene H
    (1960)
    Contributions to probability and statistics: essays in honor of harold hotelling
    In: Li R, editors. Stanford Studies in Mathematics and Statistics. Stanford, CA: Stanford Univ Press. pp. 278–292.
    1. Roselli C
    2. Chaffin MD
    3. Weng L-C
    4. Aeschbacher S
    5. Ahlberg G
    6. Albert CM
    7. Almgren P
    8. Alonso A
    9. Anderson CD
    10. Aragam KG
    11. Arking DE
    12. Barnard J
    13. Bartz TM
    14. Benjamin EJ
    15. Bihlmeyer NA
    16. Bis JC
    17. Bloom HL
    18. Boerwinkle E
    19. Bottinger EB
    20. Brody JA
    21. Calkins H
    22. Campbell A
    23. Cappola TP
    24. Carlquist J
    25. Chasman DI
    26. Chen LY
    27. Chen Y-DI
    28. Choi E-K
    29. Choi SH
    30. Christophersen IE
    31. Chung MK
    32. Cole JW
    33. Conen D
    34. Cook J
    35. Crijns HJ
    36. Cutler MJ
    37. Damrauer SM
    38. Daniels BR
    39. Darbar D
    40. Delgado G
    41. Denny JC
    42. Dichgans M
    43. Dörr M
    44. Dudink EA
    45. Dudley SC
    46. Esa N
    47. Esko T
    48. Eskola M
    49. Fatkin D
    50. Felix SB
    51. Ford I
    52. Franco OH
    53. Geelhoed B
    54. Grewal RP
    55. Gudnason V
    56. Guo X
    57. Gupta N
    58. Gustafsson S
    59. Gutmann R
    60. Hamsten A
    61. Harris TB
    62. Hayward C
    63. Heckbert SR
    64. Hernesniemi J
    65. Hocking LJ
    66. Hofman A
    67. Horimoto ARVR
    68. Huang J
    69. Huang PL
    70. Huffman J
    71. Ingelsson E
    72. Ipek EG
    73. Ito K
    74. Jimenez-Conde J
    75. Johnson R
    76. Jukema JW
    77. Kääb S
    78. Kähönen M
    79. Kamatani Y
    80. Kane JP
    81. Kastrati A
    82. Kathiresan S
    83. Katschnig-Winter P
    84. Kavousi M
    85. Kessler T
    86. Kietselaer BL
    87. Kirchhof P
    88. Kleber ME
    89. Knight S
    90. Krieger JE
    91. Kubo M
    92. Launer LJ
    93. Laurikka J
    94. Lehtimäki T
    95. Leineweber K
    96. Lemaitre RN
    97. Li M
    98. Lim HE
    99. Lin HJ
    100. Lin H
    101. Lind L
    102. Lindgren CM
    103. Lokki M-L
    104. London B
    105. Loos RJF
    106. Low S-K
    107. Lu Y
    108. Lyytikäinen L-P
    109. Macfarlane PW
    110. Magnusson PK
    111. Mahajan A
    112. Malik R
    113. Mansur AJ
    114. Marcus GM
    115. Margolin L
    116. Margulies KB
    117. März W
    118. McManus DD
    119. Melander O
    120. Mohanty S
    121. Montgomery JA
    122. Morley MP
    123. Morris AP
    124. Müller-Nurasyid M
    125. Natale A
    126. Nazarian S
    127. Neumann B
    128. Newton-Cheh C
    129. Niemeijer MN
    130. Nikus K
    131. Nilsson P
    132. Noordam R
    133. Oellers H
    134. Olesen MS
    135. Orho-Melander M
    136. Padmanabhan S
    137. Pak H-N
    138. Paré G
    139. Pedersen NL
    140. Pera J
    141. Pereira A
    142. Porteous D
    143. Psaty BM
    144. Pulit SL
    145. Pullinger CR
    146. Rader DJ
    147. Refsgaard L
    148. Ribasés M
    149. Ridker PM
    150. Rienstra M
    151. Risch L
    152. Roden DM
    153. Rosand J
    154. Rosenberg MA
    155. Rost N
    156. Rotter JI
    157. Saba S
    158. Sandhu RK
    159. Schnabel RB
    160. Schramm K
    161. Schunkert H
    162. Schurman C
    163. Scott SA
    164. Seppälä I
    165. Shaffer C
    166. Shah S
    167. Shalaby AA
    168. Shim J
    169. Shoemaker MB
    170. Siland JE
    171. Sinisalo J
    172. Sinner MF
    173. Slowik A
    174. Smith AV
    175. Smith BH
    176. Smith JG
    177. Smith JD
    178. Smith NL
    179. Soliman EZ
    180. Sotoodehnia N
    181. Stricker BH
    182. Sun A
    183. Sun H
    184. Svendsen JH
    185. Tanaka T
    186. Tanriverdi K
    187. Taylor KD
    188. Teder-Laving M
    189. Teumer A
    190. Thériault S
    191. Trompet S
    192. Tucker NR
    193. Tveit A
    194. Uitterlinden AG
    195. Van Der Harst P
    196. Van Gelder IC
    197. Van Wagoner DR
    198. Verweij N
    199. Vlachopoulou E
    200. Völker U
    201. Wang B
    202. Weeke PE
    203. Weijs B
    204. Weiss R
    205. Weiss S
    206. Wells QS
    207. Wiggins KL
    208. Wong JA
    209. Woo D
    210. Worrall BB
    211. Yang P-S
    212. Yao J
    213. Yoneda ZT
    214. Zeller T
    215. Zeng L
    216. Lubitz SA
    217. Lunetta KL
    218. Ellinor PT
    (2018) Multi-Ethnic genome-wide association study for atrial fibrillation
    Nature Genetics 50:1225–1233.
    https://doi.org/10.1038/s41588-018-0133-9
  2. Book
    1. Rosenthal N
    2. Harvey RP
    3. Bodmer R
    4. Frasch M
    (2010)
    Development and aging of the Drosophila heart
    In: Rosenthal N, editors. Heart Development and Regeneration. Elsevier. pp. 47–86.
    1. Shah S
    2. Henry A
    3. Roselli C
    4. Lin H
    5. Sveinbjörnsson G
    6. Fatemifar G
    7. Hedman ÅK
    8. Wilk JB
    9. Morley MP
    10. Chaffin MD
    11. Helgadottir A
    12. Verweij N
    13. Dehghan A
    14. Almgren P
    15. Andersson C
    16. Aragam KG
    17. Ärnlöv J
    18. Backman JD
    19. Biggs ML
    20. Bloom HL
    21. Brandimarto J
    22. Brown MR
    23. Buckbinder L
    24. Carey DJ
    25. Chasman DI
    26. Chen X
    27. Chen X
    28. Chung J
    29. Chutkow W
    30. Cook JP
    31. Delgado GE
    32. Denaxas S
    33. Doney AS
    34. Dörr M
    35. Dudley SC
    36. Dunn ME
    37. Engström G
    38. Esko T
    39. Felix SB
    40. Finan C
    41. Ford I
    42. Ghanbari M
    43. Ghasemi S
    44. Giedraitis V
    45. Giulianini F
    46. Gottdiener JS
    47. Gross S
    48. Guðbjartsson DF
    49. Gutmann R
    50. Haggerty CM
    51. van der Harst P
    52. Hyde CL
    53. Ingelsson E
    54. Jukema JW
    55. Kavousi M
    56. Khaw KT
    57. Kleber ME
    58. Køber L
    59. Koekemoer A
    60. Langenberg C
    61. Lind L
    62. Lindgren CM
    63. London B
    64. Lotta LA
    65. Lovering RC
    66. Luan J
    67. Magnusson P
    68. Mahajan A
    69. Margulies KB
    70. März W
    71. Melander O
    72. Mordi IR
    73. Morgan T
    74. Morris AD
    75. Morris AP
    76. Morrison AC
    77. Nagle MW
    78. Nelson CP
    79. Niessner A
    80. Niiranen T
    81. O’Donoghue ML
    82. Owens AT
    83. Palmer CNA
    84. Parry HM
    85. Perola M
    86. Portilla-Fernandez E
    87. Psaty BM
    88. Rice KM
    89. Ridker PM
    90. Romaine SPR
    91. Rotter JI
    92. Salo P
    93. Salomaa V
    94. van Setten J
    95. Shalaby AA
    96. Smelser DT
    97. Smith NL
    98. Stender S
    99. Stott DJ
    100. Svensson P
    101. Tammesoo ML
    102. Taylor KD
    103. Teder-Laving M
    104. Teumer A
    105. Thorgeirsson G
    106. Thorsteinsdottir U
    107. Torp-Pedersen C
    108. Trompet S
    109. Tyl B
    110. Uitterlinden AG
    111. Veluchamy A
    112. Völker U
    113. Voors AA
    114. Wang X
    115. Wareham NJ
    116. Waterworth D
    117. Weeke PE
    118. Weiss R
    119. Wiggins KL
    120. Xing H
    121. Yerges-Armstrong LM
    122. Yu B
    123. Zannad F
    124. Zhao JH
    125. Hemingway H
    126. Samani NJ
    127. McMurray JJV
    128. Yang J
    129. Visscher PM
    130. Newton-Cheh C
    131. Malarstig A
    132. Holm H
    133. Lubitz SA
    134. Sattar N
    135. Holmes MV
    136. Cappola TP
    137. Asselbergs FW
    138. Hingorani AD
    139. Kuchenbaecker K
    140. Ellinor PT
    141. Lang CC
    142. Stefansson K
    143. Smith JG
    144. Vasan RS
    145. Swerdlow DI
    146. Lumbers RT
    147. Regeneron Genetics Center
    (2020) Genome-wide association and mendelian randomisation analysis provide insights into the pathogenesis of heart failure
    Nature Communications 11:163.
    https://doi.org/10.1038/s41467-019-13690-5
    1. van Setten J
    2. Brody JA
    3. Jamshidi Y
    4. Swenson BR
    5. Butler AM
    6. Campbell H
    7. Del Greco FM
    8. Evans DS
    9. Gibson Q
    10. Gudbjartsson DF
    11. Kerr KF
    12. Krijthe BP
    13. Lyytikäinen LP
    14. Müller C
    15. Müller-Nurasyid M
    16. Nolte IM
    17. Padmanabhan S
    18. Ritchie MD
    19. Robino A
    20. Smith AV
    21. Steri M
    22. Tanaka T
    23. Teumer A
    24. Trompet S
    25. Ulivi S
    26. Verweij N
    27. Yin X
    28. Arnar DO
    29. Asselbergs FW
    30. Bader JS
    31. Barnard J
    32. Bis J
    33. Blankenberg S
    34. Boerwinkle E
    35. Bradford Y
    36. Buckley BM
    37. Chung MK
    38. Crawford D
    39. den Hoed M
    40. Denny JC
    41. Dominiczak AF
    42. Ehret GB
    43. Eijgelsheim M
    44. Ellinor PT
    45. Felix SB
    46. Franco OH
    47. Franke L
    48. Harris TB
    49. Holm H
    50. Ilaria G
    51. Iorio A
    52. Kähönen M
    53. Kolcic I
    54. Kors JA
    55. Lakatta EG
    56. Launer LJ
    57. Lin H
    58. Lin HJ
    59. Loos RJF
    60. Lubitz SA
    61. Macfarlane PW
    62. Magnani JW
    63. Leach IM
    64. Meitinger T
    65. Mitchell BD
    66. Munzel T
    67. Papanicolaou GJ
    68. Peters A
    69. Pfeufer A
    70. Pramstaller PP
    71. Raitakari OT
    72. Rotter JI
    73. Rudan I
    74. Samani NJ
    75. Schlessinger D
    76. Silva Aldana CT
    77. Sinner MF
    78. Smith JD
    79. Snieder H
    80. Soliman EZ
    81. Spector TD
    82. Stott DJ
    83. Strauch K
    84. Tarasov KV
    85. Thorsteinsdottir U
    86. Uitterlinden AG
    87. Van Wagoner DR
    88. Völker U
    89. Völzke H
    90. Waldenberger M
    91. Jan Westra H
    92. Wild PS
    93. Zeller T
    94. Alonso A
    95. Avery CL
    96. Bandinelli S
    97. Benjamin EJ
    98. Cucca F
    99. Dörr M
    100. Ferrucci L
    101. Gasparini P
    102. Gudnason V
    103. Hayward C
    104. Heckbert SR
    105. Hicks AA
    106. Jukema JW
    107. Kääb S
    108. Lehtimäki T
    109. Liu Y
    110. Munroe PB
    111. Parsa A
    112. Polasek O
    113. Psaty BM
    114. Roden DM
    115. Schnabel RB
    116. Sinagra G
    117. Stefansson K
    118. Stricker BH
    119. van der Harst P
    120. van Duijn CM
    121. Wilson JF
    122. Gharib SA
    123. de Bakker PIW
    124. Isaacs A
    125. Arking DE
    126. Sotoodehnia N
    (2018) Pr interval genome-wide association meta-analysis identifies 50 loci associated with atrial and atrioventricular electrical activity
    Nature Communications 9:2904.
    https://doi.org/10.1038/s41467-018-04766-9

Decision letter

  1. Detlef Weigel
    Senior and Reviewing Editor; Max Planck Institute for Biology Tübingen, Germany

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

[Editors' note: this paper was reviewed by Review Commons.]

https://doi.org/10.7554/eLife.82459.sa1

Author response

Overall, we acknowledge referee’s careful reading of the paper and comments that we think have helped further improvement of the manuscript.

On the attached pages are our detailed point by point responses to the referees’ comments along with a description of how the manuscript was modified in accordance.

New data included:

In response to the comments and suggestions of both reviewers 1 and 3, we conducted new experiments to test genetic interactions between different actors of the BMP and activin pathways. These new results confirm and complement the analyses described in the original manuscript. Furthermore, as suggested by reviewer 2, we have further studied the phenotypes of hiPSC-CM, by analyzing gene expression profiles and by analyzing the morphological changes induced as a result of PAX9 knockdown.

NB: The title has been slightly modified, to highlight the conserved features of the genetic architecture of cardiac performance revealed in the study

Former title: Genetic architecture of natural variation of cardiac performance in flies.

Novel title: Genetic architecture of natural variation of cardiac performance: From flies to humans.

Reviewer 1

1. The authors utilized the RNAi-mediated knockdown approach in their functional validation studies. It is not clear how each genetic variation (SNP) affects its associated genes. Could some of the SNPs activate the candidate gene expression? For the 4 candidate genes that failed to show cardiac defects, could the overexpression of these 4 genes alter cardiac performance?

Of course, we cannot predict direction of the effect of the variants on the function of the genes. In this context, loss-of-function experiments are subjected to a risk of false negatives. It is indeed possible that in the case of a lack of effect of the loss of function, a gain of function could reveal an effect. But gain-of-function experiments are difficult to control, and often subjected to non-specific effects because it is complicated to control the level of over-expression compared to endogenous expression. This did not seem suitable for an extensive analysis of a large number of genes. We therefore chose to test only for loss of function.

In addition, our approach to testing heart-specific RNAi aims to assess the quality of the association results by comparing RNAi for genes identified by GWAS to randomly selected genes. It is not intended to describe precisely the involvement of each gene individually.

(See also answer to reviewer 2 comment n°2 and the modifications to the manuscript that have been made and which address these criticism).

2. babo is the type I activin receptor, not type 2.

Thank you, we have corrected this error.

3. The authors show BMP and activin pathway genetically interacts to affect cardiac performance. But it is interesting to find that these interactions are in a trait-dependent manner. For example, it seems that babo and dpp epistatically interact to regulate FS, while they additively regulate HP and DI. The authors need to discuss the complex genetic interaction further.

See reply to reviewer 3, comment N°2 below.

4. Both snoo and sog are identified from GWAS. How about babo and dpp? Are there any identified SNPs associated with babo and dpp?

Considering GWAS for mean phenotypes, there is no variant in dpp that are within the 100 best ranked SNPs nor within the variants identified using fast epistasis. But given the size of the DGRP population we are far from being exhaustive, as we do not reach saturation. It is therefore difficult to comment on these ‘negative’ results. However, we do identify one variant in babo using fast epistasis (see figure 2B and Table S3).

5. It is unclear why the mad KD behaves oppositely to dpp mutant, although both proteins are involved in BMP pathway. In Figure S5, the mad KD shows reduced FS and HP, but dpp LOF mutant shows increased FS and HP (Figure S4). Can the authors perform RNAi to knockdown dpp specific in the heart to reexamine the role of dpp in the regulation of cardiac function. The whole body LOF mutant dpp-d14 might not target cardiac tissue directly to control heart performance like mad KD.

(See also answer to reviewer 3 comment n°2) We did perform heart specific dpp RNAi experiments together with other tests for interactions using new allelic combinations of activin and BMP pathways and therefore can compare heart specific knock down to heterozygotes for amorphic mutations for both dpp and mad.

Regarding dpp, congruent effects on HP, DI, SI, ESD and EDD were observed between mutant and RNAi, while RNAi had opposite effects on FS compared to heterozygotes dppd14 mutants (decreased and increased FS compared to control, respectively). In the case of mad, heterozygous mutants had no effect on FS, EDD and ESD, but similarly to dpp mutants it increased SI, DI and HP. mad RNAi uniquely decreased HP, DI and SI and increased AI. However, similarly to dpp RNAi, it induced a decrease of FS.

Thus, systemic versus heart specific knockdown of genes induce specific effects, suggesting cardiac non-autonomous interactions. This complex picture of TGFb involvement is now discussed in the result section (see below, Reviewer 3, major comment 2).

6. The authors selected two novel genes to study the conversed regulation in both flies and human iPSC cells. Besides testing these novel genes, the authors should also verify whether the conserved pathways, like TGF-β, regulate heart performance in human iPSC cells similar to the flies.

We focused on poxm/Pax9 and sr/Egr2 because none of these TFs were known to have cardiac function in fly nor in mammals. Our paralleled analyses in fly and hiPS-CM illustrates how the description of the genetic architecture of cardiac traits in flies can accelerate discovery in mammals.

There is extensive literature describing the involvement of TGF B /BMP and Activin pathways in heart development and diseases in humans, hence the choice not to focus on these pathways in iPS-CM.

Reviewer 2:

1. It will be interesting to compare this fly GWAS to human heart disease GWAS data (for example, cardiomyopathy, arrhythmia, heart failure) from patients. Such cross comparison could make the data set more valuable.

We actually did make this comparison (Table 2, Table S11) and we agree it significantly validates our approach. This identified a set of orthologous genes associated with cardiac traits both in Drosophila and humans, supporting the conservation of the genetic architecture of cardiac performance traits, from arthropods to mammals.

2. RNAi is the only experimental approach in this manuscript to validate the functional significance from data analyses. Authors may consider using genetic mutations such as deficiency lines or P-element lines to offer an alternative approach. This is simply a suggestion to improve the rigor and reproducibility, not absolutely required.

In an attempt to provide a consistent analysis of loss of gene function, our strategy was to concentrate our analysis on the effects of heart specific knock down. This allows us to compare -in a global way- the effects of the knock down of genes identified by GWAS to those of randomly selected genes.

Our objective was to provide a global view of the heart specific effects of the identified genes, and not to characterize precisely the involvement of each of them, using a combination of mutant alleles, RNAi and gain of function. Given the experimental burden of analyzing cardiac function, such a strategy would have indeed required us to concentrate only a very small number of genes.

We however recognize that this strategy has limitations:

– Some variants may lead to gain-of-function effects of genes, and our strategy is not able to test for these effects.

– Some variants may come from non-cell-autonomous effects, which would not be replicated by our targeted RNAi strategy in the heart.

Therefore, the false negative rate of our experiments is difficult to estimate.

We have tried to put this into perspective and to highlight the limitations of our analysis in the Results section describing RNAi validation of GWAS results.

“To assess in an extensive way whether mutations in genes harboring SNPs associated with variation in cardiac traits contributed to these phenotypes […] These results therefore supported our association results. It is important to emphasize that our approach is limited to testing the effect of tissue-specific gene knock down. Since some of the variants may lead to increased gene function and/or expression, this can lead to a false negative rate that is difficult to estimate. In addition, some of the associated variants may influence heart function by non cell-autonomous mechanisms, which would not be replicated by cardiac specific RNAi knock down.”

3. In order to validate the roles of predicted TF binding sites, the best approach would be introducing point mutations using CRISPR/Cas9 within the binding motif then testing out molecular and physiological outcomes. Rather authors chose to test indirectly to knock down those TFs. If so, authors need to at least acknowledge the potential caveats of such approach and the limitation in related data interpretation.

The reviewer is right, the definitive proof of the involvement of a potential TF binding site on the regulation of a gene located in cis requires to mutate the binding site and to analyze the effect on the expression of the corresponding gene. But this may not be sufficient to definitely demonstrate that the potential TF is indeed a regulator of that gene (the binding motif may be target of yet another TF): definitive proof may require motifs/TF DNA binding domain swaps. This would have been out of the scope of the present study. In addition, the effects on heart performance of mutating one TFBS at a time (among several dozens) may be too weak to allow their characterization with available tools and approaches.

We acknowledge however that our approach provides an indirect validation of transcription factors binding sites predictions. This was, in our opinion, the most efficient way to evaluate the potential effect of predicted transcription factors.

We clarify this in the result section:

“We did not test individually the effects on cardiac performance of mutations in predicted TFBSs located near the SNPs because any individual effect would probably be too small to be detectable by the available methods. Rather, we tested the potential involvement of their cognate TFs by cardiac specific RNAi mediated KD”

4. hiPSC-CM data is somewhat limited by only showing the HR and AP duration data. It is recommended to include some immunocytochemistry data to show the morphology, sarcomere structure of these hiPSC-CMs. Gene expression data generated by qPCR or RNA-seq in particular focusing CM structure and function genes would be helpful too.

As suggested by referee 2, we have now performed gene expression analysis and immunostaining of PAX9 KD which gave the strongest phenotype in iPSC-CM (Figure 4 J-M). This unraveled increased expression of Na+ and K+ channels, which is in line with APD shortening phenotype, as well as down regulation of CASQ2, consistent with calcium transient shortening. Expression analysis also revealed increased sarcomeric genes and NPPA/B expression, which was consistent with increased CM size as quantified by the area of TNNT2 staining per nuclei.

These new data are described at the end of the result section:

“APD shortening for PAX9 KD was coincident with increased expression of Na+ and K+ ion channels (SCN5A, KCNH2 and KNCQ1) (Figure 4J), supporting the APD shortening phenotype. In this context, the AP kinetics also correlated with shorter calcium transient duration (Figure S8A-D and H-K), including faster upstroke and downstroke calcium kinetics and increased beat rate (peak frequency) (Figure S8E-G and L, M), consistent with decreased expression of Calsequestrin 2 isoform (CASQ2) associated with PAX9 KD (Figure 4J). Finally, assessment of the PAX9 KD effect on sarcomeric content revealed an increase in sarcomeric gene expression (Figure 4K), and an upregulation of genes associated with an hypertrophic response (NPPA, NPPB and NPR1) (Battistoni Et al. Circulating biomarkers with preventive, diagnostic and prognostic implications in cardiovascular diseases, Int J Cardiol, 2012, vol. 157) which was coincident with increased CM size as quantified by the area of TNNT2 staining per cardiac nuclei (Figure 4 L, M).

Collectively, these data illustrate conserved functions for poxm/PAX9 and sr/EGR2 in setting the cardiac rhythm and identify PAX9 as a novel and key regulator of cardiac performance at the cellular level, via the integrated regulation of expression of genes controlling electrophysiology, calcium handling and sarcomeric functions in hiPSC-CMs.”

Reviewer 3

Major Comments:

1. There is an assumption in the use of RNAi knockdown to validate the genes identified in the quantitative analysis, and that is that natural variants are themselves hypomorphic. It is possible that among the variants identified some are hypermorphic, or among the transcription factor binding sites that variants lead to increased factor binding. While RNAi knockdown is an excellent choice to begin validation, I do not think the authors can rule out that a gene not functionally validated by their RNAi tests does not have a role in cardiac function.

Please see our answers to reviewer 1 comment n°1 and reviewer 2 comment n°2.

2. After performing RNAi knockdown to validate genes identified by GWAS the authors focus on the TGFbeta signaling pathway for downstream analysis. To do so they examine heterozygotes for sog, a repressor of BMP signaling, and snoo, an activator of Activin pathway. The data from the snoo/sog heterozygote is compelling in its disruption of heart phenotypes, and the authors conclude a "coordinated action of activin and BMP." snoo, however, also works as a transcriptional repressor in the BMP pathway, so it's possible that the effects the authors are seeing here could be confined to an increase in BMP signaling. Unlike snoo and sog, mutations in babo and dpp are both expected to have negative effects on Activin and BMP signaling, respectively. The babo/dpp interaction is not as quantitatively convincing as the snoo/sog data, despite the integral roles both babo and dpp play in their respective pathways. If both pathways are connected, why do snoo/sog heterozygotes affect SI phenotypes, while babo/dpp heterozygotes affect fractional shortening? I think the authors data suggest an interesting potential interaction between these pathways, which could be confirmed by examining further mutant combinations, knockdowns or increased expression transgenes, but falls short of a "confirmed synergystic genetic interaction." It does, however, underscore the value of the data in the paper for opening up new avenues for future study.

(and reviewer 1 comments 3 and 5).

These comments led us to reconsider the analysis of the phenotypes associated with loss of function of the TGFb pathway, and to analyze other pathway components combinations.

We acknowledge reviewer 3 criticisms on snoo/sog experiments, which are difficult to interpret given the broad action snoo may have on both BMP and activin pathways. We have addressed this in the result section.

We have also analyzed other allelic combinations of BMP and activin pathways components, which strengthen the analysis performed on dpp/babo. Indeed, we tested babo/tkv heterozygotes (respectively specific activin and BMP receptors) and found significant genetic interactions for ESD and EDD. Albeit non-significant, babo/tkv double heterozygotes display a tendency to non-additive effects on FS (p= 0,054). mad/smox heterozygotes (respectively specific downstream TFs of BMP and activin pathways) display interactions (non-additive effects) on HP, SI, DI, ESD and EDD. These new results (Supplemental Figure 4) are thus supporting the hypothesis of genetic interactions between the pathways, but also reveal, as suggested by reviewer 3, a complex relationship between both pathways since interactions are revealed for specific traits in each of the mutant combinations analyzed.

The phenotypes related to the individual loss of function of each of the actors of these pathways (dpp, tkv and mad for BMP; babo and smox for activin) are however very similar. When they have an effect, heterozygous amorphic alleles of these genes display increased phenotypes related to rhythmicity (HP, DI, SI, AI) and FS, but decreased cardiac diameters (ESD and EDD).

Finally, as pointed out by reviewer 1, the picture is certainly even more complex since the phenotypes of RNAi mediated heart specific loss of function are not always similar to those of systemic loss of function. Indeed, mad RNAi causes a reduction of HP, DI, SI and FS (Figure S5) whereas heterozygotes for mad12 have either no or opposite effect on these phenotypes, and mad RNAi causes a significative increase in AI whereas mad12 has no effect (Figure S4). The discrepancy between tissue specific RNAi and heterozygous background was also found in the case of dpp, but specifically for the FS. Indeed, as suggested by reviewer 1 we have analyzed the loss of function of dpp by heart-specific RNAi. dpp RNAi results in a reduction of the FS (like mad RNAi) whereas the loss of function in the whole-body results in an increase of the FS.

We therefore re-wrote the whole corresponding section of the results and modified Figure S4 to include babo/tkv; smox/mad and dpp RNAi data.

“We further focused on the TGFb pathway, since members of both BMP and activin pathways were identified in our analyses. We tested different members of the TGFb pathway for cardiac phenotypes using cardiac specific RNAi knockdown (Figure 2C), and confirmed the involvement of the activin agonist snoo (Ski orthologue) and the BMP antagonist sog (chordin orthologue). Notably, Activin and BMP pathways are usually antagonistic (Figure 2D). Their joint identification in our GWAS suggest that they act in a coordinated fashion to regulate heart function. Alternatively, it may simply reflect their involvement in different aspects of cardiac development and/or functional maturation. In order to discriminate between these two hypotheses, we tested if different components of these pathways interacted genetically. Single heterozygotes for loss of function alleles show dosage-dependent effects of snoo and sog on several phenotypes, providing an independent confirmation of their involvement in several cardiac traits (Figure S4). Importantly, compared to each single heterozygotes, snooBSC234/ sogU2 double heterozygotes flies showed non additive SI phenotypes (two-way ANOVA p val: 2,1 10-7) suggesting a genetic interaction (Figure 2E and Figure S4A). It is worth noting however that snoo is also a transcriptional repressor of the BMP pathway (PMID: 16951053). The effect observed in snooBSC234/ sogU2 double heterozygotes can therefore alternatively arise as a consequence of an increased BMP signaling without affecting the activin pathway. We thus tested other allelic combinations for loss of function alleles of BMP and activin pathways. babo/tkv heterozygotes (respectively activin and BMP type 1 receptors) displayed non additive ESD and EDD phenotypes (Figure S4C). Synergistic interaction of BMP and activin pathways was also suggested by the analysis of fractional shortening in loss of function mutants for babo and dpp, the BMP ligand (Figure S4B). Of note, babo/tkv double heterozygotes also displayed a tendency to non-additive effects on FS albeit non-significant (two-way anova p= 0,054). In addition, mad/smox heterozygotes (specifc downstream TFs of BMP and activin pathways) displayed non-additive effects on several traits, including phenotypes related to rhythmicity (HP, SI, DI) and contractility (ESD and EDD) (Figure S4D). Altogether, cardiac performance in response to allelic combinations of activin and BMP supported a coordinated action of both pathways in the establishment and/or maintenance of cardiac activity. This was further supported by the observation that simple heterozygotes for the tested loss of function alleles displayed similar trends with respect to cardiac performance, irrespective of the pathway considered (dpp, tkv and mad for BMP; babo and smox for activin). Indeed, they displayed either no effect or increased fractional shortening and rhythmicity phenotypes (HP, DI, SI, AI), and decreased cardiac diameters (ESD and EDD). This suggests coordinated activity of both pathways. Importantly, the genetic interactions were tested using amorphic alleles that lead to systemic loss of function. The observed phenotypes may thus not unravel cardiac specific effects of the pathways. In support of this, mad cardiac specific RNAi knock down was tested (see below, Figure S5) and lead to a decreased HP, DI, SI and FS whereas heterozygotes for mad12 have either no (FS) or opposite (HP, DI, SI) effect on these phenotypes (Figure S4D). Inversely, mad RNAi caused a significant increase in AI whereas mad12 had no effect. However, heart specific dpp RNAi knock down (Figure S4E) lead to similar phenotypic trends compared to dppd14 (increased HP, DI, SI, decreased EDD and ESD) with the notable exception of FS which was reduced following cardiac specific KD (Figure S4E), but increased in dppd14 heterozygotes (Figure S4B). Taken together, these data point to a complex picture of TGFb pathway activity in regulating cardiac performance, involving both the activin and the BMP pathways as well as gene specific effects with both systemic and tissue-specific contributions.”

Minor Comments:

There is an enormous amount of data in this paper, but there are places where things are summarized a little too briefly. For example, there are no definitions given at the beginning of the Results section for traits like "Heart Period" or "Systolic Interval," which would make this work significantly more accessible for other Drosophila researchers. (They do touch on this when they explain later in the paper that certain variants are "associated with quantitative traits linked to heart size and contractility" but more background earlier would be helpful.) When we consider heart performance traits, what is the baseline from known mutants? In other words, where is the line between variation and defect?

– We have detailed the description of the traits analyzed at the beginning of the result section. We hope this improves the ease of reading in the direction suggested by the reviewer.

“7 cardiac traits were analyzed across the whole population (Dataset S1 and Table 1). As illustrated in Figure 1A, we analyzed phenotypes related to the rhythmicity of cardiac function: the systolic interval (SI) is the time elapsed between the beginning and the end of one contraction, the diastolic interval (DI) is the time elapsed between two contractions and the heart period (HP) is the duration of a total cycle (contraction + relaxation (DI+SI)). The arrhythmia index (AI, std-dev(HP)/mean (HP)) is used to evaluate the variability of the cardiac rhythm. In addition, 3 traits related to contractility were measured. The diameters of the heart in diastole (End Diastolic Diameter, EDD), in systole (End Systolic Diameter, ESD), and the Fractional Shortening (FS), which measures the contraction efficacy (EDD-ESD/EDD).”

– With respect to the baseline of cardiac performance, there is no simple answer. The baseline is influenced by the genetic background and the experimental conditions. This is the reason why any analysis of mutants or RNAi is conducted in comparison with its own control, analyzed at the same time. Concerning the DGRP lines, no baseline can be defined, since the objective is to measure the diversity of cardiac performance traits within a natural population.

https://doi.org/10.7554/eLife.82459.sa2

Article and author information

Author details

  1. Saswati Saha

    Aix-Marseille University, INSERM, TAGC, Turing Center for Living systems, Marseille, France
    Contribution
    Data curation, Formal analysis, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Lionel Spinelli
    Competing interests
    No competing interests declared
  2. Lionel Spinelli

    Aix-Marseille University, INSERM, TAGC, Turing Center for Living systems, Marseille, France
    Contribution
    Conceptualization, Formal analysis, Investigation, Visualization, Writing - review and editing
    Contributed equally with
    Saswati Saha
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9228-8141
  3. Jaime A Castro Mondragon

    Centre for Molecular Medicine Norway (NCMM), Oslo, Norway
    Contribution
    Formal analysis, Investigation, Methodology
    Competing interests
    No competing interests declared
  4. Anaïs Kervadec

    Development, Aging and Regeneration Program, Sanford Burnham Prebys Medical Discovery Institute, La Jolla, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  5. Michaela Lynott

    Development, Aging and Regeneration Program, Sanford Burnham Prebys Medical Discovery Institute, La Jolla, United States
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  6. Laurent Kremmer

    Aix-Marseille University, INSERM, TAGC, Turing Center for Living systems, Marseille, France
    Present address
    INRAE, Institut Sophia Agrobiotech, Université Côte d’Azur, CNRS, Sophia Antipolis, France
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  7. Laurence Roder

    Aix-Marseille University, INSERM, TAGC, Turing Center for Living systems, Marseille, France
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  8. Sallouha Krifa

    Aix-Marseille University, INSERM, TAGC, Turing Center for Living systems, Marseille, France
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  9. Magali Torres

    Aix-Marseille University, INSERM, TAGC, Turing Center for Living systems, Marseille, France
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  10. Christine Brun

    1. Aix-Marseille University, INSERM, TAGC, Turing Center for Living systems, Marseille, France
    2. CNRS, Marseille, France
    Contribution
    Resources, Data curation, Methodology
    Competing interests
    No competing interests declared
  11. Georg Vogler

    Development, Aging and Regeneration Program, Sanford Burnham Prebys Medical Discovery Institute, La Jolla, United States
    Contribution
    Resources, Software
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8303-3531
  12. Rolf Bodmer

    Development, Aging and Regeneration Program, Sanford Burnham Prebys Medical Discovery Institute, La Jolla, United States
    Contribution
    Conceptualization, Funding acquisition
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9087-1210
  13. Alexandre R Colas

    Development, Aging and Regeneration Program, Sanford Burnham Prebys Medical Discovery Institute, La Jolla, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Writing - review and editing
    For correspondence
    acolas@sbpdiscovery.org
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8489-0570
  14. Karen Ocorr

    Development, Aging and Regeneration Program, Sanford Burnham Prebys Medical Discovery Institute, La Jolla, United States
    Contribution
    Conceptualization, Supervision, Funding acquisition, Writing - review and editing
    For correspondence
    kocorr@SBPdiscovery.org
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2593-0119
  15. Laurent Perrin

    1. Aix-Marseille University, INSERM, TAGC, Turing Center for Living systems, Marseille, France
    2. CNRS, Marseille, France
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Investigation, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    laurent.perrin@univ-amu.fr
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4139-8743

Funding

Fondation de France (00071034)

  • Laurent Perrin

Aix-Marseille Université (A*MIDEX)

  • Laurent Perrin

National Institutes of Health (R01 HL148827-03)

  • Alexandre R Colas

National Institutes of Health (R01 HL149992-02)

  • Alexandre R Colas

National Institutes of Health (R01 AG071464-01)

  • Alexandre R Colas

U.S. Department of Defense (R01 AG071464-01)

  • Alexandre R Colas

National Institutes of Health (R01 HL132241)

  • Karen Ocorr

National Institutes of Health (R21 AG061598)

  • Karen Ocorr

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The project leading to this publication has received funding from Fondation de France (N°00071034) and from Excellence Initiative of Aix-Marseille University - A*MIDEX (A*Midex International) to LP, from the National Institute of Health (R01 HL148827-03, R01 HL149992-02, R01 AG071464-01) and from the Department of Defence (W81XWH-21-1-0104-01) to AC and from National Institute of Health (R01 HL132241 and R21 AG061598) to KO. We thank the Bloomington Drosophila Stock Center (BDSC) and the Vienna Drosophila Research Center (VDRC) for fly stocks. Centre de Calcul Intensif d’Aix-Marseille is acknowledged for granting access to its high-performance computing resources.

Senior and Reviewing Editor

  1. Detlef Weigel, Max Planck Institute for Biology Tübingen, Germany

Publication history

  1. Preprint posted: June 8, 2021 (view preprint)
  2. Received: August 4, 2022
  3. Accepted: October 25, 2022
  4. Version of Record published: November 16, 2022 (version 1)

Copyright

© 2022, Saha, Spinelli et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 212
    Page views
  • 30
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Saswati Saha
  2. Lionel Spinelli
  3. Jaime A Castro Mondragon
  4. Anaïs Kervadec
  5. Michaela Lynott
  6. Laurent Kremmer
  7. Laurence Roder
  8. Sallouha Krifa
  9. Magali Torres
  10. Christine Brun
  11. Georg Vogler
  12. Rolf Bodmer
  13. Alexandre R Colas
  14. Karen Ocorr
  15. Laurent Perrin
(2022)
Genetic architecture of natural variation of cardiac performance from flies to humans
eLife 11:e82459.
https://doi.org/10.7554/eLife.82459
  1. Further reading

Further reading

    1. Genetics and Genomics
    2. Plant Biology
    Simon Snoeck, Bradley W Abramson ... Adam D Steinbrenner
    Research Article Updated

    As a first step in innate immunity, pattern recognition receptors (PRRs) recognize the distinct pathogen and herbivore-associated molecular patterns and mediate activation of immune responses, but specific steps in the evolution of new PRR sensing functions are not well understood. We employed comparative genomic and functional analyses to define evolutionary events leading to the sensing of the herbivore-associated peptide inceptin (In11) by the PRR inceptin receptor (INR) in legume plant species. Existing and de novo genome assemblies revealed that the presence of a functional INR gene corresponded with ability to respond to In11 across ~53 million years (my) of evolution. In11 recognition is unique to the clade of Phaseoloid legumes, and only a single clade of INR homologs from Phaseoloids was functional in a heterologous model. The syntenic loci of several non-Phaseoloid outgroup species nonetheless contain non-functional INR-like homologs, suggesting that an ancestral gene insertion event and diversification preceded the evolution of a specific INR receptor function ~28 my ago. Chimeric and ancestrally reconstructed receptors indicated that 16 amino acid differences in the C1 leucine-rich repeat domain and C2 intervening motif mediate gain of In11 recognition. Thus, high PRR diversity was likely followed by a small number of mutations to expand innate immune recognition to a novel peptide elicitor. Analysis of INR evolution provides a model for functional diversification of other germline-encoded PRRs.

    1. Chromosomes and Gene Expression
    2. Genetics and Genomics
    Evgeniya N Andreyeva, Alexander V Emelyanov ... Dmitry V Fyodorov
    Research Article

    Asynchronous replication of chromosome domains during S phase is essential for eukaryotic genome function, but the mechanisms establishing which domains replicate early versus late in different cell types remain incompletely understood. Intercalary heterochromatin domains replicate very late in both diploid chromosomes of dividing cells and in endoreplicating polytene chromosomes where they are also underrelicated. Drosophila SNF2-related factor SUUR imparts locus-specific underreplication of polytene chromosomes. SUUR negatively regulates DNA replication fork progression; however, its mechanism of action remains obscure. Here we developed a novel method termed MS-Enabled Rapid protein Complex Identification (MERCI) to isolate a stable stoichiometric native complex SUMM4 that comprises SUUR and a chromatin boundary protein Mod(Mdg4)-67.2. Mod(Mdg4) stimulates SUUR ATPase activity and is required for a normal spatiotemporal distribution of SUUR in vivo. SUUR and Mod(Mdg4)-67.2 together mediate the activities of gypsy insulator that prevent certain enhancer-promoter interactions and establish euchromatin-heterochromatin barriers in the genome. Furthermore, SuUR or mod(mdg4) mutations reverse underreplication of intercalary heterochromatin. Thus, SUMM4 can impart late replication of intercalary heterochromatin by attenuating the progression of replication forks through euchromatin/heterochromatin boundaries. Our findings implicate a SNF2 family ATP-dependent motor protein SUUR in the insulator function, reveal that DNA replication can be delayed by a chromatin barrier and uncover a critical role for architectural proteins in replication control. They suggest a mechanism for the establishment of late replication that does not depend on an asynchronous firing of late replication origins.