The AB strain of wild-type zebrafish was purchased from the Institute of Hydrobiology, Chinese Academy of Sciences, and raised in a standardized zebrafish breeding facility (Beijing Aisheng Technology Development Co., Ltd.) at Hunan Normal University. The animal experimental protocol was approved by the Institutional Ethics Committee of the Guangdong Academy of Medical Sciences (KY2024-847-01) and was performed in accordance with the relevant guidelines and regulations.

The CRISPR/Cas9 gene-editing system was used to generate zebrafish bmpr2a and bmpr2b knockout. Exons 8 and 9 of the bmpr2a and bmpr2b genes were selected as potential target sites, and the website http://crispor.tefor.net/crispor.py was utilized to design Guide RNA. The Guide RNA sequence was then linked to the pUC57 sgRNA backbone plasmid through homologous recombination, yielding pUC57-bmpr2a-sgRNA1 and pUC57-bmpr2a-sgRNA2, along with pUC57-bmpr2b-sgRNA1 and pUC57-bmpr2b-sgRNA2, respectively. The Guide RNAs are listed in Supplementary Table S1. Guide RNA was amplified in vitro and subjected to in vitro transcription experiments (Riboprobe® System-T7 Translation Kit (Promega, P1440) to produce sgRNA. The sgRNA (20 ng/μL) was mixed with Cas9 protein (TrueCut Cas9 v2, Thermo Fisher Scientific, A36499, 300 ng/μL) and injected into zebrafish at the one-cell phase. The positively knocked-out zebrafish were screened in the F0 generation and sequenced to verify the knockout. F1 were obtained from F0 zebrafish mated with wild-type (WT) zebrafish, and F1 were partially sequenced to verify the knockout strain. The sequence and genotype primers are listed in Supplementary Table S2. bmpr2a and bmpr2b double knockout zebrafish were hybridized using bmpr2a and bmpr2b knockouts and identified simultaneously using bmpr2a and bmpr2b primers. In this study, we used four main genotypes, namely, the WT (bmpr2a ^+/+ ;bmpr2b ^+/+ ), bmpr2a homozygotes (bmpr2a ^−/− ;bmpr2b ^+/+ ), bmpr2b homozygotes (bmpr2a ^+/+ ;bmpr2b ^−/− ), and double homozygotes (bmpr2a ^−/− ;bmpr2b ^−/− ), which were obtained from the self-crossing of double heterozygotes (bmpr2a ^+/−;bmpr2b ^+/− ).

Quantitative real-time polymerase chain reaction (qRT-PCR) was performed as previously described (Shi et al., 2019). For qRT-PCR analysis of whole embryos, genotyping was conducted using tail biopsies from 48 hours post-fertilization (hpf) zebrafish embryos, while the remaining embryonic tissues were immediately stored at −80 °C for subsequent RNA extraction (each group contained six zebrafish embryos). For qRT-PCR analysis of the heart tissues, intact heart tissues (including the outflow tract and inflow tract) were microdissected from 48 hpf zebrafish embryos. The residual tissues were retained for genotyping, and a minimum of 20 heart tissues were collected from each experimental group.

Following genotyping, samples from the same genotype were grouped, and total RNA was extracted using TRIzol (Invitrogen). The cDNA library was then synthesized according to the manufacturer’s instructions (TransScript® One-Step gDNA Removal and cDNA Synthesis SuperMix, AT311-03). Finally, qRT-PCR was performed under standard PCR conditions using SYBR Green PCR Master Mix (TaKaRa). All gene expression levels were standardized to GAPDH expression and analyzed using the 2^−ΔΔCT Livak method. Each genotype was represented by at least three independent biological replicates. All qRT-PCR primers are listed in Supplementary Table S2.

The samples were prepared according to the qRT-PCR analysis method for whole embryos. After sample preparation, they were transported to Majorbio for total RNA extraction, cDNA library construction, and RNA-seq analysis. All data analyses were performed on the Majorbio Cloud platform (www.Majorbio.com). Differentially expressed genes (DEGs) were identified using screening criteria with a significance threshold of |log₂FC| >2.0 and p < 0.05. The original data were submitted to the NCBI Sequence Read Archive (SRA) database (SRA number: PRJNA1336002).

Reverse transcription-PCR was used to amplify the mRNA sequence of the gene for probe preparation, and the reverse primers were added to the T7 promoter sequences. Digoxigenin-labeled antisense RNA probes were synthesized through in vitro transcription using the Riboprobe® System-T7 Transcription Kit (P1440, Promega) and ROCHE DIG RNA Labeling Mix (REF 11277073910, Roche), according to the manufacturers’ instructions. The zebrafish embryos were fixed in 4% paraformaldehyde, treated with a gradient methanol series (25%, 50%, 75%, 85%, 95%, and 100%), and stored in 100% methanol.

Whole-embryo in situ hybridization (WISH) was performed as previously described (Oxtoby and Jowett, 1993). In brief, stored embryos (30–50 embryos/tube) were rehydrated in a graded methanol/PBST series, digested in 10 mg/mL proteinase K (PBST), fixed again in 4% paraformaldehyde, and pre-hybridized in hybridization buffer. Subsequently, the pre-hybridization buffer was replaced with fresh hybridization buffer containing digoxigenin-labeled RNA probe (300 ng) and incubated overnight at 65 °C. Furthermore, the embryos were washed in different buffers, blocked with 2 mg/mL BSA and 2% sheep serum, and incubated with pre-adsorbed antibody overnight. The embryos were then stained using MABT and AP substrate chromogenic solution until an optimal signal was obtained (approximately 10 min–50 min). The staining reaction was terminated by washing in several changes of PBS and PBST. Finally, embryos were kept at 4 °C and photographed in 6% methylcellulose using a Leica TL 5000 microscope (Leica, Germany). After the embryos were photographed, they were collected for genotyping.

Cardiac function analysis in zebrafish embryos was performed as previously described (Fink et al., 2009). In brief, M-mode was conducted using a high-speed EMCCD camera to capture 10-s movies of zebrafish heart activity at 48 hpf under a ×20 microscope objective. The recorded cardiac motion videos were analyzed using custom heart analysis software (SOHA software) to derive the functional parameters: heart rate (HR), heart period (HP), diastolic interval (DI), systolic interval (SI), diastolic diameter (DD), systolic diameter (SD), and fractional shortening (FS). All data were visualized as scatter point histograms. Following imaging, embryos were collected for genotyping.

Embryos were incubated at 28.5 °C in Petri dishes containing fish water. To prevent the formation of melanin pigments, PTU (Sigma) was added to fish water at a final concentration of 0.003% at the end of gastrulation. Then, the zebrafish embryos at 48, 72, and 96 hpf were immobilized using 6% methylcellulose and positioned with the abdomen facing upward, and the pericardial cavity and cardiac phenotypes were imaged using an Axiocam (Zeiss) for subsequent analysis.

Intact heart tissues (including the outflow tract and inflow tracts) were microdissected from 48 hpf zebrafish, and the residual body tissues were retained for genotyping. Following genotyping, heart tissues from the same genotype were pooled, with a minimum of 30 hearts per experimental group, for subsequent protein extraction. Each genotype was represented by at least three independent biological replicates.

Total protein was extracted by homogenizing pooled heart tissues in 100 μL radioimmunoprecipitation assay (RIPA) buffer (Beyotime). Proteins were then separated by sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) using Future PAGE^TM 4%–20% 15-well gels (ACE). Separated proteins were transferred to polyvinylidene fluoride (PVDF) membranes (Millipore), which were subsequently blocked with 5% skim milk (CST) in Tris-buffered saline with Tween-20 (TBST) for 2 h at room temperature. After blocking, membranes were incubated overnight at 4 °C with the following primary antibodies: GAPDH (1:50,000, 60004, Proteintech), Vim (1:1,000, T55134, Abmart), Myl7 (1:5,000, GTX128346, GeneTex), Cdh2 (1:1,000, CST, 13116), Cdh1 (1:1,000, 3195, CST), and Fn1 (1:1,000, ab268020, Abcam). The proteins were then incubated with horseradish peroxidase (HRP)-conjugated anti-rabbit immunoglobulin G (IgG) or anti-mouse IgG (1:10,000; Abmart) for 2 h at room temperature. After washing with TBST, protein bands were visualized using an Immobilon@ Western Chemiluminescent HRP Substrate kit (Millipore). The relative signal densities of the protein bands were quantified using ImageJ (1.51 version, NIH) and normalized to GAPDH to account for variations in protein loading.

Immunostaining was performed as previously described (Yang and Xu, 2012). The tail biopsies from 48 hpf zebrafish embryos were used for genotyping, while the remaining embryonic tissues were fixed in 4% paraformaldehyde (PFA) at 4 °C and embedded in paraffin. Paraffin-embedded tissues were then sectioned into 2 μm–4 μm cross-sections, which were subjected to immunostaining with a primary antibody against myosin light chain 7 (Myl7, 1:200, GTX128346, GeneTex), goat anti-rabbit IgG (H + L) cross-adsorbed secondary antibody (Alexa Fluor™ 594, A-11012, 1:1000, Thermo Fisher Scientific), and 4′,6-diamidino-2-phenylindole (DAPI, 1:5000, 28718–90-3, Proteintech) for nuclear counterstaining.

The heart sarcomere structure was imaged using a Nikon confocal microscope (AX-NIS-Elements). The width and length of the sarcomere bands were measured based on fluorescence images. Three zebrafish embryos with well-preserved and clearly stained heart tissues were selected for analysis for each experimental group. For each heart sample, the width or length of the Z-disc was measured in at least three distinct myofibrils to ensure statistical reliability.

The Tg (flia:GFP) line (Shi et al., 2020) was crossed with bmpr2a ^+/−;bmpr2b ^+/− to generate double heterozygotes, producing flia:GFP;bmpr2a^+/−;bmpr2b^+/−  offspring, which were then subjected to self-crossing. Zebrafish embryos were maintained at 28 °C until 144 hpf, and fluorescent images were captured and analyzed using a Nikon confocal microscope (AX-NIS-Elements). After obtaining the images, embryonic tissues were collected for genotype identification, and the same genomes were grouped. The left and right valves were delineated with yellow and red dotted lines, respectively, and their diameters were measured using Digimizer software. At least six samples were used per group.

For comparisons among more than two groups, the Shapiro–Wilk test was utilized to determine the normality of data distribution, and the Levene’s test was implemented to examine the homogeneity of variances. When parametric assumptions were met, one-way analysis of variance (ANOVA) was performed, followed by Tukey’s post hoc test for multiple pairwise comparisons. When parametric assumptions were violated, the non-parametric Kruskal–Wallis test was implemented, followed by Dunn’s post hoc test combined with the Benjamini–Hochberg correction to account for multiple comparisons biases. All statistical analyses were performed using biological replicates. Data are presented as the mean ± standard deviation (SD) and were visualized using GraphPad Prism software. Statistical significance was defined as follows: ns, p > 0.05; *, p < 0.05; **, p < 0.01; and ***, p < 0.001.

To characterize the function of BMPR2 in zebrafish, we used CRISPR/Cas9 genome technology to generate bmpr2a and bmpr2b knockout strains. Bioinformatics analysis using the NCBI and Ensemble databases revealed that both bmpr2a (NM_001039817.1, ENSDART00000056764.5) and bmpr2b (NM_001039807.1, ENSDART00000125961.3) consist of 13 exons. Guide RNAs were designed to target exons 8 and 9 of bmpr2a and bmpr2b, respectively, for gene disruption (Figures 1A,B). PCR screening and Sanger sequencing of the F0 generation offspring revealed a 202 bp deletion in bmpr2a, including 131 bp within exon 8 and 71 bp in intron 8 (Figure 1C). For bmpr2b, mutagenesis resulted in a 171 bp deletion with a 3 bp insertion and a 146 bp deletion in exon 8 and a 25 bp deletion in exon 9 (Figure 1D). Both mutant zebrafish lines exhibited frameshift alterations that were predicted to result in premature translational termination (Figure 1E). Consequently, these mutants lost most of the functional domains of the protein kinase, along with the amino acid residues D482 and D485, which directly interact with ACVRL1, a ligand of the BMP signaling pathway (Iwasa et al., 2023) (Figure 1E).

F1 heterozygotes (bmpr2a ^+/− and bmpr2b ^+/− ) were generated by crossing F0 founders with wild-type zebrafish, and F2 mutant offspring were obtained by self-crossing of F1 heterozygotes. The results showed that, in both bmpr2a and bmpr2b mutant zebrafish strains, homozygous offspring from heterozygote self-crosses did not conform to Mendelian inheritance (1:2:1), indicating that the homozygous offspring were partially developmentally lethal (Supplementary Figures S1A, C). The sex ratio of homozygotes of both bmpr2a and bmpr2b showed severe imbalance; 2 of 57 bmpr2a ^−/− individuals were female in bmpr2a ^−/− , and all 71 adult bmpr2b ^−/− individuals were male (Supplementary Figures S1B, D). Zhang et al. (2020) reported that the loss of bmpr2a impairs the formation of mature follicles in female zebrafish and spermatogonia meiosis in male zebrafish, whereas bmpr2b deficiency disrupts folliculogenesis (resulting in infertility in mutant females) in female zebrafish but does not affect the formation of mature spermatogonia in male zebrafish (Zhang et al., 2020). Notably, their study described gonadal hypertrophy and dysfunction in bmpr2a/b mutant zebrafish without significant shifts in the sex ratio. However, we observed marked sex-ratio imbalances in both bmpr2a and bmpr2b single mutants. This discrepancy may be attributed to differences in the sgRNA-targeting regions. In the study of Zhang et al. (2020), sgRNAs targeting bmpr2a and bmpr2b were designed for exons 2 and 1, respectively, both of which are located near the ATG start codon. In contrast, in our study, the sgRNAs for bmpr2a and bmpr2b target exons 8 and 9, respectively, with both regions located far from the ATG start codon. However, the molecular mechanisms underlying these phenotypic differences require further investigation.

Owing to the pronounced sex-ratio imbalance in bmpr2a ^−/− and bmpr2b ^−/− , female bmpr2a ^−/− zebrafish were crossed with male bmpr2b ^−/− zebrafish to generate double heterozygotes (bmpr2a ^+/−;bmpr2b ^+/− ). Self-crossing of these double heterozygotes and genotyping analyses showed that while all other genotypes adhered to Mendelian inheritance, double homozygotes (bmpr2a ^−/− ;bmpr2b ^−/− ) were absent in adult populations (0.00%, Supplementary Figure S1F). At 48 hpf, however, bmpr2a ^−/− ;bmpr2b ^−/− individuals comprised 7.50% of the offspring, approaching the expected Mendelian ratio of 6.25% (1/16; Supplementary Figure S1E). These results confirmed the developmental lethality of the double mutants, which is consistent with the findings of Zhang et al. (2020). In their study, a significant increase in mortality was observed in bmpr2a ^−/− ;bmpr2b ^−/− mutants at 30 days post-fertilization (dpf).

Subsequent functional analysis utilized offspring from bmpr2a ^+/−;bmpr2b ^+/− self-crossing, focusing on four genotype groups, namely, WT (bmpr2a ^+/+ ;bmpr2b ^+/+ ), bmpr2a homozygotes (bmpr2a ^−/− ;bmpr2b ^+/+ ), bmpr2b homozygotes (bmpr2a ^+/+ ;bmpr2b ^−/− ), and double homozygotes (bmpr2a ^−/− ;bmpr2b ^−/− ), for phenotypic characterization.

qRT-PCR was performed to assess bmpr2a and bmpr2b. Compared to that of bmpr2a ^+/+ ;bmpr2b ^+/+ , bmpr2a transcript levels were reduced by ∼95% in both bmpr2a ^−/− ;bmpr2b ^+/+ and bmpr2a ^−/− ;bmpr2b ^−/− (Figure 1F). Notably, bmpr2a expression was also reduced by ∼55% in bmpr2a ^+/+ ;bmpr2b ^−/− , indicating cross-regulation between the two paralogs. Conversely, bmpr2b expression was reduced by ∼90% in both bmpr2a ^+/+ ;bmpr2b ^−/− and bmpr2a ^−/− ;bmpr2b ^−/− and by ∼65% in bmpr2a ^−/− ;bmpr2b ^+/+ (Figure 1G). These results confirmed the successful generation of double-knockout zebrafish and revealed the mutual regulatory interactions between bmpr2a and bmpr2b. This regulatory crosstalk is consistent with prior observations that individual knockdown of bmpr2a or bmpr2b induces comparable heart laterality defects, whereas simultaneous knockdown does not increase these phenotypes (Monteiro et al., 2008).

Through gene structure analysis and by demonstrating consistency with previous findings, we successfully generated three mutant zebrafish lines, including two single-knockout strains (bmpr2a ^−/− and bmpr2b ^−/− ) and a bmpr2a ^−/− ;bmpr2b ^−/− double-knockout strain. Phenotypic analysis revealed that homozygous mutants for either bmpr2a or bmpr2b exhibited partial developmental lethality, while all bmpr2a ^−/− ;bmpr2b ^−/− double-homozygous embryos failed to survive to adulthood, highlighting the indispensable and synergistic roles of these genes in zebrafish development.

Given that the loss of bmpr2b leads to embryonic development lethality in zebrafish, we investigated whether bmpr2a and bmpr2b affect heart development and cardiac function. Morphological analysis of cardiac development revealed that cardiac looping abnormalities and pericardial edema occurred not only in the double homozygotes bmpr2a ^−/− ;bmpr2b ^−/− but also in single homozygous mutants (bmpr2a ^−/− ;bmpr2b ^+/+ and bmpr2a ^+/+ ;bmpr2b ^−/− ) at comparable frequencies (Supplementary Figure S2). The phenotype of heart looping abnormalities is consistent with that of prior studies, in which bmpr2a/b knockdown affected the establishment of left–right asymmetry in zebrafish (Monteiro et al., 2008).

Furthermore, using M-mode (Fink et al., 2009), we analyzed the key parameters of ventricular function at 48 hpf. The results showed that compared to WT controls, bmpr2a ^+/+;bmpr2b ^+/+, bmpr2a ^−/−;bmpr2b ^+/+, bmpr2a ^+/+;bmpr2b ^−/−, and bmpr2a ^−/−;bmpr2b ^−/− mutant zebrafish exhibited no significant differences in ventricular diastolic interval and systolic diameter (Figures 2A,B,G). However, mutant zebrafish showed prolonged ventricular systolic intervals, altered heart-rate ratios, prolonged cardiac cycles, and reduced ventricular diastolic diameters and fractional shortening (FS) (Figures 2C–F,H,I), indicating impaired contractility and cardiac dysfunction. Collectively, these data indicate that bmpr2 signaling is critical for cardiac contractile and pacing functions, with loss-of-function mutant zebrafish developing heart failure at 48 hpf. Furthermore, altered heart-rate ratios, prolonged cardiac cycles, increased ventricular systolic diameters, and decreased ventricular diastolic diameters and FS were also observed in mutant zebrafish at 72 hpf (Figures 2J–N). This is consistent with previous research showing that BMP signals, such as BMP2 and BMP4, can regulate the differentiation and formation of sinusoidal node cells, thereby affecting cardiac pacing and the heart rate (Liang et al., 2021; Linscheid et al., 2019; Wang et al., 2023). An irregular heart rate impairs cardiac function through intracellular production of reactive oxygen species (Bergau et al., 2022) and altered expression of sarcomere structure genes and heart failure markers (Lee and Cha, 2021; Sossalla and Vollmann, 2018). Therefore, this may be one of the reasons for diastolic and systolic dysfunction and reduced FS in the three mutant zebrafish strains (bmpr2a ^−/− ;bmpr2b ^+/+ , bmpr2a ^+/+ ;bmpr2b ^−/− , and bmpr2a ^−/− ;bmpr2b ^−/− ).

The sarcomere structure is crucial for maintaining the cardiac architecture and enabling myocardial contraction (Zhang et al., 2023). Changes in the sarcomere structure can alter cardiac function (Crocini and Gotthardt, 2021). A previous study demonstrated that the thick filament network of the sarcomere can be visualized via myosin immunostaining (Yang and Xu, 2012). Consistent with this approach, we utilized Myl7 immunostaining to examine the sarcomere structure at 48 hpf in our experimental models. As shown in Figure 3A, sarcomeres in the WT group, bmpr2a ^+/+ ;bmpr2b ^+/+ , exhibited an ordered arrangement, and the same structural order was observed in the three mutant groups, namely, bmpr2a ^−/− ;bmpr2b ^+/+ , bmpr2a ^+/+ ;bmpr2b ^−/− , and bmpr2a ^−/− ;bmpr2b ^−/− (Figure 3A). Furthermore, comparative analysis revealed that compared with that of the bmpr2a ^+/+ ;bmpr2b ^+/+ group, the width of the bmpr2a ^−/− ;bmp2b ^+/+ , bmpr2a ^+/+ ;bmpr2b ^−/− , and bmpr2a ^−/− ;bmpr2b ^−/− embryos was significantly decreased, while the length of the bands was not affected (Figures 3B,C). Although the width of the sarcomere was affected in the three mutant strains, the changes were not sufficient to affect the contraction and relaxation functions of the nervous system; instead, together with the heart ratio, they induced diastolic and systolic dysfunction and reduced FS.

In zebrafish, the cardiac valve matures into a functional structure at 144 hpf, and the Tg (flia:EGFP) mutant line, which specifically labels endothelial cells, has been validated as a reliable tool for visualizing the valve architecture (Duchemin et al., 2019). To assess valve structural integrity, we analyzed the Tg (flia:EGFP) line at 144 hpf. Compared to that in WT (bmpr2a ^+/+ ;bmpr2b ^+/+ ) zebrafish embryos, the atrioventricular valve of bmpr2a^−/− ;bmpr2b ^+/+ , bmpr2a ^+/+ ;bmpr2b ^−/− , and bmpr2a ^−/− ;bmpr2b ^−/− zebrafish embryos exhibited a significantly thickened morphology (Figures 4A–C).

Cardiac valve development in zebrafish is initiated at 52 hpf and is governed by key markers, including nfatc1, fn1b, has2, and klf2a (Huang et al., 2018; Steed et al., 2016). To explore the molecular basis of valve developmental defects in bmpr2-knockout mutant zebrafish, we performed WISH to assess gene expression at 52 hpf. Compared to that in bmpr2a ^+/+ ;bmpr2b ^+/+ embryos, nfatc1, fn1b, and has2 showed significantly upregulated expression in bmpr2a^−/− ;bmpr2b ^+/+ , bmpr2a ^+/+ ;bmpr2b ^−/− , and bmpr2a ^−/− ;bmpr2b ^−/− mutant strains (Figure 4D). Conversely, klf2a expression was downregulated in all mutant strains (Figure 4D). qRT-PCR validated these findings, demonstrating consistent upregulation of nfatc1, fn1b, and has2 and downregulation of klf2a in all mutant strain embryos (Figures 4E–H). Collectively, these results indicate that bmpr2 loss-of-function mutation disrupts the transcriptional program governing valve development in zebrafish.

These results indicated that bmpr2-knockout mutant zebrafish exhibited abnormal cardiac function and valve development at 48 and 52 hpf, respectively. To uncover the underlying molecular mechanisms, we performed transcriptome sequencing (RNA-seq) on embryos at 48 hpf. DEGs were identified using strict criteria (p < 0.05 and FC ≥ 2 or ≤0.5). Compared to the bmpr2a ^+/+ ;bmpr2b ^+/+ group, the bmpr2a ^−/− ;bmpr2b ^+/+ group showed 753 upregulated and 126 downregulated genes, the bmpr2a ^+/+ ;bmpr2b ^−/− group had 1,197 upregulated and 234 downregulated genes, and the bmpr2a ^−/− ;bmpr2b ^−/− group contained 435 upregulated and 153 downregulated genes (Supplementary Figures S3, 4). In contrast, comparisons between double-knockout strains, bmpr2a ^−/− ;bmpr2b ^−/− , and single-knockout strains revealed minimal DEGs: bmpr2a ^−/− ;bmpr2b ^+/+ had nine upregulated and four downregulated genes, whereas bmpr2a ^+/+ ;bmpr2b ^−/− had one upregulated and 0 downregulated genes (Supplementary Figures S3, 4). Only 21 upregulated and 17 downregulated genes distinguished the two single-knockout lines (Supplementary Figure S3A).

Gene Ontology (GO) enrichment analysis of DEGs showed significant clustering in the cellular component (CC), molecular function (MF), and biological process (BP) categories, including small molecule metabolism, extracellular region, peptidase regulation, and oxidoreductase activity (Supplementary Figures S5A–C). These indicate roles in cellular metabolism, the encoding of secreted proteins or extracellular matrix components, extracellular signal transduction or structural support, and protein degradation. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis identified enrichment in PPAR signaling, phagosome, vascular smooth muscle contraction, and ECM–receptor interaction (Supplementary Figures S5D–F), linking multiple biological processes, including lipid metabolism, immune response, cardiovascular regulation, cell adhesion/migration, and signal transduction. Surprisingly, in the TGFb signaling pathway, except for the downregulation of bmp1b gene expression, other members were upregulated; the target genes of the TGFb signaling pathway, except for gatad2b, were also upregulated (Supplementary Figure S6).

Given the consistent cardiac dysfunction and valve development defects across all three mutant strains, Venn analysis identified 356 commonly dysregulated genes (Figure 5A), which were compiled into the Venn_356 gene subset (Supplementary Table S3).

Heatmap clustering confirmed consistent expression patterns (Figure 5B), with KEGG enrichment in “Cardiac muscle contraction” and “ECM–receptor interaction” (Figure 5C), which are directly relevant to contractile dysfunction and valve development defects, respectively. GO analysis further validated the overrepresentation of the CC, MF, and BP categories (Figure 5D).

Functional and molecular analyses converged to show that all three mutant lines (bmpr2a ^−/− ;bmpr2b ^+/+ , bmpr2a ^+/+ ;bmpr2b ^−/− , and bmpr2a ^−/− ;bmpr2b ^−/− ) exhibited heart failure at 48 hpf because of impaired contractility (Figure 2). Transcriptomic analysis revealed consistent upregulation of cardiac contraction genes, validated by qRT-PCR, in both whole embryos and isolated cardiac tissue for cardiac contraction genes myl9a, myl9b, tnnc1a, cmlc1, and myl7 and the heart failure marker nppa (Figures 6A,B; Supplementary Figure S7A). WISH for nppa and myl7 showed elevated expression and ventricular dilation trends in the mutant strain compared to those in the WT (Figure 6C). Furthermore, Western blotting analysis using protein extracts from zebrafish cardiac tissue at 48 hpf confirmed that the protein level of Myl7 was upregulated in the bmpr2a/b-knockout zebrafish line (Figure 6D). Collectively, these data demonstrate that the bmpr2a/b loss-of-function mutation disrupts cardiac contraction through transcriptional dysregulation of contraction-related genes.

The results presented in Figure 4 demonstrate that bmpr2 insufficiency disrupted valve morphogenesis in zebrafish. This was evidenced by significant alterations in the expression of key valve markers, where fn1b, has2, and nfatc1 were upregulated, while klf2a was downregulated in all three mutant lines (bmpr2a ^−/− ;bmpr2b^+/+ , ;bmpr2a^+/+ ;bmpr2b ^−/− , and bmpr2a ^−/− ;bmpr2b ^−/− ). These changes in expression were further confirmed by heatmap analysis (Figure 7A).

The ECM is critical for valve remodeling (Kern, 2021), and it is associated with valve disease (Huang et al., 2022). The KEGG enrichment analyses of DEGs between WT and bmpr2-knockout zebrafish highlighted significant changes in the ECM–receptor interaction pathway genes across all three mutant zebrafish groups (Figure 5C, red box). Heatmap analysis and qRT-PCR assay conducted on both whole embryos and isolated cardiac tissue consistently demonstrated the upregulation of multiple ECM–receptor interaction-related genes, including vtnb, vtna, ambp, itga1, lamb2, fn1a, tpbg, cspg4, and vcanb, along with the downregulation of aggf1 (Figures 7A–C, Supplementary Figures S7B, C). Furthermore, Western blotting analysis using protein extracts from zebrafish cardiac tissue at 48 hpf confirmed that the protein level of Fn1 (a key ECM component) was upregulated in the bmpr2-knockout zebrafish line (Figure 7D).

Zebrafish valve development initiates at 48 hpf, coinciding with the relative upregulation of ECM (Steed et al., 2016). During ECM remodeling, endocardial cells undergo endothelial-to-mesenchymal transition (EndMT) and subsequent post-EndMT processes, in which mesenchymal cells differentiate into valve interstitial cells, thereby driving valve elongation (Coram et al., 2015). Consistent with these developmental dynamics, our heatmap data showed that bmpr2a/b-knockout zebrafish exhibited upregulated EndMT-related genes (snail, cdh1, and cdh17) and downregulated mesenchymal cell-related genes (cdh2 and vim) (Figure 7A). These results demonstrate that the EndMT process is potentially dysregulated during valvular development in bmpr2a/b-knockout zebrafish.

Although Cdh1 (E-cadherin) is not a canonical EndMT marker, its downregulation is a well-recognized hallmark of cell–cell adhesion loss during both epithelial-to-mesenchymal transition (EMT) and EndMT, reflecting the shift of the cells from a tightly connected epithelial or endothelial phenotype to a highly migratory mesenchymal state (Ma et al., 2020; Singh et al., 2024). Therefore, Cdh1 can be considered an epithelial/endothelial-related gene whose downregulation indicates the progression of EndMT, and we used it to evaluate EndMT dynamics. The above expression pattern was observed in both whole embryos and isolated cardiac tissue (Figures 7E, G, Supplementary Figures S7D, E). Furthermore, Western blotting analysis of protein extracts from zebrafish cardiac tissue at 48 hpf confirmed the transcriptional trends at the protein level; the protein level of the EndMT-related protein Cdh1 was upregulated, whereas the protein levels of Cdh2 and Vim (mesenchymal cell-related genes) were downregulated in the bmpr2-knockout zebrafish line (Figures 7F,H). Protein–protein interaction (PPI) network analysis predicted that ECM–receptor genes (vtnb, vtna, ambp, itga1, lamb2, fn1a, tpbg, and cspg4) act upstream of the EndMT regulators (snail2, cdh17, and cdh1) and mesenchymal-related genes (cdh2 and vim), with the valve marker nfatc1 positioned downstream of the EndMT effectors (Figure 8A).

Collectively, these findings revealed that bmpr2a/b inactivation upregulated ECM–receptor interaction signaling, which may regulate the EndMT process to facilitate early valve development in zebrafish (Figure 9). However, the specific molecular mechanisms require further experimental verification.