Effects of extended pharmacological disruption of zebrafish embryonic heart biomechanical environment on cardiac function, morphology, and gene expression

Biomechanical stimuli are known to be important to cardiac development, but the mechanisms are not fully understood. Here, we pharmacologically disrupted the biomechanical environment of wild‐type zebrafish embryonic hearts for an extended duration and investigated the consequent effects on cardiac function, morphological development, and gene expression.


| INTRODUCTION
Congenital heart malformations (CHMs) are the most prevalent form of human birth anomalies. Studies suggested that biomechanical forces play an important role in modulating cardiac morphogenesis, and abnormal biomechanical stimuli may cause CHM. The disruption of normal physiological blood flow dynamics was found to cause cardiac defects in chick and zebrafish embryonic models. For example, the ligation of the vitelline vein in the chick embryo led to ventricular septal defects, valve abnormalities, and pharyngeal aortic arches malformations, 1 while ligation of the left atrium in the chick embryo led to hypoplastic left heart, among other defects, 2 even without genetic abnormalities in the subjects. Absence of mechanical loading in the chick embryonic heart also resulted in disorganized trabeculation and impaired conduction phenotype differentiation. 3 Furthermore, genetically engineered zebrafish embryonic hearts lacking early blood flow did not develop cardiac cushions. 4 Injection of a glass microbead into an embryonic heart impaired cardiac looping and ultimately resulted in severely defective heart function. 5 Biological mechanisms for how abnormal blood flow dynamics can lead to malformations were further demonstrated. Removal of cardiac contractility and flow via both tnnt2a morpholino and 2,3-butanedione monoxime (BDM) were observed to down-regulate expression of nrg2a, 6 miRNA21, 7 and miRNA143. 8 Nrg2a was shown to be necessary for the proper formation of cardiac trabeculations, 6 while miRNA21 and miRNA143 were shown to be important for proper valve and ventricular development. 7,8 Several other mechanosensitive factors such as klf2 9 and piezo1 10 were similarly found to be important for proper embryonic heart and valve development. These studies significantly enhanced our understanding of the role of mechanical forces in cardiac development.
In investigating cardiac developmental mechanobiology, pharmacological disruption of the heartbeat and blood flow with agents such as BDM is commonly used. BDM has a rapid and reversible effect of inhibiting cardiac contractions, 11 by inhibiting actinmyosin interaction. 12 It has been widely used as cardioplegic agent and shown to protect myocardium during cardiac surgery 13,14 and is unlikely to be toxic.
However, the global effects of prolonged disruption of the hemodynamic environment on morphology and cardiac function have not been fully evaluated. Such analysis can help to improve our understanding of the role of biomechanical environments in cardiac development. Moreover, our understanding of the mechanobiological mechanisms of development is incomplete, and there could be important expression and mechanisms that are yet discovered. It would thus be useful to conduct whole transcriptomic analysis of the embryonic heart under heartbeat inhibition to identify candidate expression.
Thus, in the current article, we add on to previous literature by studying the consequences of administering BDM to 2 dpf zebrafish embryos for a prolonged 3-day duration. We characterized the cardiac morphology and function caused by the treatment, including using computational fluid dynamics (CFD) simulations to evaluate details of the biomechanical environments during and after the treatment. We also performed whole transcriptomic RNA sequencing of embryonic cardiomyocytes, qPCR, and histological validation to identify potential pathways for disruption of normal biomechanical environments on the heart development.
2 | RESULTS 2.1 | Embryonic phenotype and cardiac morphology of BDM-treated embryos BDM treatment of embryos from 2 to 5 dpf resulted in pericardial effusion and upward-curved tails ( Figure 1A). Embryos generally displayed shorter body lengths ( Figure 1F), and yolk sizes were approximately 40% larger ( Figure 1G). This suggested that the disruption of cardiac function resulted in reduced nutrient supply from the yolk for embryonic growth, leading to growth retardation. Hearts in treated embryos also exhibited different morphology, having a noticeably larger looping angle ( Figure 1H), indicating impaired cardiac looping. Ventricular sizes were about 30% smaller than that of the controls ( Figure 1I), suggesting slower growth. As expected, prolonged treatment with blebbistatin, another drug commonly used to induce cardiac arrest by reducing actin-myosin affinity, 15 also induced aberrant features similar to those in BDM-treated hearts ( Figure 2). Despite these impairments, the survival rate of treated embryos was high at 5 dpf (88%, 44/50) before BDM was washed out. The survival rate dropped slightly to 82% (41/50) 2-3 hours after washout ( Figure 1E).

| Cardiac function of BDM-treated embryos
Cardiac function was assessed during the BDM treatment to assess direct effects of the treatment, and also after BDM washout at 5 dpf, to determine if the dysfunction persisted. At 2 dpf, when BDM was added to the embryos, the heart rate decreased by 25% within 30 minutes of treatment and stayed at this lower rate until 5 dpf for embryos that maintained a heartbeat ( Figure 3B). By 3 dpf, contractility was drastically reduced from the control cohorts in the minority of hearts that continued to beat, and continued to decline until 5 dpf ( Figure 3C,D). Thus, an extended duration of biomechanical environment disruption persisted in this group of samples. At 5 dpf, before BDM washout, most of the embryos (64%) exhibited no heartbeat ( Figure 3B). The remaining embryos retained very weak contractility ( Figure 3C) with a drastically reduced heart rate ( Figure 3B). Subsequently, 2 hours after the removal of the BDM, a number of embryos regained heartbeat (15% of all embryos, Figure 3A), adding to the 36% that retained their heartbeat during BDM treatment, and there was a modest gain F I G U R E 1 (A) All BDMtreated embryos featured pericardial effusion (red arrowheads) and some manifested an upward-curved tail (yellow arrowhead). (B) A 5 dpf treated embryo was positioned below a wild-type (control), showing an array of phenotypic distinctions, including a shorter body length, yolk retention, and pericardial effusion. (C-C 0 ) The crosssectional area (red dashed line) of the control and BDM-treated ventricles was estimated. (D-D 0 ) The measurement of cardiac looping angle from a control and a BDM-treated heart. A, atrium; V, ventricle. (E) The number of treated embryos which died or survived from 2 to 5 pdf, evaluated for control and for treated embryos at pre-and post-washout. (F) Growth retardation was observed in treated embryos at 5 dpf. (G) The size of the yolk sac among treated embryos was about 40% larger than that of control. (H) Comparison of cardiac looping angle between control and treated group. (I) Relative end-diastolic ventricular size was obtained by normalizing the cross-sectional area to that of control. For -F-I, mean and SDs were presented and two-sample t-test was used, where *P < 0.05 in contractility ( Figure 3C) and heart rate (about 19% gain, Figure 3B), but this was still much lower than controls. Thus, extended BDM treatment led the heart to have reduced heart rate and cardiac contractile function, which could be related to the abnormal morphology observed.

| Cardiac flow dynamics
The mean ventricular outlet velocity of 5 dpf embryonic hearts was measured via high-speed microscopy ( Figure 4A-C). 16 In control embryos, substantial ejected flow and minimum retrograde flow was observed, indicating the development of a fully functional bulboventricular valve. In the majority of BDM-treated hearts, no blood flow was observed. Out of those that retained a heartbeat, greatly diminished ejected flow velocities and substantial retrograde flows were observed, suggesting depressed myocardial contraction and incompetent bulboventricular valve. With the subsequent removal of BDM, approximately half of the treated embryos exhibited heart beats. However, flow velocities were still less than half of that of control embryos, and noticeable retrograde flow continued to occur during late diastole. Prolonged inhibition of the cardiac function by BDM was thus likely to have caused bulboventricular valve maldevelopment and incompetence.
The ejected stroke volume and cardiac output were calculated from velocity measurements via the assumption of a circular outflow channel and parabolic flow profile. Due to observed regurgitation (Figure 4), BDMtreated embryos with heart beats were severely depressed in both measures, and BDM washout only improved them marginally, to about 10% of that in the control group ( Figure 4E,F). The cardiac output was more severely depressed than ejected stroke volume because treated embryos also had lower heart rate ( Figure 4D).
Thus, consequent to prolonged BDM treatment, the embryonic heart became severely depressed in its cardiac flow dynamics and flow function, which were likely due to undersized ventricles, maldeveloped and incompetent valves, low heart rate, and impaired myocardial contractility.
Microscopy-based CFD flow simulations were conducted for four 5 dpf treated hearts, two before BDM washout (hearts I and II), and two after (hearts III and IV), using previously established methods. 17 Simulations of flow was conducted only for the fluid volume in the ventricle, with the atrioventricular canal as the inlet and the outflow tract as the outlet. Simulations were only conducted for hearts with a heartbeat, as those without will be entirely devoid of biomechanical features. Biomechanics evaluation before BDM washout was indicative of biomechanical stimuli exerted on the heart that might contribute toward its final morphology, and evaluation after BDM washout indicated biomechanical functions associated with the final morphology. Flow features were compared with three healthy ventricles from our previous work, 17 which acted as control samples.
3D wall shear stress (WSS) color contours are shown in Video S1, Supporting Information. The WSS vs. time waveforms are shown in Figure 5. WSS contours demonstrated that the WSS spatial pattern appeared similar to control embryos, with elevated WSS occurring at the inlet and outlet of the ventricle.
During BDM treatment (hearts I and II), systolic WSS was significantly depressed. For example, peak forward flow WSS at the outlet region was only around 80% of that found in controls ( Figure 5B vs. Figure 5F). This was likely due to the depressed cardiac contractility, as the F I G U R E 2 (A) Morphology of blebbistatin-treated hearts at 5 dpf. Blebbistatin-treated hearts exhibited a significantly ($27%) smaller end-diastolic ventricular size (A) and an abnormally large cardiac looping angle (B). Mean and SD were presented and twosample t-test was used, where *P < 0.05 severity of the WSS depression was approximately similar to the severity of contractility depression ( Figure 5E). However, during diastole, time averaged WSS (TAWSS) showed that the outlet regions suffered greater depression in WSS than the inlet region ( Figure 5H). However, it should be noted that in these embryos that retained their heartbeat, significant magnitudes of WSS and flow velocities were still observed, despite being reduced from controls.
With the BDM washout, however, WSS recovered, and the peak WSS in the outlet region reached approximately 80 dyne/cm 2 , which was close to peak magnitudes in control hearts ( Figure 5D vs. Figure 5G). TAWSS were also observed to be close to those in control hearts ( Figure 5H). This recovery suggested that the maldeveloped hearts could still generate substantial fluid forces, despite having depressed function.
However, WSS temporal waveforms demonstrated that WSS were oscillatory in all treated embryos, whether before BDM washout or after, having several reversals in the WSS directions over the cardiac cycle ( Figure 5E-G). This was likely due to valvular incompetence. The oscillatory shear index (OSI) is a good way to quantify the extent of WSS oscillatory nature, and has a range of 0 (effectively unidirectional) to 0.5 (fully oscillatory). 18 Valvular incompetence led to elevated OSI values of between 0.07 and 0.43 ( Figure 5I), and there was little difference from before BDM washout to after. In comparison, control cases had little regurgitation ( Figure 4A) and the OSI were close to zero. The exact OSI for normal hearts could not be obtained from the flow simulations, as their boundary conditions enforced the condition of no valvular regurgitation. However, from the velocity waveforms in Figure 4A, we estimated the oscillatory index of the outflow velocities and found them to be low, around 0.02, which effectively approximated the unidirectional flow condition.

| Cardiomyocyte gene expression
The gene expression data of the GIS and Novogene samples were collected in a single data frame and were jointly analyzed. PCA analysis of the top 3,000 variable genes ( Figure 6A) showed that the samples are separated by both the biological conditions and the technical effect of the different sequencing centers. Further PCA analysis after normalization ( Figure 6A) to remove technical effects showed that it was successfully minimized, and the biological conditions could be isolated. The differential expression and gene set enrichment analysis results of the RNA sequencing of the embryonic cardiomyocytes are shown in Tables 1 and 2, and detailed results are available in Supporting Information spreadsheet. Comparing the BDM-treated myocardium to control myocardium, a total of 2,620 differentially expressed genes F I G U R E 3 (A) The number of zebrafish embryos that survived with no heartbeat or mild contractions was assessed. Number of embryos with mild contraction increased from 36% (16/44) to 51% (21/41), suggesting that the washout partially restores heart contractility. (B) The heart rate (n = 7) dropped by 25% within 30 minutes of treatment at 2 dpf. At 5 dpf, the heart rate of treated hearts before washout (n = 15) was far lower than that of control (n = 20) but washout increased heart rate by 19% (n = 20). Note that the treated embryos with no heartbeats were excluded in this calculation. (C) Ventricular fractional area changes from systole to diastole of control (n = 10) and treated hearts (n = 10), calculated from 2D microscopy, evaluated based on the end-systole which was normalized to the end-diastole at 2, 3, and 5 dpf, respectively. For B-C, mean and standard deviations were presented and a twosample t-test was used, where *P < 0.05. (D) Fractional shortening, a metric of myocardial contractility, measured from cardiac motion tracking of 3D LSFMM microscopy images, taking the end-diastole state as the zero stretch reference, at 3 and 5 dpf (DEGs) were identified, based on the criteria that the FDR ≤0.05 and that the log 2 FC ≥ 1 or ≤À1. Out of these, 922 were significantly up-regulated genes and 1,698 were significantly down-regulated genes, and some of them are listed in Table 1. All DEGs were mapped to KEGG databases to investigate the predicted functions of the DEGs. Table 2 shows the list of significant KEGG pathways with less than 5% false discovery rate. A small selection of genes was selected for RT-qPCR validation and provided positive validation (Figure 7), having a Pearson correlation coefficient of 0.927 with the RNA-seq data ( Figure 8), and thus suggesting high quality and robustness of the RNA-seq. The primer sets used for qPCR can be found in Table 3.
From the results, a strong down-regulation of calcium-dependent signaling pathways was observed. This was expected, given that BDM could alter calcium currents while inhibiting cardiac contractility. 19 A second observation was that treated hearts had up-regulated inflammatory responses, with a significantly up-regulated TLR KEGG pathway, and up-regulation of several inflammatory genes such as tlr5a, tlr5b, tnf-α, and various interleukins. Further, apoptosis responses appeared to be up-regulated as well, with the apoptosis and p53 KEGG F I G U R E 4 Controls (n = 5) and treated hearts before (n = 5) and after (n = 5) washout were quantified for (A-C) ventricular outlet velocity, (D) heart rate, (E) ejected stroke volume, and (F) ejected cardiac output. For (D-F), mean and standard deviations were presented and a two-sample t-test was used, where *P < 0.05. Note that the heart rate plot (D) differs from Figure 3B as the chosen embryos were specifically assessed for other parameters and the sample size was limited to five F I G U R E 5 Legend on next page. pathways, and a number of apoptosis related genes such as fas, casp7, casp8, and casp9 being up-regulated. On the other hand, genes known to be linked to cardiomyocyte proliferation and differentiation were down-regulated, including gli1, sox6, and nrg1. There was also evidence for down-regulated ECM remodeling, with cell-adhesion molecules and ECM-receptor interaction KEGG pathways being down-regulated, and a substantial number of F I G U R E 5 (A) LSFMM image and spatial WSS contour of 5 dpf treated heart I (pre-washout) during systole and diastole, respectively. A, atrium; V, ventricle. Spatial-temporal WSS assessment of controls at 5 dpf (hearts I, II, and III), and treated hearts at 5 dpf before and after BDM washout. Treated hearts I and II were collected before BDM washout and treated hearts III and IV were analyzed after BDM removal. Area-averaged WSS experienced at the inlet, mid-ventricular, and outlet were measured over a course of cardiac cycle for controls and treated hearts (B-G), respectively. ECM fiber production genes being down-regulated, including many types of collagen, tenascin, laminin, and thrombospondin. On the other hand, ECM degradation genes such as mmp9 and mmp13a were up-regulated. There also appeared to be a down-regulation of contractile apparatus synthesis, with the cardiac muscle contraction KEGG being down-regulated, and with genes such as myhz1.3, myh11b, tnnc1b, tpm1, and tpm2 downregulated. A number of genes related to mechanosensitive factors were also down-regulated, including itgb1b.2, klf2a, and rho. Finally, a number of genes related to congenital heart malformations were also observed to be down-regulated. These included hey1 and heyl, the combined loss of which were known to cause impaired epithelial-to-mesenchymal transition leading to congenital heart defects, 20 and gli3, which is known to be related to ventricular septal defect. 21 Previous studies reported that reduced biomechanical stimuli in the zebrafish embryonic heart can prevent trabeculation formation through inhibiting the notch1 and Erbb2 pathways. 22,23 Consistently, we found that treated hearts did not develop obvious trabeculations ( Figure 9A-C), both with BDM or blebbistatin treatment. Interestingly however, RNA-seq results showed that Erbb2 was not significantly differentially expressed, but Erbb4a was significantly down-regulated.
A TUNEL assay ( Figure 10) generally corroborated the RNA analysis results, as the treated embryonic hearts demonstrated elevated apoptosis both in the myocardium and the endothelium, while there was little apoptosis found in control hearts. Apoptosis signals were found in both BDM-and blebbistatin-treated hearts, even though both drugs acted by different mechanisms, suggesting that the lack of mechanical stimuli was likely the common cause of apoptosis. Further, immunostaining for Tnf-α ( Figure 11) demonstrated elevated expression of this inflammatory markers in both the endothelium and the myocardium, corroborating the RNA sequencing results that indicated elevated inflammatory responses.

| DISCUSSION
In the current article, we investigated the effects of extended disruption of the cardiac biomechanical T A B L E 1 Selected cardiomyocyte genes which were significantly up-or down-regulated, based on the criteria that significance change should have FDR ≤0.05 and that the log 2 FC ≥ 1 or ≤À1 (data can be found in Supporting Information) F I G U R E 7 Real-time RT-qPCR validation of a subset of differentially expressed genes between control and BDM-treated hearts.
Expression of genes of interest were first normalized to a housekeeping gene (actb) and then plotted as relative to expression levels in control. Data were presented as mean ± SE. Sample size was three for (A-I) but limited to two only for (J, K) as the expression level in some treated samples for col15a1b and erbb4a were too low to be detected. For (A-I), statistical test was one-way ANOVA test, where *P < 0.05 BDM has a strong depressing effect on zebrafish embryonic heart function and was found to substantially reduce the contractile myocardial deformations and flow stresses in the heart, making it a suitable intervention to test the role of biomechanical stimuli on cardiac development. At the end of the 3-day treatment, we found that hearts were maldeveloped, having smaller sizes, impaired looping and valve incompetence. After the washout of BDM to remove its inhibitory effects, we found that the maldevelopment resulted in a significantly reduced heart rate, contractile function, and due to valve incompetency, severely reduced cardiac output and abnormally oscillatory flow and wall shear stresses.
We believe that these impairments are in part due to the lack of biomechanical stimuli, as there is much evidence on the role of mechanics in enabling proper cardiac development, as listed in our Introduction section. It is possible that the maldevelopments were due to extended pharmacological disruption of myocyte cell function, such as disrupted calcium cycling 24 and crossbridge dynamics, 25 rather than having a biomechanical origin. However, evidence of reversibility of BDM's effects 11 reduces this possibility. It is further possible that these cardiac maldevelopments were the manifestation of growth retardation due to reduced nutrient supply from a disrupted cardiac function. However, the heart is situated next to the yolk, and is likely to still obtain sufficient nutrients, and due to the small size of the embryo, oxygen diffusion limits are unlikely to pose significant restriction on development.
It is important to note that under the influence of BDM, reduced cardiac function did not mean an abolishment of WSS. Not all treated embryos had a complete stoppage of the heart's motion, and in embryos that retained mild contractions and had little forward flow, substantial WSS magnitudes were observed. This can explain why in earlier investigations, seemingly abolished blood flow by BDM can still enable endocardial ring formation, suggesting that WSS is still an important factor to cardiac development. 4 Transcriptomic analysis provided a deeper insight into changes in expression profile of the cardiomyocytes due to an altered hemodynamic environment and provided clues to possible mechanobiological pathways for maldevelopment of the heart due to abnormal biomechanical stimuli. We summarized several possible such pathways in Figure 12 and explain it below. We propose that the biomechanical stimuli primarily came from two sources, stresses from intracardiac blood flow and tissue deformations from myocardial contraction.
From our investigations, BDM either abolished WSS or induced low magnitude and oscillatory WSS on the embryonic endothelium. This pattern of disturbed WSS is known to cause inflammation and atherosclerosis in adult endothelia, 26 and in fetal endothelial cells from umbilical veins. 27 We thus propose that the inflammatory indicators found in the myocardium from our RNA sequencing analysis could reflect the myocardial response to endothelial inflammation, which was, in turn, due to prolonged abnormal WSS stimuli. Since up-regulation of inflammatory markers such as tlr and tnf-α often indicates activation of apoptotic pathways, 28,29 we propose that both endothelial and myocardial inflammation led to the observed elevation in TUNEL signals, and in the up-regulation of myocardial apoptosis genes. Elevated apoptosis could then be playing a role in the underdevelopment of the heart. Interestingly, recent studies showed that inflammation markers including tnf-α, il-6, and crp, 30 and inflammation associated with the activation of nf-κβ were found in children with congenital heart disease, 31 further suggesting inflammation as playing an important role in disease pathophysiology.
We further propose that cardiac inflammation led to down-regulation of ECM production and elevation of ECM degradation. We propose this because MMP9 was reported to be robustly stimulated by increased expression of pro-inflammatory cytokines, 32 and collagen synthesis could be inhibited with increasing activity of proinflammatory cytokines such as Il-1β. 33 At present, however, there is little data on the remodeling responses of the fetal heart to inflammation, and much future work is warranted. In adult hearts, ECM remodeling is also often closely related to inflammation, such as during heart failure or infarction in adults. 34,35 In the adult heart, inflammation typically leads to fibrosis rather than reduced ECM, but the immaturity of the embryonic immune F I G U R E 8 Correlation of log 2 FC of nine genes ( Figure 7A-I) between RNA-seq and quantitative real-time PCR (RT-qPCR). A strong, statistically significant Pearson correlation (R = 0.9269) was shown between the expression levels detected using RNA-seq and real-time PCR. Each data point represented the mean value of three biological replicates system could play a role in the difference observed. Interestingly, remodeling during diseases are known to involve reactivation of fetal genes. 36 Reduced myocardial deformation/strain, on the other hand, are likely to reduce cardiac growth. Enhanced strain in diseased adult hearts led to enhanced growth and remodeling, 37 while in the chick embryo, disrupted left ventricular deformation/strain due to left atrial ligation led to a hypoplastic left ventricular phenotype. 2 We propose that these effects occurred through a reduction of ECM-receptor interaction as suggested by our RNA sequencing results, and through the down-regulation of itg-β, which is the beta subunit of a well-known mechanosensitive receptor 38 with essential intracellular signaling function, 39 and which is crucial to the form and function of the developing myocardium. 40,41 It is well known that integrins and focal adhesion complexes interact intricately with ECM fibers and with intracellular filament proteins such as actin and actinin. A downregulation in itg-β can explain the down-regulation in these intra-and extracellular fibers in our samples. We further propose that the lack of stress stimuli on cytoskeletal stress fibers led to a down-regulation in contractile genes.
Our sequencing results showed that nrg1 and erbb4a were down-regulated, but notch1 and erbb2 were not differentially expressed in the myocardium of treated embryos. A previous article showed that the absence of flow-induced forces activated notch1 and nrg1 pathways in the endothelium to lead to a loss of trabeculations, but similarly did not find a decrease in erbb2, 22 even though erbb2 inhibition was found to inhibit trabeculation formation as well. 23 Studies in mouse embryos suggested that activation of Erbb2/Erbb4 expression by Nrg1 in the myocardium is important to induce trabeculation muscle differentiation and that Nrg1 could be a soluble signal between the endocardium and the myocardium. 42 We thus propose that the lack of flow biomechanical stimuli had down-regulated the notch1 and nrg1 pathways in the endocardium, which signaled a downregulation of nrg1 and Erbb4a in the myocardium, leading to the loss of the trabeculation.
In our current work, we identified potential genetic pathways for abnormal biomechanical conditions to influence morphological development of the embryonic heart. It would be interesting to investigate whether similar pathways are found in congenital heart diseases with disturbed biomechanical conditions, such as hypoplastic left or right heart syndromes. If the pathways are found to be clinically relevant, they could be good molecular targets for pharmacological therapy.

| LIMITATIONS
Several limitations are noted here. Most importantly, we have used BDM to induce cardiac function depression as a means to suppress biomechanical stimuli to the embryonic heart. However, BDM could have pharmacological effects on the cardiac biology, and the possibility that the observed morphology and expression changes were due to pharmacological instead of biomechanical effects cannot be ruled out. Further, these morphological and biological responses could also be affected by a reduction in nutrients supply due to the lack of heartbeat, and this effect could not yet be distinguished from the effects of biomechanical stimuli suppression. For example, BDM has a significant impact on calcium signaling, and this could have contributed to some maldevelopment. In addition, our transcriptomic analysis had focused on the myocardium and did not include the endocardium, and some of the biological pathways proposed in Discussion section remain to be tested in future work. Moreover, our work did not yet cover validation of the proposed pathway in Figure 12, and such future work is warranted. We also did not analyze cell proliferation responses, but this may be relevant. Further, we did not conduct transcriptomic analysis at 2 dpf, but such an analysis could be useful when compared to 5 dpf expressions to understand expressions before and after BDM treatment. It could also be fruitful to conduct transcriptomic analysis and cardiac morphology and function analysis at a later time point after BMD washout, after longer exposure to normalized flow stimuli.

| CONCLUSION
We studied the effects of suppressing cardiac contractility for an extended 72 hour duration via BDM and investigated the resulting morphology, cardiac function, and gene expression. The treatment was found to result in maldeveloped hearts with smaller size, defective looping, incompetent valves, and an absence of trabeculations. These hearts exhibited poor function, low contractility, reduced heart rate, and low cardiac output. Incompetent valves further led to oscillatory F I G U R E 1 1 Immunofluorescence staining of Tnf-α on cryo-sectioned embryonic tissues. GFP-positive myocardium (A-C) and Alexa Fluor 568 dye which indicates Tnf-α (A 0 -C 0 ) are displayed. Relative fluorescence intensity plot of heart region (D) indicates that Tnf-α was significantly up-regulated in treated hearts (n = 5 for each group). Mean and standard deviations are shown, and a two-sample t-test was used, where *P < 0.05. Scale bars: 50 μm F I G U R E 1 2 Schematic diagram showing several possible mechano-responsive pathways of a malformed heart after a 3-day BDM treatment flow and WSS on the embryonic ventricle. Gene expression analysis revealed that myocardium of treated hearts had up-regulated inflammatory and apoptotic responses, and down-regulated ECM-receptor interaction, ECM remodeling, and contractile genes. The absence of trabeculation was found to be accompanied by a down-regulation in myocardium nrg1 and erbb4a, but not notch1 or erbb2. Potential mechanobiology pathways to explain the maldevelopment were proposed.
6 | EXPERIMENTAL PROCEDURES 6.1 | Zebrafish transgenic line A transgenic line, Tg(phiC31.attP.2A, À0.8myl7: EGFP) 43 was bred and raised in compliance with relevant guidelines and regulations approved by the National University of Singapore (NUS) Institute of Animal Care and Use Committee (IACUC; protocol number BR19-0120). This line expresses green fluorescent protein (GFP) in the myocardium, which enables imaging of cardiomyocytes. Fertilized eggs were collected and incubated at 28 C and 0.2 mM 1-phenyl 2-thiourea (PTU) was added to the fish medium to inhibit pigmentation.

| Drug treatment
Ten millimolar of BDM (Sigma, B0753) was used to depress myocardial contractility from 2 to 5 dpf. BDM was washed out at 5 dpf for 2 hours, and both control and treated embryos were observed directly thereafter. For blebbistatin (Sigma, B0560), 5 μM was used under the same experimental conditions.

| Body length, yolk size, and looping angle measurement
The body length, yolk size, cardiac looping angle, and ventricular size of control and treated embryos were measured using ImageJ software. Ventricular size was estimated based on its cross-sectional area.

| Blood flow measurement
Ventricular outlet velocity was estimated according to a previously described method. 44 Brightfield images were taken with a high speed camera, pco.edge sCMOS (PCO AG, Germany) at 1,000-2,000 frames per second. 6.5 | CFD simulation CFD simulation of 5 dpf treated ventricles were performed based on the images taken using line-scan focal modulation microscopy (LSFMM) according to methods established in our previous publication (Method 1). 17 In the current article, zero-pressure reference was set at the inlet throughout the cardiac cycle while measured outlet velocity was set as the boundary condition of the outlet.

| Ventricular contractility quantifications
Ventricular contractility was quantified in two ways, namely fractional area change and fractional shortening. The former was a quantification of systole-to-diastole ventricular area changes from 2D microscopy images using ImageJ software. The latter was a 3D strain quantification from 3D LSFMM images, via a previously established cardiac motion tracking technique. The enddiastolic phase was chosen to be the reference state.

| RNA sample preparation and sequencing
For each RNA sample preparation, myocardial cells were pooled from 50 to 70 larvae from the transgenic zebrafish line, Tg(phiC31.attP.2A, À0.8myl7: EGFP) that were dissociated using the Worthington Papain Dissociation system (Worthington Biochemical) according to the manufacturer's instructions. Embryos were selected regardless of whether a heart rate was present, to ensure that sufficient cells can be harvested as we intended to quantify gene expressions for all treated embryos. GFPpositive cardiomyocytes were flow-sorted using a BD FACSaria II cell sorter (BD Bioscience). About 3,000-15,000 GFP+ cardiomyocytes were sorted directly into lysis buffer solution (Qiagen). RNA was extracted on the same day or cells were stored at À80 C. All RNA isolation experiments were completed within 1 week after FACS-sorting. Total RNA was extracted using RNeasy Micro Kit (Qiagen) and taken through library preparation to be subjected to Illumina sequencing. Three pairs of biological replicates were prepared for sequencing at the Genome Institute of Singapore (GIS) and another two pairs were sequenced by NovogeneALT. To ensure high RNA quality, RNA integrity number (RIN) was measured and confirmed to be greater than 7.0 (details in Table 4).
For sequencing at GIS, cDNA conversion and preamplification was performed using a SMARTer II A oligonucleotide and template switch primer (Clontech, 634833). Only RNA strands with poly-A-tails were converted to cDNA. Library preparation was carried out using Nextera XT DNA Sample Preparation Kit (Illumina, FC-131-1096). "End-repaired" fragments were ligated with dual unique Illumina adapters. All individually indexed samples were subsequently pooled together and multiplexed for sequencing. Libraries were sequenced using an Illumina Hiseq 4000 sequencing system and paired-end 150 bp reads were generated for analysis, with a sequencing output of 50 million reads per sample. For sequencing with NovogeneALT, cDNA synthesis and amplification were performed using SMART-Seq v4 Ultra Low Input RNA Kit (Clontech). Libraries were sequenced with the same specifications as above, but the sequencing output was around 150 million reads for each sample.

| Bioinformatic analysis
For quality control, RNA-Seq pair-end FASTQ reads from both sequencings were assessed using the FastQC's software (https://www.bioinformatics.babraham.ac.uk/ projects/fastqc/). Initially, the metrics "per base sequence quality" and "per sequence quality score" were used to determine the overall base pair and read quality. All reads with quality scores less than Q20 (<1% error rate) or shorter than 30 bps were removed. After inspecting the FastQC's "Adapter Content," adaptor sequences (3 0terminal) were trimmed from the raw sequences by trimmoMatic (version 0.30) using the base quality encoding phred33. 45 The fastq files were PCRdeduplicated by FastUniq. 46 The reads were mapped to the GRCz11 (Genome Reference Consortium Zebrafish Build 11, May 2017) downloaded from Ensembl (http:// m.ensembl.org/Danio_rerio/Info/Annotation), and the transcript quantification was performed with Kallisto v.43. 47 The conversion of transcript to gene counts was performed by the tximport package in R 48 that generated the raw gene counts, the transcript length corrected raw counts and the TPM-normalized counts. The normalized profiles of all available samples were used for principal component analysis (PCA). A PCA was performed based on the top 3,000 variable genes. Figure 1A illustrates a PCA plot for five pairs of biological replicates sequenced by GIS and NovogeneALT highlighting the separation of the biological conditions and need for data adjustment to remove the technical effects of the sequencing centers. The differentially expressed genes were estimated by DESeq2 on the raw gene counts. The model included a sequencing center and a biological condition factor and separated the technical from the biological effects of interest. The differentially expressed genes were detected at a 5% false discovery rate (FDR ≤0.05) and an absolute log fold change (FC) of 1 (jlog 2 FCj ≥1). The TPMnormalized counts were adjusted to remove the effect of the sequencing center from the normalized counts with the removeBatchEffect function of the Limma R package. Figure 1B shows the number of detected genes in each sample from five independent experiments. Hypergeometric test and Fisher exact test were used in statistical tests and the Benjamini-Hochberg procedure was performed for FDR correction.
6.9 | Quantitative real-time polymerase chain reaction Three additional sets of biological replicates were separately prepared for qPCR validation. RNA samples were reverse transcribed to cDNA using RevertAid First Strand Note: RNA concentration of the samples sequenced by NovogeneALT was significantly higher as the FACS condition was made less stringent and the protocol was further improved, including the use of QIAshredder (Qiagen) for better lysate homogenization during the process of RNA isolation and reduced elution volume. These steps resulted in a considerable increase in the harvested cell number and thus, a higher amount of RNA. Each qPCR reaction was prepared by combining 1 μL of cDNA, 10 μL of Power SYBR Green PCR mix (Applied Biosystems, CA), 5 μL of H 2 O, 0.25 μL of 100 μM forward primer (IDT) and 0.25 μL of 100 μM reverse primer (IDT) for a final volume of 20 μL. qPCR reactions were analyzed on a CFX96 Touch Real-Time PCR Detection System (BioRad). Cycling conditions were as follows: 95 C Â 3 minutes, followed by 40 cycles of denaturing (95 C Â 10 seconds) and annealing (60 C Â 30 seconds). Each sample was run in at least triplicate in optically clear 96-well plates or qPCR strips.
The relative expression of each gene of interest was calculated by the delta-delta Ct method on Bio-Rad CFX Maestro. Expression of the gene of interest was first normalized to a housekeeping gene (actb) 49 and was later plotted as relative expression normalized to control. Oneway ANOVA test was used for the statistical comparison of differences between control and treated group.

| Immunofluorescence
Cryostat sections were blocked with 10% sheep serum in PBST for 2 hours at room temperature (RT). Sections were incubated with primary antibodies overnight in a humidified chamber at 4 C (anti-TNF alpha antibody [Abcam, ab1793], 1:400, diluted in 10% sheep serum/ PBST). The next day, sections were washed with PBST (4Â 15 minutes) and incubated with secondary antibodies (goat anti-mouse Alexa Fluor 568 [Abcam, ab175473], 1:500, diluted in 10% sheep serum/PBST) for 1 hour at RT. Stained sections were washed with PBST (4Â 15 minutes) and mounted with Mowiol before imaging under an Olympus FV3000 confocal microscope. Color images were exported using the Olympus FV3000 Fluoview software. Fluorescence intensity of the heart region was measured using MATLAB.