Novel cardiovascular gene functions revealed via systematic phenotype prediction in zebrafish

De novo gene function discovery has been greatly facilitated by systematic gene deletion and observation of resulting phenotypes in scalable model organisms. Indeed, systematic gene disruptions in S. cerevisiae (Costanzo et al., 2010; Giaever et al., 2002), C. elegans (Kamath et al., 2003) and Drosophila (Boutros et al., 2004) have each revealed molecular functions for thousands of genes. However, given the breadth and complexity of observable phenotypes in vertebrates, comprehensive assessment of gene function through serial observation of all possible phenotypes following gene disruption remains infeasible. A more efficient alternative would be to use gene function prediction to prioritize gene candidates for more detailed phenotypic testing based on a variety of known gene and protein properties and relationships.

Computational prediction of molecular function has been effective in assigning physiological roles to genes across eukaryotic model organisms (Deng et al., 2004; Guan et al., 2008; Huttenhower et al., 2009; Karaoz et al., 2004; Lee et al., 2004; Mostafavi et al., 2008; Tasan et al., 2012; Tasan et al., 2008; Troyanskaya et al., 2003). Similar prediction frameworks have been applied to predict associated phenotypes in yeast (King et al., 2003; Saha and Heber, 2006) and worm (Lee et al., 2008) and to identify putative human disease gene candidates (Linghu et al., 2009; Woods et al., 2013) but have not been systematically coupled with in vivo validation in a vertebrate model organism.

Predictions of gene function or phenotype have generally used at least one of two basic strategies. The first strategy, termed guilt-by-profiling (GBP), begins by identifying gene properties (features, e.g. ‘gene is expressed in the brain’) that are common to the genes currently associated with a particular function of interest (training set, e.g. ‘genes known to be important for cognition’). Additional, untested genes (test set) exhibiting the hallmarks of the function of interest are then predicted to have that function. A second general strategy, termed guilt-by-association (GBA), identifies gene-gene relationships (e.g. a physical interaction between the corresponding proteins) that tend to connect genes that share a function. Untested genes that are ‘well-connected’ to the set of genes known to hold a particular function are then prioritized as top candidates.

Although there are many effective prediction algorithms, a recent benchmarking comparison (Peña-Castillo et al., 2008) ranked the Funckenstein approach (Tian et al., 2008) highly, based on estimates of precision of top-ranked function predictions. Funckenstein uses both GBA and GBP, combining results to output a single confidence score for each potential gene-function or gene-phenotype association. These predicted relationships can then be directly validated experimentally.

Zebrafish embryos are transparent, small, fast growing and are developmentally similar to higher vertebrates (Howe et al., 2013), making them ideal for large-scale study of vertebrate developmental gene functions. Zebrafish are amenable to forward genetics via chemical (Driever et al., 1996; Haffter et al., 1996) or retroviral-based (Amsterdam et al., 1999; Golling et al., 2002) mutagenesis, and to reverse genetics via knockdown (Nasevicius and Ekker, 2000), small molecules (Burns et al., 2005; Peterson et al., 2000) or, more recently, targeted gene deletion (Huang et al., 2011; Hwang et al., 2013; Meng et al., 2008). Morpholinos, which are oligomers that block translation through steric transcript inhibition, are particularly efficient (Nasevicius and Ekker, 2000), offering scalable, specific gene product inhibition lasting for up to 72 hours.

Here we predicted the effects of gene knockdown on 338 zebrafish embryonic anatomical processes (Bradford et al., 2011), with the results of these predictions being made publicly available as a community resource through two searchable online gene browsers: FuncBase and GeneMANIA. We deliberately focused on the prediction of morphant rather than mutant phenotypes so as to allow direct confirmation of prediction results at a reasonable scale. Because extant data indicated high performance for predictions of cardiac phenotypes, we aimed to experimentally validate this framework for cardiovascular phenotype prediction. Expecting a broad range of potential cardiac phenotypes, we used a five-parameter manual evaluation system to assess cardiac function, and applied it in a blinded fashion to predictions and to positive and negative controls. Subsequent quantitative assessments of cardiac physiology largely confirmed the predictions, identifying novel myocardial effects for 11 surveyed genes. These 11 include tmem88a, a known attenuator of Wnt signaling (Lee et al., 2010), shown here to influence the physiological coupling of embryonic cardiomyocytes.

We implemented a machine learning procedure based on a diverse collection of large-scale datasets to predict phenotypes resulting from loss of gene function across the zebrafish genome (Fig. 1). We then used each of the generated computational models (one for each of the 338 anatomical process terms) to score the 15,103 zebrafish genes with available gene feature information. The over 5 million assigned scores each reflect the likelihood of observing a given phenotype (disruption of a given anatomical process) following loss of function of a specific gene.

Gene features most important for the prediction process differed for guilt-by-association (GBA) and guilt-by-profiling (GBP) predictors. For GBP, feature importance varied based on the phenotype being predicted, but tissue expression and phylogenetic relationships appeared to consistently be the most relied upon features (see supplementary material Table S1 for gene feature rankings). In GBA, where only a single random forest classifier was used to create a gene-gene functional linkage network, phylogenetic profile similarity was the most important feature, followed by similarity in GO terms (especially those describing cellular compartments), and then similarity in protein family IDs. These feature importance results are consistent with the observed usefulness of phylogenetic profiles in predicting gene function and phenotype in other species (Levesque et al., 2003; Pellegrini et al., 1999).

Examining the training data for the 338 phenotype predictors, the broadest phenotypic terms assigned to the most genes (e.g. anatomical system, whole organism) had poor performance, either because of the limited number of negative training examples or because the predictive features we used simply held little value in predicting function at this level of abstraction. Predictions also appeared highly variable when only a small number of positive training examples were present (see supplementary material Table S2 for a full list of prediction scores by phenotype). Filtering out anatomical process terms on both extremes of the specificity spectrum (supplementary material Table S2) resulted in 242 retained terms.

Examining remaining predictions by anatomical process term, the prediction process appeared to be most effective for neuronal, sensory and cardiac phenotypes (Fig. 2A), and performed most poorly for hematopoietic phenotypes (supplementary material Table S2). Additionally, within each phenotypic category there was substantial variation in the precision of predictors. For example, among cardiovascular phenotypes (Fig. 2B), predictors of cardiac structures (e.g. heart tube, cardiac muscle cell, myocardium) showed greater performance than terms abstractly describing vascular structure and the endocardium. This did not appear to be a ‘rich-getting-richer’ scenario of terms with many training examples showing greater performance, as the top five scoring cardiovascular phenotype terms had 15 or fewer training examples. We hypothesize instead that performance stems from the fidelity of observation of anatomical structures in the transparent embryo, which might have allowed these defects to be more rigorously and consistently characterized. For example, in the case of hematopoietic phenotypes, abnormal blood flow or blood accumulation is often reported as a secondary consequence of another phenotype, which might have led to inconsistency among associated genes and poor performance for the term ‘blood’. Conversely, the term ‘nucleate erythrocyte’, which is more stringently defined, showed better performance (supplementary material Table S2).

Examining individual gene-phenotype predictions in the test set, 8742 gene/term predictions involving 2459 genes and 93 terms were above an 80% estimated precision cut-off. These high-confidence predictions (see supplementary material Table S3) are searchable online in the context of the FuncBase Genome Browser [http://zfunc.mshri.on.ca (Beaver et al., 2010)] and the GeneMANIA functional annotation browser [http://www.genemania.org (Mostafavi et al., 2008)]. The majority (87%) of these 2459 genes were predicted to affect five or fewer anatomical processes at this precision level, with genes encoding ribosomal constituents frequently being most pleiotropic. Indeed, using gene set enrichment analysis (Subramanian et al., 2005) we found that genes annotated as being ribosomal were significantly enriched for participation in a greater number of phenotypes (q<0.001). This is consistent with the observed phenotypic importance of ribosomal proteins in other organisms (Giaever et al., 2002), as well as in prior genetic screens in zebrafish (Amsterdam et al., 2004).

As substantial literature support exists for cardiovascular phenotypes, we examined the literature surrounding our cardiovascular predictions in more detail. The phenotype term broadly describing the cardiovascular system exhibited a precision at 20% recall value of 0.72 (Fig. 2C). Examination of independent lines of evidence for top hits also suggested prediction results for ‘cardiovascular system’ to be of high quality. Indeed, six of the ten top-scoring predictions were already known to be associated with cardiovascular defects. These genes had not been used as positive training examples either because the phenotype was observed through genetic mutation rather than morpholino, because a morpholino had been injected into a zebrafish line not used for training purposes, or because the publications were curated in the Zebrafish Information Network (ZFIN) database (Bradford et al., 2011) subsequent to the download of phenotypic data to establish a training set. Examples of literature-validated predictions include vmhc, in which a premature stop codon causes ventricular failure (Auman et al., 2007), and apobec2a, which causes decreased heart rate upon knockdown (Etard et al., 2010).

To more systematically determine whether the predictive value of our framework extends beyond morpholino-mediated gene knockdown, we obtained all mutant-phenotype association data from ZFIN (using an identical approach to that used above for morpholino-phenotype data). In total, 818 genes in our test set had associated mutants with phenotype information. We used our ‘cardiovascular system’ phenotype predictor, trained exclusively using morpholino experimental data, to rank these genes based on likelihood to elicit a cardiovascular phenotype. Of the 50 highest ranked genes, 41 had associated mutants with a cardiovascular phenotype (P<1×10^-9, Fisher’s exact test).

To more directly test the results of our cardiovascular phenotype predictions, we selected all genes exceeding an estimated 95% precision cut-off for targeted morpholino knockdown, excluding genes for which morpholinos targeting the transcription start site could not be designed (as a result of high GC content or self-complementarity). To increase the potential for novel discovery, we also excluded genes that had existing deletion mutants or morpholino experimentation confirming cardiac function, or which, despite not having experimental evidence, had an apparent cardiac function (e.g. vmhcl). This resulted in a list of 16 target genes (Table 1). To determine how well our procedure could discriminate against genes unlikely to cause a phenotype, the five genes with the lowest prediction score were screened as well. Finally, a morpholino targeting a gene with a known severe cardiac developmental phenotype [slc8a1a, which encodes a sodium calcium exchanger (Langenbacher et al., 2005; Stainier et al., 1996)] was screened as a positive control (see supplementary material Table S4 for a list of all morpholinos used).

We used a simple five-parameter categorical system to score cardiac effects (‘S score’ hereafter; see Materials and methods and supplementary material Fig. S1), comparing S scores for all five cardiac parameters between test genes and negative controls. We reasoned that a manual, blinded evaluation of cardiac function would be more sensitive in identifying a range of potential cardiac defects than other comparable quantitative approaches. Using a cut-off value of S>0.5 to identify cardiovascular defects, 12 of 16 test genes were found to have a defect in at least one quantified parameter of heart function following knockdown, as compared with zero of five negative control genes (P<0.05, Fisher’s exact test; Table 1). Our results were not sensitive to the S score threshold used, in that test genes scored significantly higher than negative controls at a range of cut-offs spanning 0.3 to 0.8 (supplementary material Fig. S2A). Additionally, all five parameters of heart function scored significantly higher in test genes than negative controls (P<0.05; supplementary material Fig. S2B). Together, these data confirm that a machine learning procedure can be successfully used to differentiate genes for which loss of function elicits a cardiovascular phenotype.

One area of experimental concern given our deliberate focus on morphants was whether the cardiac effects observed could be due to morpholino toxicity or non-specific binding effects. Morpholino toxicity occurs randomly for reasons thought to be related to p53-mediated activation of apoptosis (Robu et al., 2007) and can induce a phenotype with a cardiovascular component. Additionally, although morpholino sequences were vetted for binding fidelity, off-target effects remain a possibility. The phenotypes observed appeared diverse with respect to the cardiac parameters affected, suggesting that these responses were not due to a shared general morpholino effect (supplementary material Fig. S2C). Based on the lack of phenotype observed in any of our five negative control morpholinos (each targeting a characterized gene and thus potentially subject to both morpholino toxicity and non-specific binding effects), we estimated the prevalence of cardiac phenotypes due to non-specific effects to be below 20%, which is in keeping with prior work (Ekker and Larson, 2001; Heasman, 2002; Nasevicius and Ekker, 2000). However, even if three of 12 hits (25%) were due to non-specific effects, highly ranked test genes would remain significantly differentiated from low-ranked controls (P<0.05, Fisher’s exact test).

We next used more detailed physiological measures to further characterize observed cardiac defects. We first examined alterations in cardiac output (CO) following injection of all test morpholinos. Here we followed a previously described approach to further confirm that observed defects were not due to non-specific toxicity: all CO experimentation was performed following co-injection of the test morpholino with a morpholino targeting p53 (also known as tp53) (Robu et al., 2007).

Of the 12 test morpholinos identified through the S score (S>0.5) as inducing a cardiac phenotype, eight caused a significant reduction in CO (supplementary material Table S5). Among these eight was tnni1b, which has known cardiac expression (Fu et al., 2009) and high similarity to human troponin I type 1 (skeletal, slow), but no known phenotype. Knockdown of tnni1b caused virtually complete asystole (supplementary material Movies 1, 2) with no obvious effects on zebrafish motility, a phenotype similar to that observed upon knockdown of cardiac troponin T (tnnt2a), the gene underlying the contractile defect in the silent heart mutant (Sehnert et al., 2002). A second morpholino targeting tnni1b splicing also demonstrated near asystole (supplementary material Movie 3) with no apparent motility defects, supporting this as a gene-specific effect.

Also among the set of genes causing reduced CO was hspb7, the human ortholog of which has been associated with heart failure through large-scale genome-wide association screening (Stark et al., 2010; Villard et al., 2011). A reduction in CO was confirmed using an additional morpholino targeting an mRNA splice site in the hspb7 transcript. However, unlike tnni1b, the phenotype following hspb7 knockdown did not appear exclusively cardiac. Most obviously, at the highest dose nearly half of all splice morpholino-injected embryos had a shortened tail, abridged at the trunk region (a phenotype not seen with any other morpholino injection; supplementary material Fig. S3). Supporting the generality of the systematic gene prioritization strategy, hspb7, in addition to its predicted cardiac effects, was associated with the phenotypic term ‘trunk’ at an estimated precision of 0.89.

Although they were confirmed as having a cardiac phenotype according to the S score, knockdown of nppa, itpr3, tmem88a and adprhl1 did not cause a significant reduction in CO. However, valvular regurgitation was observed in 55% of nppa morphant embryos examined at the highest injected dose (supplementary material Movie 4), as compared with 3% of wild-type embryos (as assessed through blinded manual evaluation; P<1×10^-5, Fisher’s exact test). A regurgitation phenotype was confirmed using a second morpholino targeting the nppa transcription start site through an alternative sequence (31% of injected embryos versus 4% in controls, P<0.05).

Next, itpr3 (inositol 1,4,5-triphosphate receptor, type 3) knockdown caused a notable alteration in atrial contraction, without a corresponding significant reduction in CO (CO is largely a metric of ventricular function). Because itpr3 is known to mediate the release of intracellular calcium in epithelial cells (Maranto, 1994) we went on to examine calcium dynamics in isolated hearts following itpr3 knockdown. We found a significant decrease in diastolic calcium concentration in the atrium of morphants following injection of both the initial itpr3 morpholino and a second morpholino targeting an itpr3 mRNA splice site (Fig. 3).

Finally, unlike nppa, itpr3 and tmem88a (tmem88a knockdown caused quantifiable alterations in cardiac conduction; more detail below), the adprhl1 knockdown phenotype was largely abrogated by co-injection with p53 morpholino. This suggests that the initially observed adprhl1 phenotype is either p53 dependent or is due to morpholino toxicity.

Thus, of 12 test genes that were both computationally predicted and subsequently identified through the S score as having a cardiac phenotype, 11 were further validated to have underlying alterations in atrial function, valvular function, CO or cardiac conduction.

The transmembrane protein Tmem88 is a recently reported inhibitor of the Wnt signaling pathway in multiple species, acting through interaction with Dishevelled (Dvl) via a C-terminal tri-peptide VWV interaction motif (Lee et al., 2010). Knockdown of the zebrafish tmem88 ortholog tmem88a caused observable alterations in heart function at the injected dose (Table 1), albeit with no significant corresponding change in CO. Similarly, mosaic overexpression of a tmem88a transcript lacking the VWV motif (see Materials and methods) resulted in severe, systematic growth defects for 17% of embryos and observable cardiac defects for an additional 20% of embryos. By contrast, mosaic expression of the full tmem88a transcript resulted in no observable phenotypes.

Given the known effect of the Tmem88 VWV motif on Wnt signaling (Lee et al., 2010), the described role of wnt11 in establishing physiological myocardial electrical polarities in the developing zebrafish ventricle (Panáková et al., 2010), and the apparent cardiac defects following tmem88a knockdown, we next sought to identify any functional synergy between wnt11 and tmem88a. Co-injection of either the tmem88a ATG morpholino tested above or a second morpholino blocking tmem88a splicing [both independently validated as targeting tmem88a (Cannon et al., 2013)] with a previously validated wnt11 morpholino (Panáková et al., 2010) showed sensitization of the wnt11 cyclopia phenotype (Fig. 4). This suggests a synergistic genetic interaction between tmem88a and wnt11. Quantifying cell coupling, we found an increase in impulse propagation velocity in the body of the ventricle following tmem88a knockdown by both the ATG and splice-blocking morpholinos (Fig. 5A-D). This finding is consistent with an effect that enhances the intercellular coupling gradient induced by physiological levels of Wnt11 (Panáková et al., 2010). We found no appreciable difference in epithelial morphogenesis or ventricular connexin 43 expression following tmem88a inhibition that could otherwise account for this alteration in impulse propagation (Fig. 5E-L).

We next examined the expression of several key developmental markers following tmem88a knockdown to confirm that the observed effect of tmem88a on cardiac cell coupling is not a secondary consequence of a generalized effect on development (Fig. 6). fli1 served as an in situ hybridization marker for the integrity of endothelial development, myoD (also known as myod1) enabled monitoring of muscle formation, and beta-globin E3 (hbbe3) provided a readout for functional circulation and timing of hematopoiesis, as its expression is reduced by 48 hours postfertilization (hpf). Expression of fli1 and myoD in tmem88a morphant embryos was comparable to that in the wild-type reference. beta-globin E3 showed an apparent upregulation at 48 hpf in morphants and a slightly longer persistence of expression at 72 hpf when voltage measurements were taken (Fig. 6). These results suggest that tmem88a morpholinos might induce alterations in hematopoiesis at the injected dose [consistent with recent findings (Cannon et al., 2013)], but do not cause any obvious systematic developmental abnormalities that might otherwise account for the cardiac effects observed at the time point surveyed.

In this study we used high-confidence predictions of gene-phenotype association across over 15,000 genes to prioritize the screening of novel gene candidates for developmental functions. Through rigorous experimental assay we identified genes affecting a range of observable cardiac phenotypes. This approach to large-scale phenotype prediction and in vivo validation in a vertebrate establishes a framework for de novo discovery of gene functions across a broad spectrum of vertebrate developmental phenotypes.

Using machine learning, we identified patterns in gene expression, conservation and interaction to generate over 5 million predictions of gene-phenotype association in embryonic zebrafish. We have made the resulting high-confidence predictions of gene-phenotype association available through two public browsers: FuncBase and GeneMANIA. GeneMANIA is capable of integrating existing predictions with new gene feature data as they become available, thereby acting as a dynamic resource for future experimentation. The precision of phenotype predictions appeared highest for consistently defined anatomical phenotypes, suggesting that effective gene prioritization can be achieved for any standardized phenotype with sufficient rigorous training information.

We used morpholino knockdown to both train and test our computational approach. Morpholinos enable gene knockdown at a reasonable scale, but can have appreciable phenotypic noise due to non-specific toxicity. Despite these confounders, a clear phenotypic enrichment was seen using our in silico prediction strategy.

To directly test the fidelity of our gene-phenotype association predictions, we chose to focus on the observation of cardiovascular phenotypes. From modeling basic processes such as cellular migration (Lazic and Scott, 2011; Zhou et al., 2011) and intercellular coupling (Panáková et al., 2010), to complex arrhythmogenic (Arnaout et al., 2007) and structural (Vogel et al., 2009) cardiac disorders, zebrafish are a powerful model organism for the study of cardiac development. Zebrafish have over 1000 genes with detectable cardiac expression [as annotated by ZFIN (Bradford et al., 2011)] and thus systematic screening of each gene for associated phenotypes would require tremendous experimental resources. Alternatively, our method allowed focus on a set of genes enriched for potential phenotypes.

Through testing and validation, we were able to identify 11 genes with cardiac phenotypes. Notably, despite the fact that cardiac specificity of predicted phenotypes was not included as a criterion for initial survey, the majority of these 11 genes had an apparently cardiac-specific phenotype upon knockdown (supplementary material Table S4). Eight of these genes were demonstrated to reduce CO and the remaining three (nppa, itpr3 and tmem88a) were shown to alter other quantifiable aspects of cardiac function.

nppa, which was known to have variant expression associated with alteration in valve tissue formation in zebrafish (Camarata et al., 2010), was shown here to be associated with regurgitant defects in vivo. Additionally, itpr3, recently shown to be dysregulated in human arrhythmia (Hasdemir et al., 2010), was observed to have a direct role in atrial calcium handling. Last, tmem88a, a proposed attenuator of Wnt signaling (Lee et al., 2010), was demonstrated here to regulate ventricular cell coupling.

We observed the influence of tmem88a on wnt11 in two separate, quantifiable phenotypes (cyclopia and myocardial cell coupling), suggesting a functional synergy between these two genes. Both wnt11 overexpression and knockdown are associated with a disruption of the myocardial electrical gradient (Panáková et al., 2010). tmem88a knockdown enhanced this gradient, suggesting that the relationship between tmem88a and wnt11, potentially mediated by the known interaction between tmem88a and dvl (Lee et al., 2010), appears more complex than a simple activation or inactivation. Notably, two additional studies have been published during the review of our manuscript (Novikov and Evans, 2013; Palpant et al., 2013), both confirming a Wnt-dependent role for tmem88a in cardiomyocyte differentiation.

The number of vertebrate phenotypes potentially observable even within the finite timeframe of embryonic development is vast. When combined with the enormous number of potential genotypes in a genome, it is clear that comprehensive experimental mapping from genotype to phenotype is unlikely without informed direction of experimental effort. Here we demonstrate that systematic data integration and objective gene prioritization can successfully direct limited experimental resources for high-resolution phenotyping to ranked subsets of genes for each phenotype, accelerating the potential for gene function discovery in vertebrates. With large-scale projects underway to generate and phenotype zebrafish carrying null alleles (Kettleborough et al., 2013) and other allelic classes, future iterations of our approach will be crucial in deciphering the complex relationship between genotype and phenotype.

We obtained gene-phenotype associations to train our learning procedure from the Zebrafish Information Network (ZFIN) database (Bradford et al., 2011) and included only experimentally verified gene-phenotype associations that satisfied the following two criteria: (1) gene inhibition was via a morpholino that had only one identifiable target annotated by ZFIN; and (2) the morpholino causing the observed phenotype was injected into a wild-type (i.e. non-transgenic) fish line, as we anticipated carrying out experimental verification using wild-type lines. Phenotypes used in the learning procedure were based on a hierarchical list of anatomical terms developed by a consortium of researchers together with ZFIN (Bradford et al., 2011). Through literature curation, ZFIN has captured the anatomical structures in embryonic development that are affected by injection of a given morpholino. Existing gene-phenotype associations were propagated to ancestral terms within the ontology prior to training.

Gene feature data used for training and prediction included: expression data [both tissue localization (in situ data collected by ZFIN) and microarray experimental data (more below)], phylogenetic profiles, annotated protein domains, Gene Ontology (GO) (Ashburner et al., 2000) and KEGG (Kanehisa et al., 2008) annotations, and genetic and protein interactions mapped through orthology from human, mouse and yeast. Only systematically obtained interaction data were mapped from additional model organisms so as to reduce the impact of ascertainment bias. As functional annotations can be derived directly from phenotypic evidence and thus introduce circularity into our prediction procedure, annotations derived from GO underwent stringent filtering before use. Expression data included both cell type-specific expression (annotated by ZFIN) and an experimental expression compendium generated for this purpose. To generate this compendium, we combined over 60 experimental datasets that had been published using the Affymetrix zebrafish platform from the Gene Expression Omnibus (Barrett and Edgar, 2006), after renormalizing each set using genechip robust multiarray averaging (Wu et al., 2004) via Bioconductor (Gentleman et al., 2004).

The phenotype prediction procedure used here was based on the Funckenstein approach previously used to make predictions of gene function in yeast (Tian et al., 2008), mouse (Tasan et al., 2008) and human (Tasan et al., 2012). Briefly, separate classifiers were used to infer putative gene-phenotype associations based on correlations between gene features and phenotypes (Guilt-By-Profiling; GBP), and to transfer known phenotype associations between genes (Guilt-By-Association; GBA) (Fig. 1). Both GBA and GBP used random forest classifiers for prediction. A random forest classifier is an ensemble-based learner in which collections of decisions trees are constructed, each using a randomly selected subset of training information and features.

For GBP, separate random forest classifiers were constructed for each ZFIN anatomical term, resulting in a gene by anatomical term matrix of confidence values. Positive training examples for each anatomical term were genes for which inhibition through morpholino affected that term (‘phenotype’ hereafter). Negative training genes were all other genes annotated with at least one phenotype term. We used out of bag (OOB) scores to evaluate prediction accuracy for each phenotype, such that each gene is assigned a score using only models that had not been trained using information about that gene.

For GBA, gene feature data were used to construct a functional linkage network (FLN), a graph in which gene-gene linkage weight corresponds to the strength of association between two genes, as estimated via a random forest classifier. Positive training examples were all genes with any shared phenotype, while negative examples were all other gene pairs drawn from the same set that failed to exhibit any shared phenotypes. OOB scores were used to generate performance estimates for gene-gene associations. Gene-phenotype associations were then derived from the FLN probabilistically. Specifically, edge weights within an FLN clique containing all genes annotated as causing a particular phenotype (herein referred to as the ‘core’ set) were used to derive a probability density function (PDF), an estimate of the distribution of edge weights expected of gene pairs sharing the phenotype of interest. A separate PDF was created using edge weights between the core set and all other genes not having the phenotype of interest. Edge weights between each candidate gene and the core set were then obtained, and a gene-phenotype association score was assigned to be the log of the likelihood of obtaining the observed edge-weight distribution under the core PDF model relative to the likelihood under the non-core PDF, as in previous approaches (Tasan et al., 2012). GBA and GBP scores were then combined for each term using a logistic regression model optimized to maximize the cross-validated estimates of area under the cumulative precision recall curve.

Gene feature data were used differently for GBA and GBP. For example, while phylogenetic data were used as a feature for both GBA and GBP, GBP contained all species used in phylogenetic profiling as individual features, while GBA used only the similarity in phylogenetic profile between any two genes. Similarly, correlation in expression over the 60 combined experimental datasets was used as an individual feature in GBA. However, in GBP, genes were clustered based on correlation in expression over these 60 datasets, and the presence in any of these expression clusters was used as an individual feature. In total, there were 65 features used to train the GBA predictor, while 3034 were used for GBP. Feature importance for both GBA and GBP was calculated as the root mean square error introduced by excluding the given feature.

Efficacy of each of the phenotype predictors was assessed according to the precision at 20% recall using cross-validation. This measure quantifies precision (the estimated fraction of predicted gene-phenotype associations that are true, i.e. the positive predictive value) at a given recall (the fraction of known gene-phenotype associations that were correctly predicted, i.e. the true positive rate or sensitivity). Precision-recall curves for each phenotype were generated using ROCR (Sing et al., 2005), and interpolated precision at 20% recall was calculated as previously described (Manning et al., 2008). Where the precision at 20% recall values were compared for multiple predictors, a relative precision score was calculated as the precision at 20% recall divided by the prior expectation (fraction positive training examples). This was done when comparing across multiple phenotypes to reduce the influence of differences in prevalence (i.e. the fraction of genes known to be associated with the phenotype). To generate the high confidence gene-phenotype association list, a likelihood score corresponding to a precision of 0.8 was determined for each phenotype, and all genes not in the original training data scoring above this threshold were included.

For GBA, feature data were used to construct a functional linkage network (FLN), a graph in which edge weight corresponds to the strength of association between two genes as estimated via a random forest classifier. This classifier was trained to predict edge weights based on known examples of shared phenotype (positive), versus all other gene pair combinations in the training set (negative), with all possible gene-gene pairings acting as a test set. To avoid the inflated performance estimates that can arise from overfitting, we used out of bag (OOB) scores (such that each gene is assigned a score using only models that had not been trained using information about that gene) to generate performance estimates for gene-gene associations. Gene-phenotype associations were then derived from the FLN probabilistically. Specifically, edge weights within an FLN clique containing all genes annotated as causing a particular phenotype (herein referred to as the ‘core’ set) were used to derive a probability density function (PDF), an estimate of the distribution of edge weights expected of gene pairs sharing the phenotype of interest. A separate PDF was created using edge weights between the core set and all other genes not having the phenotype of interest. Edge weights between each candidate gene and the core set were then obtained, and a gene-phenotype association score was assigned to be the log of the likelihood of obtaining the observed edge-weight distribution under the core PDF model relative to the likelihood under the non-core PDF, as in previous approaches (Tasan et al., 2012). Performance estimates for the gene-phenotype association predictions were evaluated using leave-one-out cross-validation.

For GBP, separate random forest classifiers were constructed for each phenotype, resulting in a gene by phenotype matrix of confidence values. OOB scores were used to evaluate prediction accuracy for each phenotype. GBA and GBP scores were then combined for each term using a logistic regression model optimized to maximize the cross-validated estimates of area under the cumulative precision versus recall curve.

ATG-blocking morpholinos were designed using Gene Tools oligo design service for uniformity. To ensure a lack of off-target binding, the resultant morpholino sequences were aligned against the zebrafish genome using BLAST (Altschul et al., 1990) with settings optimized for small sequences.

Male and female wild-type (AB) fish were housed and embryos bred for microinjection using standard protocols. Embryos were collected, pooled and immediately used for injection. Techniques used for morpholino injection follow those previously outlined (Westerfield, 2000). Briefly, morpholinos (Gene Tools) were resuspended in sterile water at 1 mM and diluted to working concentration with Danieau’s solution [58 mM NaCl, 0.7 mM KCl, 0.4 mM MgSO₄, 0.6 mM Ca(NO₃)₂, 5 mM HEPES]. Morpholinos were introduced into the zebrafish yolk via microinjection no later than the two-cell developmental stage. Injected embryos were then kept at 28.5°C in E3 solution (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl₂, 0.33 mM MgSO₄).

Each morpholino was injected into zebrafish embryos at a range of doses high enough to elicit an observable phenotype (each initially at three doses: 0.5, 0.25 and 0.125 pmol). If any morpholino caused frequent, systematic growth deformities, or resulted in greater than 30% mortality at the highest dose, a lower dose was used for further experimentation. Once a suitable morpholino dose was found, injections were performed in triplicate (at least 40 embryos per injection, three sets of embryos resulting from separate mate pairings).

Custom software (implemented in MATLAB; freely available upon request) was used to determine heart rate from sequential image files (obtained as above), while measurements of ventricular long and short axis in both diastole and systole were obtained manually for each video using ImageJ (http://rsbweb.nih.gov/ij/) and used to estimate chamber volume using standard geometric assumptions. Cardiac output was then calculated as diastolic minus systolic ventricular volume multiplied by heart rate, in a method analogous to previous approaches (Shin et al., 2010), for at least ten embryos per morpholino dose.

In situ hybridization was performed as previously described (Thisse and Thisse, 2008).

For β-catenin and Connexin 43 analysis, isolated hearts from 72 hpf embryos were fixed in Prefer fixative (Anatech) and incubated with primary antibodies mouse β-catenin (1:200; BD Biosciences, 610154) or rabbit Connexin 43 (1:50; Sigma-Aldrich, C6219) followed by goat anti-mouse Alexa Fluor 488 (1:1000; Invitrogen, A11029) or donkey anti-rabbit Alexa Fluor 555 (1:1000; Invitrogen, A31572), respectively. Stained hearts were mounted in ProLong Gold antifade reagent with DAPI (Invitrogen) and imaged using a Leica SP5X laser-scanning confocal microscope at 63× magnification. Images were analyzed using ImageJ.

Hearts were isolated from embryos at 72 hpf and stained with either the transmembrane potential-sensitive dye di-8-ANEPPS (Invitrogen) or the calcium-sensitive ratiometric dye Fura-2AM (Invitrogen) for measurement of voltage or calcium transients, respectively. Resulting fluorescence intensities were recorded with a high-speed charge-coupled device camera (RedShirtImaging) and images analyzed using software implemented in MATLAB (Panáková et al., 2010). Pacing at a constant rate of 80 beats per minute was used for calcium mapping to eliminate heart rate effects. For all comparisons, regions of interest were determined and compared between morpholino-injected embryos and controls from at least two separate mate pairings using a two-sided Student’s t-test.

RNA was isolated from 48 hpf wild-type embryos and used to synthesize cDNA as previously described (Peterson and Freeman, 2009). tmem88a-specific primers containing the MultiSite Gateway attB1 and attB2 recognition sites were then used to amplify tmem88a from cDNA by PCR. To generate the abridged tmem88a transcript (tmem88aΔVWV), the reverse primer carried a stop codon immediately before the C-terminal VWV motif. The verified PCR product was combined with pDONR221 in a BP reaction (BP Clonase II; Invitrogen) to generate MultiSite Gateway entry vectors pENTR_tmem88a and pENTR_tmem88aΔVWV. These entry vectors were combined with p3E-polyA [Tol2kit #302 (Kwan et al., 2007)], a Tol2 destination vector containing an alpha-crystallin:YFP selection marker, and pENTR5′_ubi [containing the ubiquitin promoter (Mosimann et al., 2011)], in an LR reaction (LR Clonase Plus II; Invitrogen) to generate the ubi:tmem88a and ubi:tmem88aΔVWV expression vectors. Zebrafish embryos were injected with 20 pg of either of these vectors and 10 pg transposase at the one-cell developmental stage, and screened for fluorescence and developmental phenotypes at 48 hpf.

Test morpholinos were injected either singly or in combination with a previously verified morpholino targeting wnt11 (Panáková et al., 2010) at low doses to identify potential synergistic effects. Resulting morphants were scored for cyclopia at 48 hpf based on the five classes previously described (Marlow et al., 1998). The evaluator (G.M.) was blinded as to the identity of the morphants at the time of survey. Percentages reflect the total embryos in each cyclopia class summed over three separate mate pairings. Results were compared using a Chi-squared test.

We thank all members of the C.A.M. and F.P.R. labs for helpful discussion, especially Dr Daniela Panakova and Dr Andreas Werdich, as well as Jason Montojo from the Q.M. lab. We also thank Dr Rahul Deo (UCSF), Dr Martha Marvin (Williams College), Dr Andrew Emili (University of Toronto) and Dr Natalie Morson (University of Toronto) for helpful discussions and critical reading of the manuscript.