Diverse yet highly selective interorgan crosstalk mechanisms shape the bodywide transcriptome landscape

The Thomas N. Sato BioMEC-X Laboratories, Advanced Telecommunications Research Institute International (ATR), Kyoto, Japan; ERATO Sato Live Bio-Forecasting Project, Japan Science and Technology Agency (JST), Kyoto, Japan; Kyoto University, Graduate School of Biostudies, Kyoto, Japan; Division of Cardiovascular Medicine, Department of Medicine, University of California, San Diego, La Jolla, CA, USA; Department of Biomedical Engineering, Cornell University, Ithaca, NY, USA; Centenary Institute, Sydney, Australia *These two authors contributed equally to this work


INTRODUCTION
The organ-to-organ crosstalk is critical for development and physiological homeostasis (Droujinine and Perrimon, 2013). Failure or dysfunction of such interorgan communication system severely hampers the development and functions of virtually all organs (Droujinine and Perrimon, 2013;Noble, 2002;Rajan and Perrimon, 2011;Stainier, 2001). Despite such importance, a comprehensive understanding of the interorgan signaling mechanism regulating the diverse organ functions is lacking.
Each organ may discriminate among the various interorgan signals. In addition, the differential combination of interorgan signals may be delivered to each organ via differential mechanisms. The molecular and cellular underpinnings of such bodywide organ-to-organ interactions are important but defining them remains challenging. Furthermore, only a small subset of interorgan interactions are currently known and the bodywide landscape of such interactions are lacking. This is, in part, due to the fact that such problems have been traditionally studied based on the previous knowledge about the interorgan crosstalks. In contrast, a bodywide and unbiased approach may provide an advantage and that comprehensive bodywide profiles representing the molecular-level changes in each organ resulting from the specific ablation and/or mutations of the interorgan crosstalk signals would be extremely useful. However, this approach is intractable, in particular for vertebrate models with closed circulatory systems, as organs do not form in the absence of the cardiovascular system, a primary mediator of the interorgan communication.
One exception is zebrafish (Stainier, 2001;Weinstein and Fishman, 1996). While zebrafish possess a closed circulatory system, organ-formation appears to initiate and proceed at least to a certain extent even in the absence of the cardiovascular system (Liao et al., 1997;Stainier, 2001;Weinstein and Fishman, 1996). Hence, zebrafish provide a model vertebrate system for 4 such an analysis. Zebrafish also provide a relatively simple vertebrate model system where genetic tools are already available for organ-specific ablations and systematic mutations of organ-specific genes (Stainier, 2001;Weinstein, 2002).
Here, we used this vertebrate model organism to generate and characterize the dynamic changes of the comprehensive bodywide transcriptome landscapes upon the ablations of the heart and the vasculature, two primary components of the cardiovascular system, and also the liver, the largest organ regulating bodywide metabolism. We also characterized the bodywide transcriptome landscapes of the mutants for 3 cardiac-specific genes, each of which exhibit differential contractile phenotype. Additionally, the effects of the mutations of 9 liver-and liver/kidney specific genes were individually characterized. These datasets were integrated into a mechanistic model describing how differential interorgan crosstalk mechanisms could regulate the distinct functions of each organ and also the distinct processes within the same organ. The experimental validation of the model unveiled unexpectedly diverse yet selective interorgan crosstalk mechanisms that differentially regulate 73 genes implicated to support and/or guide the differential organ development and function across 8 organs. The findings are summarized as a bodywide landscape depicting the interorgan genetic interaction network that serves as a platform resource for studying organ-to-organ interactions at the organismal level. This landscape is also made available at i-organ.atr.jp as a public resource.
The specific ablation of cardiac muscle, but not endocardium/endothelial cells, was confirmed by the lack of cardiac muscle fluorescent reporter expression (Fig. 1a), and the presence of 5 endothelial reporter gene expression, fli1a:efgp (Isogai et al., 2003;Lawson and Weinstein, 2002) (Fig. 1a). The "heartless" larvae at 4.5 dpf (days post fertilization) were processed for genome-wide transcriptome analysis.
The wholebody RNAseq analyses followed by the Gene Ontology (GO) enrichment analysis found differentially regulated genes including those related to stimulus via external cues, stress, immune responses, and homeostasis ( Fig. 1b, raw data available at i-organs.atr.jp).
The RNAseq data, followed by the qRT-PCR validation, confirmed the significant reductions in known cardiomyocyte-specific genes (myl7, cmlc1, tnnc1a, nppa, etc.) (Fig. 1c, Table S1), further validating the cardiomyocyte-specific ablation. The analysis identified 54 genes whose expression levels were reproducibly altered in non-cardiac organs in "heartless" larva ( Fig. 1c, Table S1). Wholemount in situ hybridization (WISH) analyses identified their expression patterns: 22 specifically in the liver, 3 in the liver and the kidney, 1 in the somite, 1 in the pancreas, 1 in the liver and the somite, 2 in the liver and undefined tissues, 1 in the somite and undefined tissues, 23 in unidentified or undetermined tissues (Figures S1&S2, Table S2), each representing mostly organ-specific gene implicated for the organ-specific physiological functions such as fatty-acid metabolism for the liver, transporters and enzymes for the kidney, a hormone for the pancreas, etc.. Many of the 23 genes found to be upregulated in unidentified or undetermined tissues encode proteins involved in immune responses, inflammation and/or tissue remodeling. The liver-specific expression was also confirmed by transcriptome analysis of the liver-specific ablation line (Curado et al., 2007;Curado et al., 2008) and also of the isolated liver tissue (see METHODS, Figure S3, Table S1).
Next, we determined whether a single cardiomyocyte-specific gene mutation is sufficient to induce the same effects as in "heartless". Three genes encoding contractile cardiac-muscle-specific proteins were individually mutated (Fig. 1d). Mutations were introduced into myl7, cmlc1 and tnnc1a genes by the CRISPR/Cas9 system (Gagnon et al., 2014;Jao et al., 2013;Thomas et al., 2014) (see METHODS, Figure S4a, b). Both myl7 and 6 cmlc1 mutations resulted in the indistinguishable contractile dysfunction of the heart (Movies S1 -S3). The tnnc1a mutation resulted in the lack of ventricular, but not atrial, contraction (Movie S4), correlating with its ventricular-specific expression. No circulation was observed in any of these three mutants (Movies S5 -S7).
The qRT-PCR analysis of the 54 non-cardiac genes differentially regulated in "heartless", 4 cardiac genes (myl7, cmlc1, tnnc1a and nppa), and 6 liver genes that are unaffected in "heartless" was performed in each mutant (Fig. 1d, Table S1). The results with the cardiac gene mutants unveil the common and distinct cardiac properties differentially regulating the gene expressions in each non-cardiac organ ( Fig. 1e): (I) Positive regulation by any one of the three cardiac-genes; (II) Positive regulation by the myl7 mutant, but not by cmlc1 or tnnc1a; (III) Positive regulation by myl7 and tnnc1a, but not by cmlc1; (IV) Negative regulation by all three cardiac genes; (V) Negative regulation by cmlc1 and tnnc1a, but not by myl7; (VI) Negative regulation by tnnc1a, but not by myl7 or the cmlc1.
The differential effects of the mutations on another cardiomyocyte-specific gene, nppa, suggest the genetic interactions among the cardiac-genes (Fig. 1e). In the myl7 mutant, the cmlc1 and tnnc1a expressions were partially reduced. In contrast, in both cmlc1 and tnnc1a mutants, the myl7 expression is not affected, despite the nppa upregulation. These results suggest that the nppa expression is under the negative regulations of both cmlc1 and tnnc1a.
Moreover, both cmlc1 and tnnc1a appear to be under the regulation of yet another identified input pathway(s), as nppa expression is unaffected by the myl7 mutation which partially reduces the cmlc1 and tnnc1a expression (Fig. 1e).
The gene expression effects induced the myl7 mutation were completely rescued (R 2 =0.95501) by the re-expression of myl7 using the myl7 cardiomyocyte-specific promoter ( Figure S4c, d; Table S1), which was accompanied by the complete rescue of contractile and circulatory dysfunctions (Movie S5). This result supports the idea that the effects of myl7 mutation are due to its functions in cardiomyocytes. The re-introduction of cmlc1 using the myl7 promoter ( Figure S4c, Table S1) rescued the cardiac contractility, but not the circulation (Movie S6). The gene expression effects were partially rescued (R 2 =0.76025) ( Figure S4d, Table S1), suggesting that cardiac contraction alone is mostly sufficient to rescue the normal gene expression pattern in the cmlc1 mutant. In contrast, the re-expression of tnnc1a under the control of the myl7 promoter failed to rescue the cardiac functions or the gene expression pattern (Movie S7, Figure S4d, Table S1). This is presumably due to the fact that tnnc1a is expressed in ventricle, but not in atrium, and the transgene tnnc1a under the control of the myl7 promoter is expressed in both ventricle and atrium, resulting in aberrant contractions (Movie S7). The results suggest that the aberrant contraction is not sufficient for the normal gene expression pattern -"the normal" contraction is required.

Liver gene mutants
Many genes influenced in "heartless" are liver-genes (Fig. 1e). Hence, a possibility of their cross-regulation was examined by analyzing the mutants for 11 liver genes ( Figure S5a, b). The expression of the 54 genes regulated in "heartless" was mostly unaffected in any of these mutants (Fig. 2a, Table S1). Only subtle, but statistically significant, influences were detected with a few genes in two mutants: sepp1b and lepa in cbln13 mutant, cbln9 in rpb2b mutant, indicating the presence of only a weak cross-regulations among them (Fig 2a, Table S1). There are other genes that are weakly influenced such as scpp8 in cbln13 mutant, and itln3 in fabp10a mutant (Fig. 2a). The expression of these genes was, however, upregulated in "heartless," while in each mutant their expression is downregulated. Such influences in the opposite directions in "heartless" and the mutant are inconsistent with the possibility of a cross-regulation among these genes.
Since only weak cross-regulations among the liver genes were found at 4.5 dpf, the gene expression was studied in these mutants at later stages where the interorgan crosstalk is presumably more complex. Comprehensive transcriptome analyses were performed using 30 8 dpf juveniles of the liver-gene mutants ( Figure S5a, b). One such candidate gene, bscl2l, was found whose expression was dependent on dpys, a gene expressed in both the liver and the kidney. In the dpys mutant, the bscl2l expression remains unchanged at 4.5 dpf or 10 dpf, but is significantly reduced at 30 dpf (Fig. 2b, Table S1). To identify the organs where bscl2l is expressed and where its expression is affected by the dpys mutation, several organs of 30 dpf juvenile fish were isolated ( Figure S5c, Table S1) and the level of bscl2l in each organ was measured by qRT-PCR. Expression of bscl2l was detected in several organs with the highest level in the pancreas ( Figure S5d, Table S1). The largest reduction was, however, found in the intestine (Fig. 2c, Table S1). The results using the liver gene mutations are summarized in Fig.   2d, illustrating the weak cross-regulations among liver genes and stronger interaction between liver/kidney (dpys) and non-liver gene (bscl2l) (Fig. 2d).

"Vesselless"
The bodywide transcriptome landscape in the absence of the vascular network was studied using cloche (Liao et al., 1997;Reischauer et al., 2016) (referred to as "vesselless" here) (Fig.   3a). The absence of the vessels and endocardium, but the presence of cardiac muscle and contraction was confirmed (Fig. 3a, Movie S8). The comprehensive transcriptome analysis of "vesselless" larva at 5 dpf revealed many vascular and non-vascular genes whose expression was differentially affected. GO enrichment analysis using the RNAseq datasets (available at i-organs.atr.jp) found both commonly and distinctly enriched GO terms between "vesselless" and "heartless" (Fig. 3b). Many of the GO terms enriched for "heartless" are also enriched for "vesselless (Fig. 3b). In addition, there are a number of GO terms that are specifically-enriched for "vesselless" (Fig. 3b.). The most prominent of these "vesselless"-specific terms are those related to sensory system and cell-cell signaling (Fig. 3b). Validation studies by qRT-PCR identified 13 non-vascular genes whose expression was significantly and reproducibly altered in "vesselless", but not or only weakly altered in "heartless" (Fig. 3c, Table S1). Four (ompa, pglyrp2, crx, rx1) were downregulated and nine (hmgcra, sqlea, lss, nsdhl, ebp, cyp51, memo1, scd, fads2) were upregulated (Fig. 3c). All of the upregulated genes encode enzymes critically involved in cholesterol biosynthesis (Lu et al., 2015;Mazein et al., 2013;Paton and Ntambi, 2009) ( Figure S6). WISH analyses indicated that all were expressed in the liver, and five (hmgcra, sqlea, lss, msmo1, fads2) were also expressed in the brain and two (hmgcra, fads2) in the intestine (Fig. 3d). Two downregulated genes, crx and rx1, are eye-genes, confirming the previous report (Dhakal et al., 2015). One newly identified downregulated gene, ompa, is expressed in the olfactory bulb (Fig. 3e). The comparison to WISH patterns of a pan-olfactory-bulb marker, ompb (Celik et al., 2002), shows that both ompa and ompb are expressed in the olfactory bulb, but the ompa expression is restricted to a subset of neuroepithelial cells of the bulb (Fig. 3e). Only weak signals were detected for another olfactory bulb-specific gene, ora1 (Ahuja and Korsching, 2014;Behrens et al., 2014), in wild type and in "vesselless" (Fig. 3e, Table S1). Another newly identified downregulated gene, pglyrp2, a cellular protein known for its important role in innate immune response, was expressed in the liver and multiple other undefined organs/tissues ( Figure S2c). That such differential expressions of the genes were due to the lack of the vessels is further supported by the characterization of etv2 morphant ( Figure S7). Etv2 morpholino injection was previously shown to reduce the vascular network in a relatively specific manner (Craig et al., 2015;Sumanas and Lin, 2006;Veldman and Lin, 2012) ( Figure S6a). Gene expression correlation analysis between "vesselless" and etv2 morphant showed a high correlation coefficient (R 2 =0.81821) (Fig. 3f, Table S1), supporting the idea that the differential gene expressions in "vesselless" is likely due to the lack of vessels, rather than the direct effects of the cloche gene mutation. The results of "vesselless" are summarized in Fig. 3g.

Modeling and experimental validation of diverse interorgan mechanisms
The findings hitherto suggest the cardiomyocytes, distinct cardiac contraction properties and the vessels impose differential influences on the bodywide transcriptome landscape (Figs. 1e & 3g).
To integrate such findings and gain a mechanistic insight, we developed a model describing how the differential gene expression in each organ could be regulated (Fig. 4a) Table S1). Based on these results, each gene is classified to a specific group (A -H) and summarized in Fig. 4d.
All findings in this paper are summarized in Fig. 5, depicting how the expression of each gene in each organ is regulated by distinct and/or common interorgan mechanisms mediated by the cardiovascular system. In the liver and kidney, more genes (13 genes) are under the distantly-acting positive regulations (A, D, E categories) of the cardiovascular system. Such regulation is differentially mediated by cardiac genes (myl7, cmlc1, tnnc1a) or the contraction properties (I, II, III, IV, V) (Figs. 1e&5). In contrast, the genes expressed in less discrete organs/tissues (i.e., "others") are under the distantly-acting negative regulation (B, C, F categories), and they are all negatively-controlled by all these (myl7, cmlc1, tnnc1a) cardiac-genes or contraction (IV). One exception is the mmp13a expression which is mediated specifically by the tnnc1a function (VI). No genes in these tissues/organs appear to be under distantly-acting positive regulation (A, D, E categories). Other organ specific responses can be found in sensory organs. In the olfactory bulb and eye, the cardiovascular system appears to locally-act on subsets of neuroepithelial cells in the positive direction (G category) (Figs. 4&5).
In the olfactory bulb, ompa is under the positive regulation of the vascular endothelial cells. In the eye, crx and rx1 also are regulated by the same mechanism (Fig. 3g). This regulatory mode could not be found in other organs, except for pglyrp2 gene which was expressed in the liver and other tissues (Fig. 5).

DISCUSSION
In the classic view, the heart generates force to deliver systemic factors such as oxygen, nutrients, hormones, etc. to distant organs via vessels to maintain homeostasis of the body (Aaronson et al., 2014;Noble, 2002). In addition, the immune cells reach to the distant organs via vessels upon inflammation (Aaronson et al., 2014;Noble, 2002). M1 (A, B) and M2 (C, D) (Figs. 4a) correspond to these classic cardiovascular mechanisms, respectively. This study unveils an unconventional mechanism, M3 (E, F), where a cardiac-derived mediator acts on distant organ(s) independently of the vessel circulation. What could be the mediator(s)? The possibility of the neural system involvement was examined. We tested this possibility by specifically ablating a subset of dopaminergic (DA) neurons, a neural type that is known to be involved in a number of homeostatic processes including the regulation of cardiovascular functions (Aaronson et al., 2014;Myers and Olson, 2012;Noble, 2002). A subset of DA neurons was genetically ablated by the double mutations of otpa and otpb genes (Fernandes et al., 2013) ( Figure S9a). The expression of the genes that belong to M3 (E, F) was quantitatively analyzed by qRT-PCR ( Figure S9b, Table S1). While the ablation of these DA neurons was confirmed by the downregulation of oxt as previously described (Fernandes et al., 2013), none of the genes in E or F categories exhibit the regulated expressions that are consistent with the possibility of otpa/otpb-positive DA neurons as the mediators (Fig. 9b). The expressions of il4r.1 and timp2b are slightly downregulated in the otpa/otpb double mutant, but the direction of the effects is opposite. Hence, it is unlikely that DA neurons mediate the cardiovascular effects on these genes in distant organs. What are then the mediators? Diffusible cardiomyocyte-derived factors? The identification of such cardiomyocyte-derived mediator(s) may provide an important step towards a more complete understanding of the cardiovascular mechanisms that influence the organ development and functions.
We also identified a vessel-dependent, but independent of cardiomyocytes or the circulation, mechanism (M4: G, H). The expressions of ompa and pglyrp2 were found to be dependent on the presence of vessels (Figs. 3a&5). The ompa expression was detected in a specific subset of the olfactory neural epithelial cells (Fig. 3e). The pglyrp2 expression appears to be in multiple tissues/organs including the liver ( Figure S2c). While the critical roles of the vessels in the development of the liver and the pancreas have previously been shown (Lammert et al., 2001;Matsumoto et al., 2001), the question of whether such local effects are independent of the circulatory force remained controversial. The present study shows that the circulation-independent local-vessel-mediated mechanism, at least in part, contributes to the ompa and pglyrp2 expression. Further, neurovascular interactions are critical in development and disease (Weinstein, 2005). The data show that the local vessel-mediated expressions of ompa & crx/rx1 in the olfactory bulb and the eye, respectively. Hence, it is possible that these neural tissues are more susceptible to the local presence of the vessels.
The expression of the genes encoding enzymes for cholesterol biosynthesis is also under the negative local regulation by the vessels (Figs. 3g & 5). Cholesterol in the circulation is taken up by the vascular endothelial cells via endocytosis (Anderson et al., 2011;Ho et al., 2004).
This endothelial mechanism is critical to maintain the homeostasis of cholesterol level in the circulation (Anderson et al., 2011;Ho et al., 2004). It is possible that such a mechanism operates as a negative feedback, suppressing the expression of the genes encoding enzymes for 13 cholesterol biosynthesis -hence, the absence of the vessels results in the upregulation of these genes. The treatment of the wild type or "heartless" zebrafish larvae by atorvastatin (D'Amico et al., 2007;Mathews et al., 2014) resulted in the upregulation of these genes ( Figure S10, Table   S1). The atorvastatin treatment of "vesselless" resulted in only a weak or no upregulation of expression of these genes ( Figure S10, Table S1), supporting the notion that the expressions of these genes are in a negative feedback loop where the vascular endothelial cells function as a suppressive interface in an circulation-independent manner.
In conclusion, the present study provides a genome-wide and whole-body level profiles unveiling the unexpected divergence and coherence in the interorgan communication mechanism regulating the gene expression in various organs. The findings provide an insight into how differential organ development and function may be regulated at the organismal level by the cardiovascular system, a primary mediator of the interorgan crosstalk. The bodywide interorgan transcriptome landscape shown here also serves as a platform resource for studying organ-to-organ interactions at the organismal level. The expression of the genes shown in Fig. 5 could also serve as organ-specific functional reporters for the functional states of the cardiovascular system for the studies of vertebrate development. Many diseases are caused by the changes of the functional states of the cardiovascular system (Noble, 2002). Hence, it is possible that the altered gene expression shown in this paper could be exploited to evaluate the effects of a therapeutic treatment on such diseases.

Fish husbandry
Zebrafish were maintained in circulation-type aquarium system (Iwaki) with 14 h/day and 10 h/night cycle at around 27 o C. The fertilized eggs were collected and raised at 28.5 o C in egg water (0.06% artificial marine salt supplemented with 0.0002% methylene blue) until around epiboly stage and subsequently in 1/3 Ringer's medium (1.67 mM HEPES, 38.7 mM NaCl, 0.97 14 mM KCl, 0.6 mM CaCl 2 , pH 7.2) containing 0.001% phenylthiourea (PTU) (Sigma) to prevent pigmentation. Embryos and larvae were staged to days post fertilization (dpf) according to Kimmel et al (Kimmel et al., 1995). Zebrafish maintenance and experiments were conducted in accordance with animal protocols approved by the Animal Care and Use Committee of Advanced Telecommunications Research Institute International (A1403, A1503, A1603).
For the liver-specific ablation, Tg(fabp10a:CFP-NTR) homozygous fish were used. Tg (fabp10a:CFP-NTR) embryos/larvae were treated with 7 mM MTZ prepared as above from 2.5 dpf to 5.5 dpf. MTZ-containing media was changed everyday during the treatment.

Organ dissection
A liver was dissected out from a larva at 4.5 dpf under stereotype fluorescent microscope. To isolate a liver free of other organs/tissues, we used Tg(lfabf:DsRed;elaA:egfp) that labels hepatocytes with DsRed fluorescence. The dissected liver and remaining body were placed in a 1.5 ml tube and quickly frozen on dry ice.
After fasting for one day juvenile zebrafish at 30 dpf were anesthetized and the heart, the liver, the intestine, the kidney and the pancreas were dissected and isolated under the stereomicroscope. The organs and remaining body were placed in a 1.5 ml tube and then quickly frozen on liquid nitrogen.

RNA extraction
To obtain total RNA, embryos and larvae were harvested in a 1.5 -2 ml tube at appropriate stage and frozen in liquid N 2 to be stored in -80 o C. To prepare total RNA for RNAseq analysis of embryonic and larval stages, 10 -20 embryos and larvae were pooled in a 1.5 ml tube, and total RNA was isolated using RNeasy Mini Kit (QIAGEN). The pooled embryos/larvae were homogenized in Buffer RLT included in the kit using 5 ml syringe and 24G needle by passing through the needle for 20 times. After homogenization, the manufacture instruction was followed. To prepare RNA for real-time PCR analysis, embryos/larvae were individually harvested in a 1.5 ml or 2 ml tube and total RNA from each individual embryo/larva was isolated by AllPrep DNA/RNA Mini kit (QIAGEN). Individual Embryo/larva was either homogenized using a syringe as described above or crashed in 700 µl Buffer RLT using approximately 50 zirconia balls with 1.5 mm diameter (YTZ balls) (NIKKATO) by centrifuging at 4260 rpm for 60 sec in Cell Destroyer PS1000 (BMS). After the homogenization or crushing, the manufacture instructions were followed. Subsequently, the isolated genomic DNA and total RNA were subjected to genotyping and reverse transcription reaction.

RNA sequencing
Total RNA was prepared from two biological replicate pools of 4.5 dpf wild type and cardiomyocyte-ablated larvae and 5.5 dpf wild type, cloche mutant and the liver-ablated larvae.
Each pool had about 15 larvae. Total RNA of 10 dpf and 30 dpf zebrafish was prepared by mixing equal amount of total RNA isolated from 5 -10 single zebrafish prior to the RNA quality check by Agilent 2100 Bioanalyzer (Agilent, CA, USA). The RNAseq analysis (read length:100 bp, total reads number per sample: about 100 million, single-end read) for Fig. 1 and  (Afgan et al., 2016). The obtained bam file was used to calculate fragments-per-kilobase-of-exon-per-million (FPKM) of transcripts and the differential gene expression data using Cuffdiff (Trapnell et al., 2010). To perform Gene ontology enrichment analysis, the enriched genes were defined as an absolute of log 2 fold-changes is greater than or equal to 1 and p-value is less than 0.05. A gene ontology enrichment analysis is performed by R package "topGO" using a root category "BP" and a reference database "org.Dr.eg.db". To prepare volcano plot graph from RNAseq data, p-value and fold-change was calculated using DESeq2 (Love et al., 2014) with default settings.

Whole-mount in situ hybridization (WISH)
To synthesize an antisense RNA probe, the template DNA was amplified by PCR using

Two-color fluorescent in situ hybridization
DNP-labeled RNA probe was synthesized using T7 RNA polymerase (Roche) by incubating template DNA with 0.35 mM DNP-11-UTP (Perkin Elmer), 1 mM of ATP, GTP and CTP, 0.65 mM UTP (Invitrogen) for 3 hrs at 37 o C. The synthesized DNP-probe was purified by LiCl precipitation.
Embryos were prepared and hybridized with DIG-labeled and DNP-labeled riboprobe as described in WISH using AP-system, except for addition of 5% dextran sulfate (Sigma) to the hybridization buffer. The hybridized embryos were washed and blocked as WISH method of for imaging under confocal microscope.

Microscopy and image process
To observe and record the heart beating and the circulation, embryos/larvae were anesthetized using 0.012% MS-222 and were mounted either laterally or ventrally in 1.0% NuSive GTG Agarose (Lonza) on glass-bottomed 35 mm dish. Imaging was performed using a 10x dry objective lens (Plan Apo, NA0.45) or 20x multi immersion objective lens (Plan Fluor, NA0.75) (Nikon) mounted on Nikon A1R confocal microscope. Time-lapse image was recorded with a resonant scanner for 30 f/s imaging, and converted to QuickTime movie using IMARIS software (BITPLANE). These movies were converted into mp4 movies using iMovie software.
To take a image of whole mount in situ hybridization, specimens were mounted in 75 -80% glycerol and imaged using 4 x (Plan Apo/NA0.20) or 10 x (Plan Apo/NA 0.45) (Nikon) objective lens mounted on Nikon eclipse inverted microscope and 1 x objective lens (Plan Apo) mounted on Leica M165 FC microscope.
To make a large image in Fig.1a, Fig. 3a, Figure S3, the tiled images of an embryo/larvae were taken using a 10x dry objective lens (Plan Apo, NA0.45) mounted on Nikon A1R confocal microscope. The obtained images were assembled using MosaicJ, a plugin of Image J. The assembled image was imported into Adobe Photoshop.

CRISPR/CAS9 mutagenesis
For CRISPR/Cas9, sgRNAs were designed using the online tool CHOPCHOP (http://chopchop.cbu.uib.no/#), CRISPR DESIGN (http://crispr.mit.edu/) and CRISPRdirect (https://crispr.dbcls.jp/). Target sequences and guide RNA (gRNA) sequences were listed in Table S2. For preparing gRNA, we followed either plasmid-based method, where the template sequence for gRNA was cloned in plasmid, or cloning-free method. For the plasmid based method (Jao et al., 2013), two complementary 20 µl base oligonucleotides corresponding to the target sequence were annealed in 20 µl solution with 1x NEBuffer3 by the following procedure: min. The resulting double strand DNA was purified using QIAquick PCR purification kit (Qiagen). The gRNAs were transcribed using MEGAshortscript kit (Ambion). gRNAs were treated with DNase, which is included in the kit, and precipitated using lithium-chloride precipitation solution (Ambion). For making Cas9 mRNA to co-inject with gRNAs to zebrafish eggs, we used pCS2-nls-zCas9-nls (Jao et al., 2013) (Addgene) as a template DNA. The template DNA was linearized by NotI (NEB) and purified using QIAquick PCR purification kit.
Capped nls-zCas9-nls mRNA was synthesized using mMESSAGE mMACHINE SP6 transcription kit (Ambion) in a volume of 20 µl. The synthesized mRNA was treated with DNase and precipitated with lithium-chloride precipitation solution (Ambion).
To assay and determine the indel mutation by gRNAs, the genome obtained from individual embryos of 1 to 4 dpf or tail clipping for direct sequencing and/or high resolution melt (HRM) analysis was used (Thomas et al., 2014). Embryo or tail fin clips were transferred into 25 to 50 µl of lysis buffer (10 mM Tris-HCl (pH 8.0), 50 mM KCl, 0.3% Tween20, 0.3% NP40, 1 mM EDTA, 0.2 mg/ml Proteinase K (Invitrogen)), and incubated at 55 o C for 2 hrs to overnight, followed by the incubation at 98 o C for 10 min. For direct sequence, the sequence spanning the gRNA target site was amplified by PCR using 1 µl of the lysed sample as a template, and purified the PCR product. Primers used for the PCR and direct sequencing are listed in Supplementary Table 2 Table   S2.

Atorvastatin treatment
Atorvastatin Calcium Trihydrate (Wako) was dissolved in 100% DMSO at a concentration of 10 mM. Drug was diluted in 0.001% PTU to make the working solution of 2 µM with 0.2% DMSO. Atorvastatin treatment was initiated at 24 hpf and replaced with fresh drug every 24 hrs to 4.5 dpf. To combine atorvastatin treatment and heart ablation, embryos were first treated with atorvastatin to 3 dpf. Then the embryos were soaked in the mixture of 2 µM atorvastatin and 10 mM MTZ for 6 hrs at 3 dpf, followed by the incubation in atorvastatin solution again to 4.5 dpf.

Data analyses and statistics
For data collection and analysis of qRT-PCR, no statistical methods were used to predetermine sample size. Embryos/larvae subjected to qRT-PCR analysis were blindly collected from cluchmates. The etv2 morphants were identified by the reduced expression of Tg(fli1a:egfp) reporter and processed for the analyses. For collecting larvae in Figure S4, those expressing GFP widely in heart (more than 70% in heart in appearance under fluorescence stereo microscope) were collected for qRT-PCR analysis. For qRT-PCR data analysis of mutant fish of myl7, cmlc1, tnnc1a and dpys in Fig. 1  and Figure S9, we considered the Ct of lepb as 46, if the signal was not detected. Ora1 gene also showed a similarly low expression in WT -the signals for one of eight WTs and one of eight cloche mutants with MTZ treatment were undetectable. Thus, for these samples, the Ct was also considered as 46. For drug treatment in Figure S10, embryos/larvae were randomly selected and processed for each treatment group. For the WISH experiments, at least 8 individual fish were processed and examined. Student t-test was performed for statistical analysis, and a p-value less than 0.05 was considered to be statistically significant (* p<0.05, **p<0.01 and ***p<0.001).
The horizontal line represents the mean.

SUPPLEMENTARY INFORMATION
Supplementary Information includes 10 figures, 2 tables and 8 movies and can be found with this article online.    and "heartless" (red  in the expression levels are set at x1.5 or x0.66, respectively. *indicates those of which expression levels change less than x1.5 -fold or more than x0.66, but exhibit statistically significant (i.e., p<0.05) changes.  Figure S1. WISH expression patterns of the genes. The patterns are shown for wild type (WT) and "heartless" larvae at 4.5 dpf for each gene. a. Liver genes. b. Liver/kidney genes. The renal expression of dpys in 4.5 dpf larva is shown by co-staining with a marker for proximal convoluted tubule, slc20a1a. c. Genes expressed in the liver, somite and others. d. Somite gene.

AUTHOR CONTRIBUTIONS
qRT-PCR result (top graph) shows the higher expression of tgm2a at 3.5 dpf. e. Pancreas gene.
d. Genes expressed in tissues/organs that are undefined (Others). Scale bars, 500 µm.  showing the liver genes (red). Their expression was significantly reduced in "liverless", except two genes (igfbp1a and hamp). However, their expression in the liver was confirmed by WISH (Extended Data Fig. 1a). It is possible that they are upregulated in the small residual liver tissues in "liverless". for the control cmlc1 -/-), WT and tnnc1a (n=8), respectively.
The purities of the heart, the kidney, the liver, the intestine, the pancreas are evaluated by measuring the expressions of myl7, cdh17, fabp10a, fabp2 and ctrb1, respectively, by qRT-PCR.
Movie S2. This video shows the lack of cardiac contraction and circulation of myl7 mutant.
Movie S3. This video shows the lack of cardiac contraction and circulation of cmlc1 mutant.
Movie S4. This video shows aberrant cardiac contraction and circulation of tnnc1a mutant.
Movie S5. This video shows normal rescued cardiac contraction and circulation of the myl7 mutant.
Movie S6. This video shows partially rescued cardiac contraction but no circulation of the cmlc1 mutant.
Movie S7. This video shows the lack of cardiac contraction and circulation of the tnnc1a mutant following the re-introduction of tnnc1a via the myl7 promoter.