Scale space detector for analyzing spatiotemporal ventricular contractility and nuclear morphogenesis in zebrafish.

In vivo quantitative assessment of structural and functional biomarkers is essential for characterizing the pathophysiology of congenital disorders. In this regard, fixed tissue analysis has offered revolutionary insights into the underlying cellular architecture. However, histological analysis faces major drawbacks with respect to lack of spatiotemporal sampling and tissue artifacts during sample preparation. This study demonstrates the potential of light sheet fluorescence microscopy (LSFM) as a non-invasive, 4D (3days + time) optical sectioning tool for revealing cardiac mechano-transduction in zebrafish. Furthermore, we have described the utility of a scale and size-invariant feature detector, for analyzing individual morphology of fused cardiomyocyte nuclei and characterizing zebrafish ventricular contractility.

Introduction

Zebrafish (Danio rero) are emerging as potent vertebrate model’s for modeling human congenital heart disorders (CHD) (Kula-Alwar et al., 2021, p. 2; Lee et al., 2018; Miura and Yelon, 2011; Vedula et al., 2017a; Yu and Hwang, 2022; Zhao et al., 2020). This is due to numerous attractive traits such as embryonic optical transparency, high fecundity, and ease in genetic or biomechanical modulation for mimicking the human CHD pathophysiology (Choi et al., 2013; Lee et al., 2018; Miura and Yelon, 2011; Rafferty and Quinn, 2018; Tu and Chi, 2012). As a result, zebrafish enable access to phenotypic screening of dynamic biomechanical stimuli such as contractility and blood flow, responsible for modulating heart maturation (Kula-Alwar et al., 2021; Lee et al., 2018; Miura and Yelon, 2011; Tu and Chi, 2012; Vedula et al., 2017a).

Previously conducted zebrafish studies have observed conserved cardiomyocyte count proportional to atrial/ventricular mass or volume per developmental stage (Kula-Alwar et al., 2021, p. 2). Moreover, recent studies suggest cardiomyocyte shape hypertrophy across three distinct ventricular regions—atrioventricular (AV) canal, outer curvature (OC), and inner curvature (IC) regions—apart from distinct atrial cardiomyocyte morphology (Kula-Alwar et al., 2021; Miura and Yelon, 2011; Tu and Chi, 2012). In addition, biologists have questioned the implications of cardiac mechano-transduction on enlarged cardiomyocyte morphology in the OC region in front of AV canal, with respect to spherical (isotropic) cardiomyocytes in the IC (Tu and Chi, 2012; Zhao et al., 2020). However, the ability to observe cardiomyocyte morphogenesis in vivo is adversely affected by tissue birefringence, hindering characterization of beforementioned cardiovascular phenotypes (Bensley et al., 2016; Bray et al., 2010; Ghonim et al., 2017; Teranikar et al., 2020). In this regard, automated feature detectors are proving to be an indispensable tool for segmenting cellular volumes without human intervention to avoid gross inconsistencies and produce refined datasets. (Bolón-Canedo and Remeseiro, 2020; Sargent et al., 2009; Torres and Judson-Torres, 2019).

Conventionally, invasive sectioning procedures have offered revolutionary insights into aberrant tissue up to the cellular scale (Bensley et al., 2016; Javaeed et al., 2021; Teranikar et al., 2022). However, histopathological analysis currently suffers from severe limitations, primarily disruption to tissue homeostasis (Kleczek et al., 2020; Teranikar et al., 2022). With respect to these drawbacks, the optical sectioning modality light sheet fluorescence microscopy (LSFM) has proved instrumental in probing dynamic organogenesis several millimeters inside tissue.(Ding et al., 2017; Fei et al., 2019; Lee et al., 2016; Teranikar et al., 2020; Vedula et al., 2017b). LSFM has tremendously benefitted embryologists to become cognizant of dynamic phenomena such as mechano-transduction and undifferentiated precursor cell signaling pathways.

However, acquisition of dynamic organogenesis reported by endogenous fluorophores is a challenging task owing to anisotropic contrast across the field of view (FOV). This is due to photon propagation through heterogeneous tissue (Teranikar et al., 2020, 2021). Hence, precise orchestration of in vivo volumes requires high sensitivity with respect to dynamic tissue motion and differing scales. As a result, optical aberrations often induce redundancy in the imaging sample space, affecting interpretability of feature attributes. Furthermore, cell studies are largely restricted to manual boundary demarcation due to the limited availability of binary classification methods impervious to heterogeneous contrast resolution (Astrakas and Argyropoulou, 2010; Marsh et al., 2018; Rajasekaran et al., 2016; Yin et al., 2014). Traditionally intensity-based segmentation methods such as the Otsu’s method, adaptive thresholding, isodata thresholding, and entropy-based thresholding have been used for automated cell tracking for their simplicity and speed (Goh et al., 2018). However, these methods are incapable of separating attenuated objects in closeness, proximity into meaningful biological regions (Xu et al., 2020). Another conventionally favored approach for biomedical image segmentation is the watershed algorithm (Beucher and Mathmatique, 2000; Koyuncu et al., 2012; Rajasekaran et al., 2016; Veta et al., 2013). However, the technique often causes over-segmentation or false detection of non-existent objects in dense tissue (Rajasekaran et al., 2016; Xu et al., 2020). Hence, there is a clear need for feature detectors impervious to low signal-to-noise ratio bioimages, for aiding accurate cell segmentation within the tissue architecture.

In this study, we propose the application of a scale space feature detector for isolating fused, myocardial nuclei blob morphology across distinct embryonic stages (Johnsen, 2000; Johnsen et al., 2011; Johnsen and Widder, 1999; Teranikar et al., 2020; Zhang et al., 2013). The proposed segmentation framework integrates Hessian difference of Gaussian (HDoG) feature detector with the watershed algorithm, for enhancing sensitivity of localizing individual nuclei. In combining these two algorithms, provides for seamless and straightforward segmentation with respect to background noise (Bharodiya and Gonsai, 2019). The proposed algorithm enabled in vivo characterization of wild-type zebrafish ventricular contractility and morphological traits such as nuclei number, area, and sphericity and in the myocardium.

Results

Isolating individual cardiomyocyte nuclei in dynamic zebrafish ventricular volumes

Distinguishing fused nuclei boundaries by manual contour segmentation or determining the intensity threshold for distorted nuclei outside the light sheet confocal parameter (focus region) is a complex task due to varying pixel intensities of overlapping nuclei (Figures 1A–1D). Furthermore, autofluorescence and dynamic cardiac motion convolutes lateral and axial imaging planes (Figures 1I and 1J). To separate nuclei clusters into disjoint regions, we implemented the difference of Gaussian (DoG) filter to enhance cell boundaries followed by the watershed algorithm to separate merged boundaries affected by anisotropic contrast. As a result, we were able to successfully quantify nuclei across different scales in a moderately dense cell environment at 48 h postfertilization (hpf) (Figures 1E–1H). More importantly, integration of the DoG scale space detector with the watershed algorithm enabled us to split longitudinally merged nuclei (Figures 1K–1L). In this regard, photon travel through heterogeneous tissue and restrictions imposed on resolvable sample depth are prone to induce sample de-focus (Figure 2A).

Figure 1

Isolating and segmenting cardiomyocyte nuclei from contracting heart using the DOG (Difference of Gaussian) filter in combination with the watershed algorithm

(A–D) 48 h postfertilization zebrafish ventricular volume was reconstructed using light sheet microscopy, in order to visualize time-dependent motion of myocardial cardiomyocyte nuclei. Raw volume comprised of fused nuclei clusters (yellow highlighted boxes), exacerbated by tissue scattering (B–D) Zoomed in regions demonstrate fused contours of nuclei, adversely affecting individual nuclei analysis (E) Ventricular volume was processed using the difference of Gaussian (DoG) edge detector in conjunction with the watershed algorithm to distinguish individual nuclei from adjacent neighbors.

(F–H) Zoomed in regions show successful separation of nuclei for aiding cell tracking and counting.

(I–J) 2D lateral and axial views illustrate tissue birefringence resulting in merging of nuclei longitudinally (K-L) Segmented lateral and axial views were reconstructed for qualitative assessment of contour separation of overlapping nuclei. (scale bar = 100 microns), a: atrium, v: ventricle.

Figure 2

Isolating individual nuclei volumes among high-density cardiomyocyte clusters, at distinct phases of ventricular contraction cycle

(A) Illustration depicting zebrafish ventricular myocardial nuclei sections, scanned by a Gaussian light sheet (blue solid line). There exists a tradeoff between the confocal parameter i.e. excitation lateral extent and beam waist (BW) i.e. light sheet axial resolution and hence, requires optimization of the Gaussian focus spot to effectively sample embryogenesis across different growth stages. The detection objective lens modulates effective field of view (FOV). Samples are scanned through the static optical section at discrete increments (dx) using mechanical transducers, to reconstruct complete in vivo 3days + time volumes from individual sections. Red arrow represents the blood flow direction of zebrafish heart.

(B and C) Raw systolic and diastolic nuclei reconstruction at 96 h (about 4 days) post fertilization, consisted of closely packed nuclei blobs as compared to 48 h postfertilization. Inaccurate nuclei localization is further exacerbated by dynamic contraction and relaxation.

(D and E) Application of difference of Gaussian (DoG) detector in conjunction with the watershed algorithm, exhibits reduced feature detection sensitivity leading to inaccurate reporting of nuclei number.

(F and G) Hessian DoG feature detector exhibits improved sensitivity to local affine transformations experienced by nuclei pixel neighborhoods during image acquisition. (scale bar = 50 micron), av: atrioventricular canal, v: ventricle, ot: outflow tract.

Integration of Hessian and difference of Gaussian (HDoG) to segment cardiomyocyte nuclei from dense environment

Compared to pinhole-based microscopy techniques, a potential cause of concern for LSFM modality manifests in the form of background contrast between adjacent cardiomyocytes (Figures 2B and 2C). As nuclei move dynamically across the field of view (FOV), undesired fluorescence emitted from fluorophore-binding sites outside the optical section beam waist affects accurate volumetric reconstruction (Figure 2A). Furthermore, intensity attenuation caused by low numerical aperture (NA) objectives aggravates poor signal-to-noise ratio (SNR).

Although we successfully separated longitudinally merged ventricular myocardial nuclei at 48 hpf using DoG feature detector, we encountered inaccurate nuclei number quantification beyond 72 hpf. We assume low pixel intensities produced by the DoG edge detector response prior binarization, resulted in under reporting of nuclei (Figures 2D and 2E). To overcome this, we applied the Hessian difference of Gaussian detector (HDoG) in combination with the watershed algorithm, for accurate contour separation and assessing the morphology of wild-type myocardial nuclei in vivo. The hessian determinant was used to localize saddle points (Marsh et al., 2018). Saddle points can be defined as neither an intensity maximum nor minimum, that represent connecting nuclei edges. This approach improved detection sensitivity in the presence of multiple intensity peaks for a single biomarker (Figures 2F and 2G). The segmented labels were further used for investigating nuclei shape and ventricular contractility, apart from nuclei counting.

Segmentation accuracy evaluation

Nuclei were detected for each distinct developmental phase: 48 hpf (Figures 3A–3C), 72 hpf (Figures 3D–3F), and 96 hpf (Figures 3G–3I), to compare segmentation robustness for sparse nuclei distribution at 48 hpf with respect to densely populated ventricle at 96 hpf. We used a segmentation ratio to evaluate segmentation accuracy, by comparing Hessian DoG nuclei images to manual nuclei segmentation, with respect to static 3D zebrafish heart confocal images. The segmentation ratio is the number of scale space segmented nuclei divided by the number of cardiomyocyte nuclei manually counted in the confocal images as ground truth. If the numerical value = 1, the segmentation is identical to the raw images. If the numerical value >1, there is over-segmentation in the segmented images. If the numerical value <1, there is under segmentation in the segmented images (Figure S1). Our analysis found that the ideal segmentation was repeated across developmental stages using the Hessian scale space (Figure S2).

Figure 3

Visualizing cmlc:GFPnuc zebrafish ventricular nuclei deformation at distinct developmental stages (48 – 96 h postfertilization), across the cardiac contraction cycle

(A–I) The Hessian DoG scale space representation was used for localizing cardiomyocyte nuclei ranging from different sizes, as a result of which we were able to assess ventricular contractility and complex nuclei morphology in vivo (scale bar for A-C = 100-micron, scale bar for D-I = 50 micron). A:atrium, v:ventricle.

Quantification of local contractility via tracking cardiomyocyte nuclei

After we processed the images to visualize individual nuclei, we performed contractility analysis by tracking cardiomyocyte nuclei across 48 to 120 hpf to quantify the local cardiac contractility based on this novel segmentation approach (Figures 4A–4G). We investigated the stretch level change of developing zebrafish heart and normalized the temporally changing stretch values for the innermost and outermost curvatures at each developmental stage (Figure S3). In addition to stretch, we calculated area ratio comparison between innermost curvature and outermost curvature areas. The area ratio is a description of local deformation of the area inside of three markers’ 2D stretch ratio. We analyzed the area ratio as a function of time, using three cardiomyocytes as markers. We found area ratio of the outermost curvature area, where the opposite side of the atrioventricular canal receiving blood directly from the atrium, has a higher area ratio than the innermost curvature area of the ventricle (Figures 4H–4J).

Figure 4

Selected markers utilized area ratio analysis

(A–F) represents the systolic reconstruction of ventricular myocytes at 48 hpf, 72 hpf, and 96 hpf, respectively, while (D–F) represents the diastolic reconstruction of myocytes at different developmental stages.

(G and H) Schematic illustrating the nuclei region of interest. Blue windows represent light sheet sections. Zebrafish ventricular volumes were sampled to compare the innermost curvature contractility (green markers), with respect to the outermost curvature (red markers) (H) Area ratio for innermost curvature by tracking three cardiomyocytes highlighted green in the blue optical plane, which elucidate increasing contractility trend observed across distinct developmental stages.

(I) The area ratio for outermost curvature calculated by tracking cardiomyocytes highlighted red in the blue optical plane, indicates the outermost curvature has higher contractility compared to the innermost curvature. (J) Outermost curvature has a significantly higher area ratio compared to innermost curvature after 72 hpf (n = 3, p = 0.05, one-tail t-test).

Quantifying zebrafish cardiomyocyte nuclei development

The average values for number of nuclei in a developing zebrafish heart were 159 ± 13, 222 ± 17, 260 ± 13, and 284 ± 10 for 48, 72, 96, and 120 hpf, respectively (Figure 5A) (n = 15). We observed cardiomyocyte nuclei in outermost curvature had larger systolic and diastolic volumes to innermost curvature area (Figures 5B–5E, Table S1). Hence, we assume the outermost ventricular curvature experiences higher mechanical deformation (Figures 4H and 4I) due to direct inflow of blood from AV canal (Figure 4G), resulting in larger nuclei volumes.

Figure 5

Zebrafish cardiomyocyte nuclei analysis

(A and B) We observed an increase in the number of ventricular cardiomyocyte nuclei for successive developmental stages. Asterisk denotes statistically significant difference with respect to previous time point. p ≤ 0.05 (B) Systolic and diastolic nuclei volume expansion observed for the inner curvature.

(C) Systolic and diastolic volume trends observed for the outer curvature. (D and E) Systolic and diastolic nuclei surface area growth observed for the inner curvature (E) Systolic and diastolic nuclei surface area observed for the outer curvature. n = 15.

Contractility effect on morphology of cardiomyocyte nuclei

Apart from area characteristics, we also quantified the circularity of myocardial ventricular nuclei (Figure 6). Isotropic/spherical nuclei in the innermost curvature region (Figures 6A and 6B) were evaluated to have an average elongation index of 0.91, while more elongated/ellipsoid nuclei in the outer curvature region nuclei had an average elongation index of 0.71, suggesting structural anisotropy. Interestingly, nuclei exhibit distinct eccentricity (major/minor axis ratio) according to their ventricular location, despite dynamic expansion and contraction across the cardiac cycle (Figures 5B–5E). In this regard, we observed spatially confined cardiomyocyte nuclei in the innermost curvature with shorter major and minor axis lengths, compared to larger outer curvature nuclei (Figure S4). Hence, we hypothesize that different cardiomyocyte shapes (Table S1) are modulated by varying contractility in different ventricular regions (Figures 4H and 4I). This is in accordance with elongated nuclei volumes for accommodating greater mechanical stress in the outer ventricular concave regions, as compared to smaller nuclei volumes in the inner convex region.

Figure 6

Systole vs diastole circularity analysis

(A) Inner curvature nuclei are observed to have a circular shape, (symmetric circle elongation = 1, ellipse <1) with slightly higher values observed for the diastole.

(B) Outer curvature cardiomyocyte nuclei are observed to have an elongated shape with higher elongation observed in the diastole.

(C and D) Volumetric reconstructions of the circular shape of inner curvature and elliptic shape of outer curvature myocytes were visually presented, respectively. In addition, the corresponding lateral and axial views are shown with binary images (scale bar = 15 um).

Discussion

Scale space theory can be understood as a hierarchal set of 2D images produced for each optical section, obtained by blurring the image from fine to coarser scale (Lindeberg, 1993; Marsh et al., 2018; Witkin, 1983). This results in suppression of all image objects equal to the size of the Gaussian kernel. Each defocused image contains a distinct number of edges obtained by blurring unresolved pixel subsets to a coarser resolution. Hence, enabling multiscale edge visualization without any knowledge of nuclei sizes a priori (Lindeberg, 1999, 2013).

As zebrafish myocardial nuclei length varies spatiotemporally across 2–6 microns (Figure S4), we propose the integration of scale space theory and watershed segmentation for robust scale-invariant edge detection. Utilizing the inherent de-focus adaptation ability of the HDoG blob detector, we successfully isolated individual centers of mass (Videos S1, Video S2) for tracking dynamic cardiomyocyte nuclei (Video S3, Video S4). As a result, we were able to successfully characterize ventricular myocardial stretch post AV valve specification to heart maturation (Kula-Alwar et al., 2021; Miura and Yelon, 2011). Although transparency was induced in zebrafish using PTU, tissue birefringence (RI∼1.3–1.5 (Jing et al., 2018)) results in changes in optical path lengths of emitted photons (Johnsen, 2000; Johnsen et al., 2011; Johnsen and Widder, 1999; Teranikar et al., 2020). Consequently, the light scatter compromises the optical modality penetration capability, resulting in fusing of nuclei situated outside the confocal region (Figure 2A) (Teranikar et al., 2020). Taking this into consideration, we sought to design an automated blob detection framework that provides high sensitivity and repeatability for a singular Gaussian intensity peak detection corresponding to each nuclei centroid. Furthermore, the proposed framework enabled in vivo quantification of morphological descriptors such as nuclei volume, surface area, and shape.

Wild-type Tg(cmlc2:egfp) zebrafish have been observed to exhibit cuboidal cardiomyocyte morphology in linear heart tube (24 hpf) and IC ventricular myocardium, with respect to elongated cardiomyocytes in the OC (Auman et al., 2007; Kula-Alwar et al., 2021; Miura and Yelon, 2011). Studies indicate regionally distinct cardiomyocyte phenotypes such as cell count, area, or sphericity are regulated by mechanical stimuli such as contractility or blood flow during heart maturation (Auman et al., 2007; Miura and Yelon, 2011). This has been validated through the mutation phenotype half-hearted (haf) mutation lacking ventricular contractility. The haf mutant exhibited elongated cardiomyocytes with increased surface area across different parts of the ventricle including IC, resulting in a distended ventricle (Auman et al., 2007). Interestingly, cardiomyocyte count was observed to be consistent between the haf mutant and wild-type zebrafish ventricle, suggesting contractility is responsible for moderating the aberrant elongation of cardiomyocytes and not anomalous proliferation (Auman et al., 2007). Furthermore, previously performed studies (Auman et al., 2007) indicate cardiomyocyte number reflects a sigmoidal growth trend (Figure 5), subsequently plateauing at later stages (>96 hpf) thereby signaling specification into the myocardium. In this regard, the proposed feature detector and cell tracking algorithm can prove extremely beneficial for gaining insights into the effects of cardiac contraction on reducing proliferation and its secondary effects on cardiomyocyte morphology in zebrafish. Unfortunately, currently we cannot conclude that contraction is key to reducing the proliferation of cardiomyocytes due to lack of statistically significant data. However, the intensity of mechanical workload experienced by cardiomyocytes in different parts of the ventricle appears to be a regulatory mechanism for maturation into distinct shapes (Figures 4, Figure 6).

Analyzing the cardiomyocyte motion (Figure S3), we quantified the outermost curvature has higher area ratio than the innermost curvature (Figures 4H and 4I), thereby experiencing greater mechanical workload. Furthermore, we observed elongated cardiomyocyte nuclei morphology in the OC with respect to spherical morphology in the IC. Although no phenotyping screening of nuclei morphology has been performed with respect to modulation of contractility, our data indicate that cardiomyocyte nuclei shape and size corresponds to deformation experienced by distinct ventricular regions. Elongated nuclei (Figure 6) in the OC suggest larger surface area is required to accommodate greater OC mechanical intensity. On the other hand, IC consists of smaller, cuboidal nuclei due to lesser deformation compared to OC. In the developmental biology aspect, researchers primarily focus on the ventricular OC where trabeculae form, but there is lack of well-documented research regarding lack of trabeculae in the IC. Thus, the question remains how different biomechanical or molecular signaling engenders a trabeculated OC and smooth IC. Our study has the potential to elucidate ventricular development in zebrafish orthologs, and aid cardiac pathophysiology diagnosis or clinical translational of cardiac regeneration for pediatric population. However, further investigations will be required to validate this assumption. In this regard, nuclear morphology observed in cardiomyocytes isolated from neonatal rat ventricles reports similar findings, regarding systolic and diastolic heterogeneous cross-sectional surface areas due to deformation experienced by the cardiac cycle (Bray et al., 2010). Hence, our novel study provides exciting avenues to characterize cell count, morphology, and intercellular forces that may be responsible for cardiomyopathy in humans. Future studies will involve modulation of contractility to characterize cardiomyocyte morphology in the IC and OC.

In summary, we have presented a scale-invariant feature detector for quantifying individual morphological characteristics of merged nuclei and biomechanical analysis of the zebrafish ventricle. Our proposed blob detection and cell tracking approach will prove to be extremely beneficial for analyzing cell count, volume, area, sphericity, proliferation, or cardiac function for characterization of cardiomyopathy phenotypes.

Conclusion

In this report, we were able to successfully interrogate dynamic zebrafish cardiac tissue non-invasively using bona fide biomarkers such as cell elongation, volume, and surface area. Moreover, we quantified the number of cells and the mechanical workload experienced by the ventricular inflow and outflow regions during the systole and diastole, respectively.

Limitations of the study

Although we successfully separated merged nuclei clusters across varying scales and densities, the reproducibility of the Hessian DoG feature detector is highly dependent on appropriate identification of Gaussian blurring weights (Figure S1). Moreover, Hessian scale space detector followed by watershed postprocessing is more prone to over-segmentation with higher variability in nuclei count, in comparison to DoG feature detection (Figure S1) if kernel weights are not selected appropriately. On the other hand, DoG scale space detector is inherently prone to erosion of boundaries due to bandpass operation, resulting in reduction of nuclei volumes and the object area affecting quantification. Other modality limitations include absence of peripheral nuclei during diastole that may be present during the systole, due to ballooning of ventricle outside the light sheet confocal region. As in vivo cardiomyocyte cell tracking and counting requires invariancy to sample translation without distortion in shape, the Hessian DoG operation was effectively used to localize individual nuclei based on pixel intensity gradients.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Dr Juhyun Lee (juhyun.lee@uta.edu).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

The animal experiments were performed in agreement with the UT Arlington Institutional Animal Care and Use Committee (IACUC) protocol (#A17.014). The transgenic zebrafish line used in this particular study is the Tg(cmlc2:nucGFP), with the cardiomyocyte nuclei labeled with GFP (Green Fluorescent Protein). The zebrafish embryos were maintained at 28.5°C in system water at the UT Arlington Aquatic Animal Core Facility. 0.0025% 1-phenyl 2-thoiurea was added to the embryo medium starting around 20–24 hpf to suppress pigmentation. Prior to imaging, embryos were anesthetized in 0.05% tricaine (MS 222, E10521, Sigma-Aldrich, St-Louis, MO) to avoid sample movement. Sex determination and segregation of zebrafish was not performed in the embryonic stage (48–120 h postfertilization).

Method details

Light sheet microscope (LSFM) implementation

Our home built light sheet microscope consists of a single-side excitation pathway and a custom water dipping lens (20x/0.5 NA UMPlanFL N, Olympus, Tokyo, Japan) detection. In the illumination pathway, a cylindrical lens (LJ1695RM, Thorlabs) coupled with a 4x objective lens (4X Plan Apochromat Plan N, Olympus, Tokyo, Japan), are used to collimate a cylindrical light-sheet with ∼4-5-micron thickness. Furthermore, a mechanical slit aperture (VA100C, Thorlabs) is modulated across distinct developmental stages (48-,72-,96- and 120 hpf), to accommodate ventricular circumferential extent across the light sheet confocal region (Figure 2). A DC servo motor actuator (Z825B, Thorlabs) is used for sample translation in the axial direction (z-step velocity and acceleration = 0.005 mm/s). The optical detection pathway consisting of the water lens, infinity corrected tube lens (TTL 180-A, Thorlabs) and sCMOS camera (ORCA flash 4.0, Hamamatsu, Japan, camera pixel size = [6.5 (um)ˆ2] = [6.5/20x = 0.325 um], camera exposure time = 30–50 ms), is used for non-gated 4D (3D + time) cardiac volume acquisition.

As the zebrafish ventricle undergoes periodic deformation during peak systole to end diastole, optical sections were acquired at varying depths in the sample, covering 4–5 cardiac cycles20,23. Since triggering of image slices is not synchronized to a particular phase in the cardiac cycle, we performed volumetric reconstruction a posteriori to ensure alignment of adjacent optical sections. For this purpose, we estimated the period of each individual cycle by minimization of the least squares intensity difference criterion and calculated the relative period shift to ensure synchronization between independent cardiac cycles20,23

Preparation of zebrafish for assessing cardiac function

The animal experiments were performed in agreement with the UT Arlington Institutional Animal Care and Use Committee (IACUC) protocol (#A17.014). The transgenic zebrafish line used in this particular study is the Tg(cmlc2:nucGFP), with the cardiomyocyte nuclei labeled with GFP (Green Fluorescent Protein)19. The zebrafish embryos were maintained at 28.5°C in system water at the UT Arlington Aquatic Animal Core Facility. 0.0025% 1-phenyl 2-thoiurea was added to the embryo medium starting at 20–24 hpf 4,32 to suppress pigmentation. Prior to imaging, embryos were anesthetized in 0.05% tricaine (MS 222, E10521, Sigma-Aldrich, St-Louis, MO) to avoid sample movement. Upon administering the anesthetic, alive embryos were embedded in 0.5% low-melt agarose gel inside a fluorinated ethylene propylene (FEP) tube (1677L, IDEX, Chicago,IL). Furthermore, the FEP tube was suspended in water within a custom 3days printed ABS (Acrylonitrile Butadiene Styrene) cuvette (designed using solid works) housing the water dipping lens, to ensure near isotropic refractive index between the water dipping lens and sample inside the tube. (Refractive index of water = 1.33, refractive index of agarose and FEP tube = 1.34). Refractive index matching is necessary to avoid distortions and intensity attenuation in the optical sections13.

Image processing framework

Haze removal using the dark channel prior (DCP) method

Introduction of haze by the ambient medium or scattering due to particulate matter, degrades the performance of computer vision tasks(Lee et al., 2016; Teranikar et al., 2020). A haze free image can be retrieved by using the image degradation model based on the Dark Channel Prior (DCP) algorithm.

(Lee et al., 2016; Teranikar et al., 2020)

where I(x) is the degraded image, J(x) is the original irradiance captured by the CMOS camera, t(x) represents the scene depth and A is the scattering introduced by the ambient light. Using the dehazing algorithm, we estimated the intensity transmission map t(x) using the imreducehaze() MATLAB function (Teranikar et al., 2020).

(Lee et al., 2016)

where represents the scattering coefficient and d represents the scene depth. We used the estimated intensity transmission map as a preprocessing step before performing the DoG operation. By estimating the contrast attenuation with respect to distance, we were able to emphasize edges.

Intensity maxima localization at nuclei centers using the difference of Gaussian (DoG) filter

The DoG filter can be effectively used to enhance edge visualization for images suffering from poor contrast. In this study, the greyscale bandpass operation is performed by subtracting a blurred version of the transmission estimate from a lesser blurred version of itself,where g1(x) and g2(x) are the Gaussian kernels having different standard deviations. Using the DoG filter, we were able to localize blobs to nuclei centers by isolating spatial frequencies correlating to the Gaussian illumination maxima.

Precise contour delineation using the hessian scale space representation and watershed algorithm

The hessian scale space representation can be described by the convolution:

(Marsh et al., 2018; Rajasekaran et al., 2016)

Where D(x,y,t) represents the family of images, derived from the original image. t represents the degree of blurring. Hence, Equation (4) can be described as the convolution of the DoG blob maxima image with the hessian blob detector Gaussian blur kernel G(x,y) at different degrees of blur (t > 0). The blurring scale selection was based on the ratio t +1 = r ∗t (Marsh et al., 2018), where r is a constant.

The workflow involved for the hessian blob involves (Marsh et al., 2018),(Rajasekaran et al., 2016),

Computing the absolute magnitude of the intensity gradient image obtained by convolving the DoG bandpass image with the derivative of Gaussian filter.

Computing the double derivative of the absolute magnitude image calculated in the previous step

Imposing boundary conditions on the hessian determinant value [det D(x,y,t) < 0] (Rajasekaran et al., 2016) at every pixel, for indicating saddle points.

The image arithmetic operation (OR – operation) results in the union of the DoG localized intensity maxima and contour information from the Hessian blob, aiding the successful splitting of nuclei.

Preprocessing strategies

Images corrupted by noise or tissue scatter, were filtered by using a Gaussian kernel with an appropriate SD followed by the background subtract operation in ImageJ. In addition, image processing code is enclosed (Data S2, related to Figure 2).

Cell counting and area measurements

After converting raw optical images to binary images, we performed to count cardiomyocyte nuclei and their area analysis by using the 3D object counter plugin33 in ImageJ. The plugin can be accessed by: ImageJ – Analyze – 3D Object Counter. After cropping the ROI (ventricle in this case), we used the plugin to quantify number of object voxels (volume), surface voxels of individual nuclei volumes and the number of 3D nuclei objects in the ventricular stack. The plugin can also be used to retrieve the centroid geometric coordinates of object volumes. The user is required configure 2 important parameters namely, (a) intensity threshold to separate background and foreground pixel populations and (b) size threshold to exclude smaller objects from the analysis. The plugin allows user to configure object counting based on the presence or absence of touching edges.

Cardiac myocyte nuclei tracking

We utilized the segmented, processed, and time synchronized images to reconstruct three-dimensional volumes through time for a cardiac cycle to perform this tracking. We then passed these images through a custom MATLAB (Mathworks) code to perform for key steps (Data S3, related to Figure 4).

This MATLAB code performed following 4 steps.

The code compiles images into easily searchable 4D matrices.

The code resolves the 4D matrices of segmented images into centers of mass based on high pixel concentration areas for each time step.

The user selects three markers to represent our plane for stretch calculations.

The code searches through the 4D stack of centers of mass to determine the closest center of mass in the next time step and stores these points in a matrix of position values.

Each stored triplet value is the x, y, and z position of a particular nucleus at a particular time. This format is easily searchable and allows for a multitude of calculations. This code assumes that there can be no erratic motion of the nucleus with a high enough sampling frequency. The location at each time step depends on the prior location. Imaging with a high sampling frequency supply data that meets this assumption requirement. Other works have utilized similar works, including Meijerling et al. (Meijering et al., 2009).Drawbacks of this method include the requirement for user interaction. To verify that the cell tracking occurs appropriately, the user must analyze each vector to ensure the vector does not violate the small motion assumption. This process can become time-consuming and increases the chance of human error. Subsequent work can expand and refine this cell tracking method to include other parameters, including a probability net for machine learning applications and size and orientation to decrease ambiguity and reduce the user input requirement.

Contractility analysis

We selected and tracked three cardiomyocyte nuclei for both the innermost and outermost curvature. After tracking the location of three cardiomyocytes through each time instance, we utilized the following method to determine the deformation gradient with normalized one cardiac cycle as 0.5 s starting from ventricular end-systolic stage. We determined the stretch ratio at each time instance into principal stretch values reported as λ1 and λ2, or the longitudinal and circumferential principal direction followed by previous methods41,42. These principal stretch vectors correspond to the first and second principal strain directions. When viewed on Mohr’s circle, they correspond to the maximum and minimum normal strain values where the shear strain is resolved to zero (Figure S6). These values are represented in the Cartesian coordinate system as the x and y direction or in polar coordinates as the zero-degree rotation and 90-degree rotation. We established the area ratio by multiplying the two principal stretch values. Area ratio provides a description of the total in plane deformation from the initial undeformed state which was selected as the start of filling.

Quantification and statistical analysis

For statistical analysis, we performed ad hoc pairwise comparisons for three morphological parameters to characterize the maturity of the heart (p value = 0.05)10. We analyzed the number of visible nuclei, the total volume, and total surface area. We estimated each of these parameters using built in functions in ImageJ (NIH, Bethesda, MD) with n = 15. Additionally, we cleaned the data in excel utilizing Chauvenet’s criterion to determine which values were outliers and should be removed. After removing outliers and cleaning the datasets in excel to reduce the chance of error due to our sampling technique, we compared the data with one-way ANOVA. If we detected a statistically significant difference for any comparison, we performed Tukey’s test for multiple comparison of means. This test inherently compensates for multiple comparisons, which allowed us to use an alpha value of 0.5. All values herein are reported as mean +/− standard deviation in the figures and respective figure legends.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, peptides, and recombinant proteins
0.0025% 1-phenyl 2-thoiurea	Sigma-Aldrich	P7629
0.05% tricaine (MS 222)	Sigma-Aldrich	E10521
Experimental models: Organisms/strains
Zebrafish: Tg(cmlc2:nucGFP)	Sharpe M et al.	Gifted by Dr Barnes at Boston children’s hospital, Harvard Medical.
Software and algorithms
ImageJ	Schneider et al., 2012	https://imagej.nih.gov/ij/
Hessian Determinant plugin	Sato, Y. et al.	https://imagescience.org/meijering/software/featurej/
Other
Cardiomyocyte nuclei tracking code	This paper	Chuong CJ, Sacks MS, Templeton G, Schwiep F, Johnson RL Jr. Regional deformation, and contractile function in canine right ventricular free wall. Am J Physiol. 1991;260(4 Pt 2):H1224-H1235. https://doi.org/10.1152/ajpheart.1991.260.4.H1224