Introduction

The COVID-19 pandemic demanded the health care system as well as the scientific community to react fast to reduce the spreading and mortality of the disease. Research investigating therapeutics and vaccines has been set in motion with unprecedented speed and the US Food and Drug Administration (FDA) launched several Emergency Use Authorization announcements. The Expanded Access (EA) Program gained importance to promote emergency use of unapproved agents to treat COVID-19 outside of clinical trials. In this regard, drug repurposing became an important research pillar. Indeed, a more general therapeutic treatment beyond vaccination will remain of great clinical importance given the observations that vaccinations provide immunity of a limited duration1 and that newly arising SARS-CoV-2 variants can escape to some extent or completely the immune response obtained through vaccination2,3. One important caveat of drug repurposing is that even though the compounds were previously approved for treatment of other diseases or are under clinical testing, the impact of many of these agents on embryonic development is yet unclear, limiting the access of treatment to pregnant COVID-19 patients.

The zebrafish plays an important role in phenotype-based screening in the field of embryonic and cardiovascular development, drug discovery, and toxicity screenings4,5,6,7,8. In comparison to in vitro, cell-based systems, the zebrafish allows to test drugs in the context of an intact vertebrate organism9. With its small size, transparency, rapid external development and high degree of genetic similarity to humans it allows fast assessment of how drugs effect defined developmental stages at tissue and cellular resolution.

The zebrafish heart starts to beat at 25 h post-fertilization (hpf). The vasculature is initially formed by the dorsal aorta and the cardinal vein from which intersegmental vessels (ISV) sprout at the somite boundaries. The sprouting of these endothelial cells can be used as a proxy for proper angiogenesis and this process is controlled by signaling cascades that are conserved in humans10. Fluorescent myocardial and endothelial reporter lines make imaging straight-forward and scalable when combined with medium- to high-throughput microscopy in 96 well-plates on multiple resolution levels using template-matching algorithms11,12. For example, such a setup has been successfully used for high-throughput identification of small molecules that affect human embryonic vascular development with direct applicability to cancer treatment13,14. Furthermore, the zebrafish was used to identify compounds with cardio-protective function in humans15 and small molecules that promote cardiomyocyte proliferation also found to be effective on mouse cardiomyocytes16. Furthermore, studying the swimming behavior of zebrafish can provide insights into sensory and motor function of developing embryo, and this approach can also be adapted for drug screening in a high-throughput setup17,18.

Here we used the zebrafish model organism to screen for possible side effects on cardiovascular and embryonic development of 162 marketed drugs that are either already used in clinics or being considered as potential treatments for COVID-19. While the effect of some compounds has been previously studied individually in zebrafish, a systematic and comprehensive analysis allowing for comparison between compounds was lacking. We included the compound Sabizabulin not yet tested in the zebrafish model and currently used as FDA-approved drugs for COVID-19 treatment. We implemented a screening pipeline that includes U-net architectures for automatic image segmentation and analysis, which can be easily scaled up for additional compounds as they become available. The results were uploaded on a straightforward web-based interface that is easily accessible and allow further mining of results. Additionally, we show that the compounds can be efficiently tested for their ability to modulate immune signaling in upon SARS-CoV-2 spike treatment. Our results not only have important implications for the treatment of COVID-19, but also provide valuable insight into the unintended effects of FDA-approved drugs on early development and cardiovascular system formation. Moreover, our pipeline can be adapted for determining the possible effects of further compounds on embryonic development or the cardiovascular system.

Results

Establishing a quantitative screening platform to identify drug effects on embryonic and cardiovascular development

One of the earliest initiatives to promote research on drug repurposing to fight COVID-19 was supported by the Medicine For Malaria Venture (MMV). The consortium assembled and freely distributed the Covid Box for research, comprising a drug library of 160 compounds with known or predicted effects on one or multiple steps of SARS-CoV-2 infection or COVID-19 disease progression. The majority of compounds in the Covid Box are categorized as anti-infective agents, that would target biological processes of pathogens (56 agents; 35% of the library), all the other types would target biological processes in the host (Fig. 1a,b). As an extension of our screening efforts, we also purchased Molnupiravir and Sabizabulin (from MedChemExpress), which are approved agents for the treatment of COVID-19 that successfully inhibit viral replication and transmission from cell to cell, respectively19,20.

Figure 1
figure 1

Overview of drugs included in zebrafish cardiovascular and behavioral screening setup. (a) Information on type of candidate drugs used. All primary drug indications are listed and a set of anti-infective agents is shown. (b) A representative table of drugs that are FDA-approved or in III Clinical trial with the indicated Cmax blood concentrations and risk during pregnancy.

We treated zebrafish embryos at 28 hpf with the compounds dissolved in DMSO (Fig. 2a). The baseline concentration for the screening of the drugs was 1 µM, as recommended by the MMV, which was approximately 10 × lower than used commonly in zebrafish screens, e.g.21,22,23. For those compounds also used in the clinical context for COVID-19 treatment we used compounds at Cmax (Fig. 1b). Treatment with each of 162 compounds on 10 zebrafish embryos was replicated two times or more. For Ivermectin and Sabizabulin, as 1 µM concentration elicited severe toxic effects, we lowered the dose. Treated embryos were then subjected to two parallel imaging platforms, that allowed us to observe embryonic and cardiovascular development on the one hand, and behavior on the other. We analyzed the data in a combined workflow in semi-automatic as well as automatic manner to obtain quantitative data for each drug (Fig. 2).

Figure 2
figure 2

Schematic representation of the experimental workflow. (a) Tg(fli1:EGFP);(myl7:mRFP) embryos were used for the morphological assay (green: cytoplasmic GFP in endothelium; magenta: membrane RFP in cardiomyocytes). After drug treatment, between 96 and 104 hpf, the anesthetized larvae were transferred to 96 well plates with pre-formed agarose beds by a 3D printed mold. This is followed by image acquisition. Treated larvae positioned in a 96 well plate are imaged using a Smart Imaging Wide field Fluorescent Microscope Platform. (b) Image analysis. (1) Image analysis through a semiautomatic and blinded workflow based on a Fiji-macro. (2) To count the intersegmental vessels (ISV) the user draws a line (red line) along all ISVs. A kymograph along the Z-axis allows detection of local maxima. (3) The larval length from anterior to posterior (AP) is measured by drawing a line from the most anterior of the head to the most posterior of the tail. (4) Effusions are identified as shown in the examples. (5) For heartbeat measurements, hearts are detected with an automated threshold to generate a mask and a kymograph along the time axis detects the local maxima. For the deep learning-based analysis, image data of the ISV and the heart are manually annotated. The annotated data are used to train two different U-Net models. The ISV-Seg model is applied to segment each ISV and calculate the area. The Heart-Seg model is used to segment ventricle and atrium. The minima and maxima of the ventricular area along the time axis are used to calculate the ejection fraction. For the behavioral assay, wild type larvae were treated with the drugs and transferred to 96 well plates. The swimming of larvae with a defined light–dark exposure was recorded and the tracking data exported. (c) The results are collected and used for meta-analysis with Python. The results were presented as heat maps and in a customized online Streamlit data app.

To observe morphological and cardiovascular features, we used the homozygous double transgenic line Tg(fli1:EGFP)y1;(myl7:mRFP)ko08 in which the fli1a promotor controls expression of the enhanced green fluorescent protein encoding gene (EGFP) and the myl7 promoter drives expression of the membrane-tagged red fluorescent protein encoding gene (mRFP)24,25. At 4 days post fertilization (dpf), larvae were transferred to 96-well plates, prepared with beds of low-melting agarose using a 3D-printed mold and imaged using a High-content Smart Imaging Fluorescent Microscope. This was followed by an image analysis pipeline (Fig. 2b) and metaanalysis and visualization (Fig. 2c).

On the other hand, we evaluated swimming behavior, an established method for testing neuroactive compounds (Basnet et al.), using DanioVision™ recording chamber. Tracking the behavior of zebrafish larvae is based on switching between bright and dark phases and can reveal anxiety, vision impairment, muscular weakness and reactivity18 (Fig. 2, Supplementary Fig. S1). These data were also integrated into the metaanalysis and web-based visualization platform.

A third of tested compounds altered cardiovascular development, embryonic growth, or behavior in the zebrafish model

First, we evaluated the survival and pericardial effusion of the embryos during treatment for both assays (Fig. 3a and Supplementary Fig. S2). The mortality rate (MR) was calculated as the average percentage of the larvae (n = 20–60 larvae per drug) that died before imaging in the respective assay. We observed that 7 out of 162 drugs (4.3%) caused a MR over 10% at 1 µM, including Apilimod, Astemizole, Ivermectin, Manidipine, Midostaurin, Niclosamide and Pimozide. Moreover, 18 drugs (about 12.5%) caused pericardial effusion in a quarter or more of the treated larvae. Among them, Astemizole (94.4%), Pimozide (100%) and Ponatinib (100%) were the compounds most frequently leading to pericardial effusion.

Figure 3
figure 3

Summary of results obtained for the compounds screening. (a) Bar plot representing the percentage of mortality and pericardial effusion of embryos tested in morphological assay. (b) Bar plot representing percentage of compounds leading to statistically significant alterations in the analyzed parameters. (c) Heat map of the results obtained for the tested compounds. Compound names are listed on the right of the heat map and are allocated to groups: ΔHeart, ΔBody length, ΔVasculature, ΔActivity or “No significant effect”, depending on the results. E.g. a compound that led to high reduction of ejection fraction compared to controls would appear in the group ΔHeart. The formula used to calculate the scores in the heatmap is given at the bottom. The p-value is derived from the Mann–Whitney U test treatment vs. control, the lower the p-value the higher the –log10(p-value). The hedge’s effect size gives an estimation of how strong an effect is. The fold change (FC) shows whether the values are higher or lower than the control and how many times compared to the control. The scores are shown in a relative color scale: green, most positive score; magenta, most negative score; grey, no significance using a Mann–Whitney U test. Shown are the mean values from at least two technical replicates, each replicate is composed of at least 10 embryos per condition. o symbol indicates novel statistically significant observations in behavioral assay. All the compounds are used at 1 μM except when it is mentioned otherwise. Asterisk in Ivermectin refers to the fact that it was used at 1 μM for the morphological assay and at 0.5 μM for the behavioral assay. Compounds mentioned throughout the manuscript are highlighted in bold.

We measured the body length and several parameters to assess cardiovascular development and behavior of larvae and clustered them according to their effect on body length, heart function (rate and ejection fraction), vasculature formation (number and area of ISV) and effect on activity (moving during accommodation, velocity under bright light and dark conditions as well as the ratio of velocity between dark and bright phases). 20% (33/162) of the drugs significantly altered the size of the larvae, 56 drugs (34.5%) altered the heart rate and 33 drugs (20%) affected the ejection fraction (Fig. 3b). For assessment of behavior in the presence of the compounds, we scored three phases of the bright/dark locomotion test—accommodation, bright and dark phases—and calculated the ratio between swimming velocity in the dark and bright. 27 drugs affected velocity in the dark, 30 in the light phase and 27 drug treatments showed alterations in activity during accommodation. To understand whether the difference between activity in the bright and dark phases was proportional to the control, we calculated the logarithmic bright/dark ratio. The bright/dark ratio was altered in 10.5% (17/162) of the drug treatments.

We clustered compounds according to the parameter affects and also generated a final cluster by those compounds that at the tested concentrations did not produce any deviation in the analyzed parameters when compared to controls (Fig. 3c).

We looked closer to those compounds leading to more pronounced defects on development (Figs. 4 and 5). Astemizole, Pimozide, Ponatinib, and Ivermectin led to the highest reduction in body length. While these compounds had been studied in the zebrafish27,28,29,30, size reduction had not been documented. Furthermore, 17α-Hydroxyprogesterone caused a previously not described larval overgrowth. Some of the observed effects on heart rate and ejection fraction had previously been reported, as was the case for Manidipine and Astemizole31,32. However, neither the increase in ejection fraction by Tacrolimus nor the effects of Ravuconazole (an antifungal agent), GSK-369796 (an antimalarial compound), Amuvatinib and Regorafenib have to our knowledge not been reported before. The Tyrosine kinase inhibitors Regorafenib, Cabozantinib, and Sorafenib revealed a strong inhibition of ISV formation, as reported before, and thus served as a positive control for the experimental set up33,34,35. However, the negative effects on ISV area by the compounds Astemizole and Pimozide had not been reported before.

Figure 4
figure 4

Combined boxplots and swarmplots showing the morphology and activity measurements for each larva. (a) Results from treatments using drugs studied in the context of Covid-19 research. (b) Results from treatments leading to strongest phenotypic alterations. (c) Results from treatments using COVID-19 drugs used at Cmax. The asterisks represent the p-values from the Mann–Whitney U test * < 0.05; ** < 0.01;*** < 0.001. Results from control larvae highlighted with red fill color. The red line indicates the median of the control as reference for treatment conditions. The plots are shown for all performed measurements.

Figure 5
figure 5

Selection of specific compounds eliciting effects on zebrafish embryonic development. (a) Results of 11 drug treatments leading to strong phenotypic alterations. Heatmap showing the highest positive score in dark green and the lowest negative score in dark magenta. Each measurement is individually scaled. Grey color indicates no significance according to the Mann–Whitney U test. Shown are the mean values from at least two technical replicates, each replicate is composed of at least 10 embryos per condition. Boxplots with individual data points for the morphology assay can be found in Fig. 4. (b) Example images of the control and treatments with visible effects are shown in brightfield, vasculature Tg(fli1:EGFP) and heart Tg(myl7:mRFP). Brightfield and GFP images are shown as sobel-projections. The dynamic range of GFP images was homogenized using the DEVILS Fiji-plugin. Scale bars: 250 μm. (c) Selected compounds at clinically relevant concentrations result in mild modulation of embryonic. A heatmap (left) showing percentage of pericardial effusion and mortality of the selected drugs at their clinically relevant concentrations. Mortality rates for morphological assay (with PTU) and behavioral assay (without PTU) are plotted separately followed by a column of mean mortality rate. Dark magenta shows the highest effect. The relative percentage effect of these selected drugs at their clinically relevant concentrations are shown in the heatmap (right). The highest positive effect is shown in dark green and the lowest negative effect is dark magenta. Grey color indicates no significance according to the Mann–Whitney U test. Shown are the percentage from median values from at three technical replicates for drugs except for Remdesivir, that was tested in two replicates, each replicate is composed of at least 10 embryos per condition during drug treatment.

Of all compounds tested, Ivermectin treatment also caused the most severe effects on swimming behavior, consistent with previous findings26. The effect on swimming behavior of Ponatinib, Regorafenib, Cabozantinib, Ravuconazole, Delanzomib, Oxyclozanide, Anidulafungin, Camostat, Amuvatinib, Hanfangchin B, Bemcentinib and Apilimod and Molnupiravir had to our knowledge not been reported before.

Effect of compounds studied in the context of COVID-19

Next, we focused at the impact of those compounds with the highest interest for COVID-19 treatment on cardiovascular development and larval activity in zebrafish.

The compounds Ivermectin and Hydroxychloroquine were repurposed for COVID-19 treatment in the initial stage of the pandemic but their use was discontinued36,37,38. Consistent with previous reports in zebrafish and the reported side effects in human patients30,36. Ivermectin was highly toxic also in this study. At 1 µM embryos died between 4.5 and 5 dpf and at 0.5 µM mortality disappeared, we still observed effects on larval behavior (Figs. 3b,4 and 5). (Figs. 3c, 4, 5a and Supplementary Figs. S2). Hydroxychloroquine treatment showed significant reductions of ejection fraction and ISV area (Figs. 3c, 4 and 5). The effects of Hydroxychloroquine are in line with mild cardiac phenotypes, observed in neonates after treatment during pregnancy and increased cardiovascular mortality39.

Clinical trials on Favipavir, Ribavirin, Umifenovir and Lopinavir have been started40 but are not FDA approved for COVID-19 treatment. Favipiravir had a mild effect on body length and Lopinavir reduced ejection fraction (Figs. 3c, 4). We found that Ribavirin and Umifenovir treatment mildly reduced the ISV area (Figs. 3c, 4). While the effects observed for these compounds showed statistical significance, the amplitude of these effects was very low.

Remdesivir, Molnupiravir19, Paxlovid (Nirmatrelvir and Ritonavir)41, and Baricitinib are all used in the clinics and were approved by the FDA for COVID-19 treatment42,43. Sabizabulin has also been recommended for COVID-19 treatment20. For these compounds we carried out a further experimental round at concentrations with higher clinical relevance, i.e. at the highest concentration detected in blood plasma in humans (Cmax)44,45,46,47 (Figs. 4c and 5c). Sabizabulin treatment still resulted in a 20% mortality at its C max (0.2 μM) (Fig. 5c). Those larvae that survived, did however not present any cardiovascular or behavioral defects (Fig. 5c). At Cmax (4.5 μM) Remdesivir led to decrease in heart rate, ejection fraction and in body size of the embryo (Fig. 5c). For Ritonavir, Cmax was considerably higher than the tested concentrations (15.25 μM). We found that this altered the ejection fraction and heart rate (Fig. 4b). Baricitinib at C max (0.144 μM) led to 18.3% mortality and in those embryos that survived caused alterations in heart rate and motility (Fig. 5c). Despite a Cmax of only 0.04 μM, at this concentration Molnupiravir still altered swimming behavior and a decreased body size.

In sum, at clinically relevant concentrations, all five compounds used in the clinical context led to alterations in embryonic development.

Drug candidates modulate zebrafish immune response to SARS-CoV-2 spike protein treatment

Zebrafish possess many similarities in the innate immune response to that of mammals48. While successful virus amplification was not observed in wild type strains, SARS-CoV-2 spike treatment causes temporal immune response in zebrafish embryos and adult fish49,50,51. We decided to explore the effect of drugs selected for Cmax analysis on the modulation of the inflammatory and immune response gene expression profile (Fig. 6a). First, we subjected 5 dpf larvae to recombinant SARS-CoV-2 spike protein and extracted RNA from the larvae and performed qRT-PCR using a marker panel of genes involved in spike protein response. We found that Spike protein treatment altered expression of several but not all of the selected marker genes (Fig. 6b). Next, we repeated the treatments but now included a group to which we added either Remdesivir, Ritonavir, Molnupiravir, Baricitinib, Sabizabulin or DMSO for 48 h (Fig. 6c).

Figure 6
figure 6

qPCR results for inflammation associated gene expression in zebrafish larvae (a) Scheme of larvae treatment with SARS-CoV-2 spike protein and selected drugs. (b) qPCR results for zebrafish larvae treated with SPIKE protein. (c) qPCR results represented as a heatmap representing fold change with respect to treatment with DMSO without spike (− SPIKE + DMSO). The upregulation is shown in dark green and the downregulation in dark magenta. Significance calculated between DMSO with SPIKE-protein treated group (+ SPIKE + DMSO) and the rest. Each color box represents a mean of 6 biological replicates of 5 embryo pools (n = 6, Mann–Whitney U-test treatment vs. + spike + DMSO; *p < 0.05, **p < 0.01, ***p < 0.001).

All compounds successfully reduced mxa, infΦ1, and ccl19 expression to the levels equivalent or even lower than the one of untreated controls. Baricitinib and Molnupiravir treatments also led to significantly reduced ccl20a.3 levels. The drug treatments also affected inflammasome pathway genes; ptg was upregulated to control levels by Remdesivir and Molnupiravir, while other drugs led to four-fold upregulation. Molnupiravir and Sabizabulin also successfully downregulated il4 expression. Of note, we also observed that in the absence of spike protein, drug treatments were already able to alter the expression of some of the immune response genes (Fig. 7). Of all five drugs, Remdesivir had the broadest effect on inflammatory genes in the absence of spike protein treatment.

Figure 7
figure 7

Remdesivir, Ritonavir, Molnupiravir, Baricitinib or Sabizabulin alter the inflammatory gene response in the absence of spike protein. qPCR results for inflammation associated gene expression in zebrafish larvae treaded with selected compounds. Remdesivir, Ritonavir, Molnupiravir, Baricitinib or Sabizabulin treatment is compared to larvae treated with DMSO. Bars represent mean and SD of 6 biological replicates. Each replicate consists of a pool of 5 larvae. Asterisks represent p-values from the Mann–Whitney U test * < 0.05; ** < 0.01;*** < 0.001. Note that in all graphs the same DMSO control samples were used.

Overall, these results indicate that the zebrafish model is useful not only for screening phenotypic and behavioral alterations mediated by drug treatment, but also to assess compound effectiveness in attenuating the immune response after SARS-CoV-2 spike protein exposure. It also highlights that at Cmax, clinically relevant compounds can lead to alterations of embryonic development in the zebrafish.

Data access via an online data app

In order to facilitate mining of this large dataset we created an online data app using the open-source Streamlit app framework (Fig. 8 and accessible via this link: https://share.streamlit.io/alernst/covasc_dataapp/main/CoVasc_DataApp.py. The online data application allows to access raw measurements as well as batch corrected data and provides access example images for each treatment according to one’s individual interests.

Figure 8
figure 8

Generation of a Web-based Data App. (a) The data app overview window shows which information can be found in the app. A heatmap with the previously described scores, results of the literature search, survival data, morphological analysis, activity analysis and images of the larvae. (b) Example tab opened, showing interactive violin plots for the larval length, which allows to see the distribution and individual data points. The checkboxes allow to view the corrected or raw data batch, and replicates can be grouped or shown individually. https://share.streamlit.io/alernst/covasc_dataapp/main/CoVasc_DataApp.py.

Here an overview of the functionality (Fig. 8a): on the left side, we allocated a side bar to select a specific tested compound or group of compounds. In the main window, tabs were created to allow visualization of an overview heat map, literature data, mortality of the treatments, morphological analysis and behavioral analysis. Additionally, we chose representative example images for each drug treatment group in each experimental replicate. An overview in brightfield and GFP fluorescence was provided in the last tab. Within the tabs we gave options to visualize the data (Fig. 8b). We provided options to show batch corrected data (check box: “Standardize to global median”) and to visualize individual experiments with the respective control or to group replicates by the treatment.

By providing the data app, we enable better insight into the data and facilitate visualization of all screening results beyond the ones shown in the figures. Furthermore, this app could be considered as a layout to be implemented for other screening projects.

Discussion

Our motivation to perform this study was to support drug repurposing efforts in the context of the COVID-19 therapy, however, given that FDA approved compounds have been tested, the assay provides information that can be useful beyond this disease. Despite of the rapid progress of drug development and many clinical trials, the quest on the search of an ideal drug therapy against COVID-19 continues. Even though several compounds have been approved for COVID-19 treatment, most of these drugs have not undergone thorough investigation for pregnant individuals52,53,54,55, or are not indicated during pregnancy without providing further information, as is the case for Remdesivir, Molnupiravir, and Ritonavir56 . Here we report the use of the zebrafish model to assess effectiveness and side effects on cardiovascular development of compounds repurposed for COVID-19 treatment. From the compounds that are currently used in the clinic Baricitinib, Remdesivir, Ritonavir had been studied for the induction of embryonic developmental defects in the zebrafish, but not in the context of cardiovascular development, and not all at concentrations corresponding to Cmax57,58. To our knowledge the effect of Molnupiravir and Sabizabulin on zebrafish development have not been studied to date.

We confirmed that Ivermectin severely affected embryonic development. Indeed, prevention and treatment of COVID-19 by Ivermectin has been warned against by the FDA due to lack of evidence for therapeutic benefits and adverse effects already during the pandemic (FDA, 2021). The compound Hydroxychloroquine caused mild defects on the cardiovascular system. While Hydroxychloroquine’s use in COVID19 treatment has been discontinued, it is still used for treatment of Malaria. Thus, further work elucidating possible side effects on embryonic development are still of relevance. Regarding compounds currently used for COVID-19 treatment, we observed about 20% mortality at Cmax for Sabizabulin. This microtubule assembly disruptor has entered in a III Clinical trial stage and was shown to reduce deaths of severe hospitalized patients by almost 25%. Our results call for further work to study contraindication of Sabizabulin during pregnancy. Molnupiravir and Paxlovid (Nirmatrelvir and Ritonavir)41 were approved by the FDA, but the approval is still based on little clinical experience in particular regarding patients with health conditions42. Up to this point, we studied only individual compounds and no combination treatments like Paxlovid. Comparing the results of two viral-replication inhibitors Remdesivir and Molnupiravir, both treatments seem to be inappropriate as severe adverse effects were observed. Ritonavir treatments caused overall little embryonic alterations, but led to significant impairment of the heart functionality and Baricitinib treatment also affected heart rate and embryonic behavior. Overall, all five compounds let to a larger or lesser extent to alterations in embryonic development.

We have also established a workflow to treat 5 dpf embryos with SARS-CoV-2 SPIKE protein to assess whether any of the relevant candidates could affect immune signaling induced by spike treatment. All five drugs tested for modulation of immune response after SARS-CoV-2 SPIKE treatment showed a significant reduction in inflammatory gene expression, namely mxa, infΦ1, and ccl19a.1. Our results show that the zebrafish model allows not only to study the immunogenic nature of Sars-Cov-2 viral particles, that can occur in the absence of infection (Tykalska et al. 2022) but also to address how compounds modulate such an inflammatory/immunogenic response. Of note, we cannot fully exclude that the compounds inhibit spike effect through external interaction within the immersion medium. We did observe some differences in the gene expression response to spike protein compared to (Tykalska et al. 2022) with no significant effect of spike treatment on il1b, tnfa, nfkb1, cxcl8a, ifnγ, caspa and il10 expression. This is likely due to differences in the embryonic stage used for the assay and due to the use of whole embryos in our study.

Even though several lines of evidence point to the validity of the zebrafish model to study physiological effect of chemical compounds as a first hint towards its effect in mammals, including humans6, species-specific effects such as differences in target conservation cannot be fully excluded. A further limitation is the route of administration. In comparison to treatments in mammals, zebrafish embryos are commonly treated by immersion. This treatment procedure may limit the estimation of the actual dose received by the organism due to limited penetration. Microinjection into the yolk of compounds with reduced water solubility might be needed to further dissect maximal effects on embryonic development. Compared to other zebrafish screenings59,60, the used standard concentration of 1 µM is relatively low and should thus cause more specific effects and less general toxicity. Additionally, less issues regarding solubility can be expected. However, the concentration of individual compounds might vary massively if compared to the concentrations used in clinical practice. Indeed, effects seen upon low drug dosage were also reported to have an effect in mammals, indicating translational potential61. Therefore, to provide more translationally meaningful results we included drug treatment on zebrafish using concentrations relevant to clinical studies (Cmax).

Automated screening setups produce large amounts of image data, which have to be analyzed in a standardized workflow. Machine learning has become an important tool in image analysis and is based on providing examples of a specific target structures to the computer62. In particular, convolutional neural networks nowadays outperform most conventional image analysis tools. In biomedical and medical image segmentation, U-Net architectures prove to be highly versatile63. Here we provide a deep learning model for fast and precise segmentation of embryonic zebrafish hearts and ISVs between 4–5 dpf. The model can be used directly for future screening assays, using the same imaging modalities or can be employed for transfer-learning to reduce the amount of necessary training data.

Besides the code for the U-Net training and analysis, we provide with this work a variety of tools to facilitate future zebrafish developmental screenings. We include multiple Jupyter notebooks and Python functions for image and data analysis, ImageJ-macros for image processing and semi-automatic analysis, the code to perform a systematic literature search, a customized template for a 3D-printed mold to position larvae in 96-well plates and a data app including the full code to allow individual data visualization.

While the herein tested compounds are overall well-studied drugs used in advanced clinical trials and already on the market, our screening assay can also be applied at earlier stages of drug development: before or in parallel to entering clinical trials. The presented screening pipeline requires nine days for a single researcher for 30–40 drugs to be tested at a given concentration. The steps included are: setup of the fish crosses until staging (three days), drug treatment (two days), imaging of morphology (one day) and behavior (one day) and analysis (two days). Future work can aim at reducing the time for analysis by fully implementing deep learning for all analyzed parameters. We avoided such an approach to correlate results from two sources, involving easy and fast counting in a semi-automatic approach and using automatic U-Net based segmentation for the more time-consuming parts. Overall, these two approaches combined with batch correction and scoring delivered highly coherent results pointing out relevant effects on embryonic development.

Many medicines are still insufficiently studied in the context of development, posing serious limitations in the treatment of pregnant women. We envisage that pipelines using the zebrafish model as the one presented in this study will be of interest to not only in the context of COVID-19 treatment but also for characterization of side effects on embryonic cardiovascular development of further treatments.

Material and methods

Zebrafish husbandry and studies on embryos

Experiments were performed with zebrafish (Danio rerio) embryos and larva at Institute of Anatomy (National License Number 35) and Institut für Infektionskrankheiten (National Licence Number 113) from the University of Bern. Adult fish needed for breeding were maintained at the Institute of Anatomy and raised and maintained at maximal 5 fish/L with the following environmental conditions: 27.5–28 °C, with 14 h of light and 10 h of dark, 650–700 μs/cm, pH 7.5 and 10% of water exchange daily. Animal studies were approved by the Animal Care and Experimentation Committee of the Canton of Bern, Switzerland (license BE27/2021). All experiments were performed in accordance with the guidelines and regulations approved by the Animal Care and Experimentation Committee of the Canton of Bern, Switzerland. The study was designed in accordance with ARRIVE guidelines. For larvae > 5 dpf used for qPCR analysis, no animals were excluded. Used larvae were not randomized. Confounders were minimized by using as controls larvae from the same clutch as those used for drug treatments. The experimenter was aware of the group allocation at all stages of the experiment. The number of treated and imaged embryos for each experiment and drug is documented in Supplementary Table 1.

Breeding and staging

Adult zebrafish (wild type or transgenic) were kept in family breeding tanks overnight, male and female fish separated by a transparent screen. For the morphology assay homozygous Tg(fli1a:GFP)y1Tg and Tg(myl7:mRFP)ko08Tg were crossed to obtain heterozygous double transgenic Tg(fli1a:GFP)y1Tg;(myl7:mRFP)ko08Tg. For behavioral studies, we used zebrafish from the wildtype AB strain. To initiate synchronous breeding the screens of all breeding tanks were removed at the same time and subsequently eggs were collected within 30 min and kept in E3 medium (5 mM NaCl, 0.17 mM KCl, 0.33 mM CaCl2, 0.33 mM MgSO4) with Methylene Blue (10–5%). Eggs from different breeding tanks were collected and then split into petri dishes at equal density. Confounders were minimized by taking controls from the same clutch as the experimental group. 30 min after egg collection unfertilized eggs and those that did not transition to two-cell stage, were removed. The next day, at 24 h after removing the screens, chorions were removed by incubating 2 mg/ml Pronase in E3 medium for about 3 min until gentle shaking of the petridish freed the larvae from the chorion.

Pharmacological treatment

Larvae were staged according to64 before starting the drug treatments. Then, 28 hpf embryos were transferred to 24 well plates, 10 embryos per well in 2 ml E3 medium. A new 24 well plate was prepared with 975 µl treatment solution (E3, 5 mM HEPES-Buffer with or without 0.003% 1-phenyl-2-thiourea, PTU, 0.25% DMSO) containing the dissolved drug. Embryos were transferred in 25 µl of medium, to reach 1 µM drug concentration. For one experiment, typically 10 drugs were tested and two wells were kept with control embryos (treatment solution as described above with only DMSO). A full list of experimental replicates can be found in Supplemental Table S1.

The well plate was covered from light and placed into the incubator at 28 °C. After 24 h of incubation the drug solution was replaced by fresh solution. After additional 24 h the embryos were washed three times with 2 ml E3 medium (with or without PTU) and kept in the respective medium until the effects are recorded.

For RNA extraction, groups 5 dpf AB wild type zebrafish were placed in 24 well plates. Each well contained 1 ml of fish water with 0.1–15.2 µM drug (fish water, 5 mM HEPES-Buffer, 0.25% DMSO) (refer to Fig. 1C) in the presence or absence of SARS-CoV-2 spike protein (5 ng/ml; GenScript). Controls were kept in the treatment solution with only DMSO. The plates were kept in the incubator at 28 °C for 48 h without medium change. Before the collection for RNA, zebrafish were euthanized with 0.16% Tricaine and washed with PBS. We have chosen to use 5 larvae per each treatment, as this enables to extract sufficient RNA for cDNA production (200–300 ng of RNA was obtained per larvae). We also used GPower analysis to predict how many biological replicates coming from different parents we would need for each treatment (n = 6). In total 360 larvae at 7 dpf were used for this experiment.

Preparation of 96 well plates with 3D printed mold

For lateral positioning of larval zebrafish in 96 well plates (Greiner) we used either a 3D printed mold provided by Acquifer Germany or a custom-made 3D print (Template: https://github.com/Alernst/CoVasc_DrugScreen/tree/main/3D_print). Using a multipipette, 70 µl of melted 1.5% low melting agarose (Promega) was injected into each well. Eventually occurring bubbles were removed using a metal needle. The template mold was carefully inserted in the plate and the plate with the mold was placed at 4 °C for 30 min. Subsequently, the mold was pulled out and the plate with agarose beds was ready to use.

Organismal and cardiovascular development assay

Up to four experiments were run per day, for this reason larvae were taken between 96 and 104 hpf. Each experiment was compared to its own control group with the same embryonic stage. Embryos that were harmed due to handling were excluded from the assay. Between 7 and 10 embryos were taken per experimental replicate per drug. Embryos were transferred from the 24-well plate into a 12 well plate for easier handling with 4 ml volume per well. For anesthesia, tricaine was added at 0.16 mg/ml (pH 7) and the larvae were immediately transferred to the prepared 96 well plate with agarose beds using a cut yellow pipette tip with 70 µl volume, resulting in a final concentration of 0.08 mg/ml tricaine.

The imaging of the larvae in the 96-well plate was performed using an automated microscopy platform (Acquifer Imaging Machine, Acquifer Germany). Initially, an overview image of each entire well was taken with the 2 × objective (pre-scan). From all the brightfield and green fluorescence channel overview images (11 Z-planes), we chose a reference image of the trunk and the head. This image served as a template. Subsequently, an implemented algorithm (template matching) was used to detect trunk and head region of the rest of larvae. After successful detection, a 4 × objective was used to acquire Z-stacks (25 Z-planes) of all trunks in brightfield and green fluorescence channel. Next, the 10 × objective was used to acquire a time series of 100 frames of the head region focused on the heart. For initial screening experiments, 300 frames in brightfield and red fluorescence channel were acquired.

Semi-automated image analysis

The analysis was performed using the Fiji-software. A customized ImageJ-macro (https://github.com/Alernst/CoVasc_DrugScreen/blob/main/macros_imagej/AIM_Semiauto_Analysis.ijm) was written to sequentially open the images of the dataset without showing the treatment group. The macro guides the user through the workflow. At each step requiring manual interaction a dialog window shows the instructions which action to perform. (1) To count ISVs draw a line crossing all ISVs in the green fluorescence channel. The macro creates a kymograph along the Z-axis and detect the local maxima in the kymograph. The user is asked to supervise the detected maxima. (2) Based on the shape of the pericardial cavity, a decision is taken whether cardiac pericardial effusion is present or not. (3) Draw a straight line along the anterior–posterior axis of the larva from head to tail to measure the linear axis length, following recommendations65. (4) The macro opens the heart time series, performs a maximum intensity projection along the time axis, detects the heart in the T-projection by thresholding. The macro suggests a region of interest (ROI) that contains the heart, if this ROI is misplaced, the user can correct the location. The macro now creates a perpendicular line in the center of the lower edge of the ROI, creates a kymograph and detects the local maxima, which needs supervision by the user. All measurements were summarized in a result table containing embryo length, heart rate, ISV number and presence of pericardial effusion. Only after the semi-automatic steps, the drugs were assigned to the analyzed data to ensure an unbiased analysis.

Automated image analysis

The second part of the analysis was performed using Python via Jupyter notebooks. The Jupyter notebook contains each step of the analysis and can be found in the repository, including all required packages (https://github.com/Alernst/CoVasc_DrugScreen/tree/main/Deep_Learning). The aim was to perform a segmentation of the ISVs, similar to66, and the heart distinguishing between ventricle and atrium, comparable to67.

The core components for segmentation are two newly trained convolutional neural network with a U-Net architecture. The first model, now called IsvSeg, was trained to perform a binary segmentation to distinguish background from ISVs, 30 maximum intensity projections (2048 × 2048) were manually anotated and a data augmentation was performed, also out-of-focus data were included in the training. The second model for multiclass annotation, named HeartSeg, was trained to segment the ventricle, atrium and background. For training, 348 images (512 × 512) were manually labelled and out-of-focus data were again also included in the training. The model was tested on 96 new set of images and was evaluated with an unweighted dice coefficient of 0.902.

The HeartSeg model was applied to 100 frames of each heart time series in the red fluorescence channel. To increase the quality of the data, small objects under a threshold area were removed and if any frame did not contain a segmentation of the ventricle, the entire dataset was removed. To calculate the ejection fraction the ventricle area at each time point of the times series was extracted from the segmentation. A peak detection algorithm was used to identify all systoles and diastoles. Due to the shape of the embryonic fish heart the 2D area can be used to calculate the ejection fraction68. A median systole and diastole area was calculated per embryo, which was used to obtain the ejection fraction for each larva, using the formula below:

$$EF \left(\%\right)=\frac{(EDV-ESV)}{EDV}\times 100=\frac{SV}{EDV}\times 100$$

Alterations from the standard experimental procedure in a subset of experiments are explained in the comment section of Supplementary Table 1.

The IsvSeg model was applied to maximum intensity projections of all 4 × images of the green fluorescence channel. We measured each ISV’s major (length) and minor axis (width), using the functions included in SciKit image, and calculated a median for each larva. The median ISV area was estimated by multiplying median width and length.

As a final step, we combined the information related to ISV area, ejection fraction, larva length, heart rate, ISV number and pericardial effusion for further meta-analysis (Fig. 2c).

Behavioral analysis

To measure the activity and potential impact of drugs on behavioral patterns of larval zebrafish we used the DanioVision™ recording chamber (Noldus Inc., Wageningen, Netherlands). Between 115–120 hpf, the previously treated wild-type larvae were placed individually in the wells of a 96 well plate with 200 µl E3 medium pre-warmed to 28 °C. No anesthesia or fixation was applied to allow free swimming in the well. The plate was mounted in the DanioVision™ recording chamber. A programmed light cycle is automatically switching on and off the light inside of the recording chamber. Then, larvae had 30 min of accommodation time in the dark, followed by 6 cycles of 10 min dark and 10 min bright periods. During these 2.5 h in total, the larvae were recorded. Larvae tracking upon light stimulation was carried out using the EthoVision XT™ software (Fig. S1). The screening analysis was performed using the 1 min binned track results calculated by the software as xlsx-file. Prior to analysis, incorrectly tracked larvae were excluded by visual inspection of the movies. A table was generated to annotate bright and dark phases as applied in the experiment. All calculations with detailed description were made in Jupyter notebooks which can be accessed on (https://github.com/Alernst/CoVasc_DrugScreen/blob/main/CoVasc_Data_Update.ipynb). From the annotated data a median was calculated for the swimming velocity, moving duration and the distance moved in dark and bright. The logarithmic bright/dark ratio was calculated by dividing velocity in the bright by velocity in the dark and calculating the natural logarithm. To handle zero values, 0.001 was added to both values.

qPCR

RNA extraction of 7 dpf larvae (5 dpf + 2 days of drug treatment) in batches of 5 individuals was performed using TriZol reagent (Sigma-Aldrich) and RNA purification kit (Zymo). cDNA was generated using Maxima First Strand cDNA synthesis kit (Thermo Fisher Scientific). RT-qPCR analysis was performed using PowerUp SYBR Green Master Mix (Thermo Fisher Scientific) and following primers:

Primer name

Sequence

References

mxa_F

TTGACCTCCCTGGCATTGCA

49

mxa_R

GATTGTCTCTTGCCTTGTAACA

il4/13b_F

GCAGGAATGGCTTTGAAGGGTAAA

il4/13b_R

AAACTCCTTCATTGTGCATTCCCC

ifnΦ1_F

CGCAAAGCCAGCACACAAGGA

ifnΦ1_R

CTCCGGATCTGCTCCCATGCT

il1b_F

CGCTCCACATCTCGTACTCA

il1b_R

ATACGCGGTGCTGATAAACC

ccl20a.3_F

TGATGGTGCTGACAATCGTG

ccl20a.3_R

CTTTGGACGGGTCTGTGCA

tnfa_F

GCGCTTTTCTGAATCCTACG

tnfa_R

TGCCCAGTCTGTCTCCTTCT

ccl19a.1_F

GCCCACGTGATGCTGTAATA

ccl19a.1_R

ACAGCGTCTCTCGATGAACC

nfkb1_F

TTCTTCTTGGTCACGTGCAG

50

nfkb1_R

ACTCTCAGCATCCGCATCTT

cxcl8a_F

GTCGCTGCATTGAAACAGAA

cxcl8a_R

CTTAACCCATGGAGCAGAGG

il10_F

AACTCAAGCGGGATATGGTG

il10_R

ATCAAGCTCCCCCATAGCTT

infg_F

CTTCAGACAACCAGCGCATA

infg_R

TTTTCCAACCCAATCCTTTG

ptgs2a_F

TGGATCTTTCCTGGTGAAGG

ptgs2a_R

GAAGCTCAGGGGTAGTGCAG

caspa_F

CGACGTCAGGGAGATAAGGC

caspa_R

TGGATACTAAGGTTTTGAACGACG

pycard_F

ATTTTGAGGGCGATCAAGTG

pycard_R

GCATCCTCAAGGTCATCCAT

rps11_F

GATGGCGGACACTCAGAAC

 

rps11_R

CCAATCCAACGTTTCTGTGA

 

qPCR was performed on an Applied Biosystems QuantStudio 6 Pro System. 6 replicates, each comprising a pool of 5 embryos were performed.

Statistics and batch correction

All results from Fiji and Python were saved as csv-files and subsequently merged to a single data table. The Mann–Whitney U test was applied to compare each treatment group to the control group of the same experiment (acquired on the same 96 well plate). Additionally, we calculated the effect size and the relative difference of the medians compared to the control. To evaluate multiple experiments and replicates together a batch correction was performed. This correction was achieved by calculating the median of all control values (global median). The difference between the global median and the experimental control median was calculated and each measurement of the experiment was increased or reduced by this difference. The results from qPCR for the heatmap and histograms were analyzed using The Mann–Whitney U test.

Effect score

The score to estimate relevant effects was calculated by multiplying the –log10(p-value) of the Mann–Whitney U test with the hedges effect size and the relative difference of the median. The Mann–Whitney U test was used to estimate, whether the difference between the control and the treatment is significant. The hedges effect size was utilized to obtain a quantitative measure of how strong an effect is considering the spreading of the data points. The relative difference of the median was used to see how different the treatment median is from its respective control.

Percentage effect

The Mann–Whitney U test estimates whether the difference between the control and the treatment is significant. Following this, the fold change of the treatment and control is calculated. From this value, percentage effect is calculated. This helps in gaining an estimation of how much treatment varies from the control and there by understand the severity of the effect.

Literature search

To focus on compounds with poor characterization of possible side effects during embryogenesis in the context of COVID-19 treatment, we performed a systematic literature screen (Fig S3a,b). We retrieved the PubMed ID of each article and counted the total number of articles found. We listed the 15 compounds with most peer-reviewed publications related to COVID-19 and analyzed to which extent these had also already previously been studied in the context of cardiovascular development and embryogenesis (Fig. S3c). We found extensive number of research articles on half of the compounds, and 73.8% of tested drugs had more than ten articles in the cardiovascular system and almost 50% had more than ten articles in embryogenesis (Fig. S3d). However, we found that for approximately half of the compounds studied in the context of COVID-19 there was no literature related to their effect on the cardiovascular system.

The systematic Pubmed search for all investigated compounds was performed using a Jupyter notebook (https://github.com/Alernst/CoVasc_DrugScreen/blob/main/CoVasc_Data_Update.ipynb). The PubchemPy package was used to search synonyms for each compound. These synonyms were combined with specific keywords by “AND” or “OR” operators. Example search term: “Covid-19” AND “Remdesivir” OR “Covid-19” AND “Veklury”. These search terms were generated for each compound combined with context keywords. These search terms were submitted to search in Pubmed using the EntrezPy package. The results were collected as counts and Pubmed IDs. The Pubmed IDs of different context searches were compared for overlapping articles in different contexts.

Online app

The app was generated in Python code. The code is accessible as a GitHub repository (https://github.com/Alernst/CoVasc_DataApp). The app is hosted via the Streamlit Cloud.