The sequences of all species (Homo sapiens, Mus musculus, Gallus, Xenopus tropicalis, Danio rerio, Latimeria chalumnae, Oryzias latipes, Oreochromis niloticus, Takifugu rubripes, and Gasterosteus aculeatus) for the phylogenetic analysis were obtained from NCBI (https://www.ncbi.nlm.nih.gov). All the sequence IDs are available in Supplementary Table S1. GSDS2.0 (http://gsds.gao-lab.org/index.php) and MEME (https://meme-suite.org/meme/tools/meme) were used to reveal the homology of gene structure among species. The SWISS_MODEL (https://swissmodel.expasy.org/interactive) was used to predict the protein tertiary structures of bmp7a and bmp7b.

The transgenic zebrafish parent labeled with green fluorescent protein for osteoblast-specific transcription factor (Osterix GFP) used in this experiment was provided by Professor Xiao Chongde from Chung Yuan Christian University, Taiwan. Fish were kept at 26–28°C under a controlled light cycle (14 h light, 10 h dark) to induce spawning. CRISPR-Cas9 genome editing technology was used to target zebrafish bmp7a and bmp7b. The sequencing of gRNAs and PCR primers are listed in Supplementary Table S2. The mutations of bmp7a and bmp7b were predicted to cause a translational frame shift and induced premature transcription termination.

The zebrafish embryos were fixed with 4% paraformaldehyde for 24 h. Fragments (length between 300 and 600 bp) of bmp7a and bmp7b were cloned using the primers shown in Supplementary Table S2. The PCR products were purified, and then the whole-mount in situ hybridization probes were synthesized in accordance with the manufacturer’s instructions using the DIG RNA Labeling Kit (SP6/T7) (Roche Diagnostics GmbH, Mannheim, Germany, SKU: 11175025910). The steps of whole-mount in situ hybridization were performed according to the method described by Thisse and Thisse (2008). Anti-DIG-AP,Fab fragments (Roche, SKU: 11093274910) from Sigma-Aldrich were used as antibodies.

The eye tissues at 30 days post-fertilization (dpf) and skin tissues of wild types and bmp7b mutants at 60 dpf (3 samples each) were fixed in 4% paraformaldehyde for 24 h. After that, the samples were dehydrated through a series of graded ethanol solutions (70–100%), embedded in paraffin, and sectioned at 5 μm thickness. These sections were stained with hematoxylin and eosin stain; then analyzed and photographed using an Olympus BX53 Upright microscope (Olympus, Beijing, China).

Bone staining with alizarin red was performed as previously described (Nie et al., 2019) for 90 dpf wild types and bmp7b mutants (10 samples each). The stained specimens were analyzed and photographed under an Olympus SZX2 Stereo microscope (Olympus, Beijing, China).

To study the change in melanin content, we stripped the skin (0.2 g) of wild type and mutant fish (20 samples each) after anesthetizing, and extracted the melanin referring to a modified method (Yuan et al., 2008). A UV spectrophotometer UV-3100 (Mapada, Shanghai, China) was used to measure the photometric value of the solution containing melanin. Melanin content was expressed as A₂₉₀ nm/g (wet weight).

Observations of feeding behavior was performed according to previously study (Miller et al., 1993). Fish were starved for 12 h before the experiments. Wild types and bmp7b mutants at 90 dpf were separately transferred to an experimental tank (length 30 cm × height 24 cm × width 12 cm filled with 2 L of clean water) and allowed to acclimatize to the environment for at least 1 h. On one side of the tank, 5 ml of brine shrimp (Artemia salina) was added at a speed of 15 ml/min to stimulate the fish. A digital video format (Sony, Shanghai, China) was used to record videos of 5 min before and after adding brine shrimp, and EthoVision XT 14 software (Noldus Information Technology, Beijing, China) was used to analyze the behavioral changes (speed, swimming route and reaction time) between the wild-type and mutant fish. The wild type and mutants each had six biological replicates. According to report on zebrafish behavior from Fero et al. (2011), we tested the photosensitivity of the wild-type and mutant fish. Wild types and mutants (all 5 dpf) were separately transferred to an experimental tank (diameter 35 mm, height 10 mm) and allowed to acclimatize to the environment for at least 15 min in the light. Then the light was turned off and the behavior of the fish was recorded. All experiments were performed in EthoVision XT 14 software and there were six samples for wild-type and mutant fish, respectively. During the experiments, the room and water temperatures were maintained at 25 ± 1°C.

Total RNA was isolated from the eyes and skin of 30 dpf wild types and bmp7b mutants with three biological replicates using RNAiso Plus Reagent (TaKaRa, Beijing, China) according to the manufacturer’s protocol. Each sample consisted of tissues from six fish. Subsequently, total RNA was quantified using a 2100 Bioanalyzer system (Agilent Technology, CA, United States) and the samples that met the quality requirements were used for RNA-seq library preparation. After RNA quality and quantity control, 12 transcriptome sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, Beijing, China) following the manufacturer’s recommendations, and index codes were added to attribute sequences to each sample. After cluster generation, the library preparations were sequenced on an Illumina Hiseq platform. After quality control of raw data, paired-end clean reads were aligned to the zebrafish reference genome (ftp://ftp.ensembl.org/pub/release-92/fasta/danio_rerio/) using Hisat2 v2.0.4 (Kim et al., 2015). HTSeq v0.9.1 was subsequently used to count the read numbers mapped to each gene, and then the FPKM of each gene was calculated (Trapnell et al., 2010). Differential expression analysis of two conditions/groups (three biological replicates per condition) was performed using the DESeq R package (1.18.0) (Wang et al., 2009). Genes with a p-value < 0.05, adjusted using the Benjamini and Hochberg’s approach for controlling the false discovery rate, were defined as differentially expressed genes (DEGs). KEGG enrichment analysis of DEGs was performed and KOBAS software was used to test the statistical enrichment of DEGs in the KEGG pathways (http://www.genome.jp/kegg/) (Mao et al., 2005).

The NCBI primer tool (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) was used to design primers for the tested genes (rgs9a, rgs9b, rcvrn2, guca1d, grk1b, opn1mw4, gc2, wnt7ba, gna14, erbb3b, bambia, bmp7b, β-actin, gh1, ghra/b, igf-1/2a/2b, igf1ra/b, and igf2r). We used the PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, code No. RR047A) to synthesize the cDNA. The qRT-PCR experiments were performed using TB Green^® Premix Ex Taq™ II (Tli RNaseH Plus) (TaKaRa, code No. RR820A) and the QuantStudio™ ⁶ Flex qRT-PCR system (ABI, Foster City, CA, United States). Thermal cycling conditions were as follows: 5 min at 95°C, followed by 40 cycles at 95°C for 20 s, and at 60°C for 25 s. Using the β-actin gene of zebrafish as a reference gene, we used the comparative 2^−△△Ct method to determine the relative gene expression between the wild types and mutants (Livak and Schmittgen, 2001). IBM SPSS Statistics software (version 22.0) was used to perform a Student’s t-test to compare the differences between the two groups. Statistical significance was set at p < 0.05, and statistical high significance at p < 0.01. The fold-change values are the average of three biological replicates in each group. The primer sequences are listed in Supplementary Table S2.

The evolutionary origin and phylogenetic distribution of BMP7 showed that BMP7 orthologs, which evolved from the same ancestral gene family of vertebrates, had duplications in teleosts after the evolution of L. chalumnae, named bmp7a and bmp7b. The bmp7b clade was clustered into one group of humans and mammals (Figure 1A). The gene sequence and structure analysis of BMP7 genes showed that the zebrafish bmp7b gene is highly homologous to the BMP7 gene in humans and mammals. Although the exon numbers of these genes were the same, zebrafish bmp7b (100.69 kb) was closer in length to human (97.87 kb) and mouse (72.31 kb) BMP7 genes, and was significantly larger than zebrafish bmp7a (31.76 kb) (Figure 1B). In addition, these species have a common sequence encoded by exon 6 and exon 7, called the TGF-β domain (103aa). The TGF-β domain-encoding sequences were highly conserved across different species (Figure 1C).

To explore the role of BMP7 gene in fish species and the differences between orthologous genes, we constructed bmp7a and bmp7b mutant zebrafish using CRISPR-Cas9 technology. One bmp7a nonsense allele with a 19 bp deletion in the first exon and two bmp7b nonsense alleles with 49 and 25 bp deletions in the first exon were generated. The bmp7a and bmp7b nonsense alleles all led to premature translation termination (Figure 1D). Additionally, protein structure prediction showed that the TGF-β domains of bmp7a and bmp7b were destroyed because of the missing sequence and could not form a complete protein structure (Figures 1E,F). Since the phenotypic consistency of different mutant lines of the same gene was confirmed in subsequent experiments, the homozygous mutant lines of bmp7a and bmp7b were referred to as bmp7a ^-/- and bmp7b ^-/-, respectively, unless otherwise stated.

We did not find any mutants when we tested the 48 hpf embryos from the bmp7a heterozygote cross. Therefore, we suspected that bmp7a deletion is lethal and would cause early embryonic death in the bmp7a ^-/- mutants. To confirm this, we detected embryos at different stages of development, and all genotypes containing a homozygous knockout of bmp7a were found to be lethal at 12–24 hpf. After all the surviving embryos were tested, no homozygous embryos were found to survive at 48 hpf (Table 1). To elucidate how a homozygous bmp7a deletion causes early embryonic death, we examined the morphology of the wild-type and bmp7a ^-/- embryos. After a large-scale screening of the 501 embryos from three pairs of heterozygous parents, we counted the number of three genotypes. The proportions (bmp7b ^+/+: bmp7b ^+/-: bmp7b ^-/- = 24.55%: 48.70%: 26.75%) of genotypes met Mendel’s laws of inheritance. Moreover, we simultaneously observed the morphology of these embryos and found that bmp7a ^-/- embryos had defects in the dorsal-ventral pattern. The bmp7a ^-/- embryos were morphologically distinguishable as distinctly ovoid-shaped embryos at 11 hpf. In the initial stage of somite formation (13–15 hpf), all somites of the mutants extended to the ventral side, whereas the tail did not develop to the ventral side. At the 10-somite stage (16 hpf), the ventral fusion of somites caused the yolk cells to rupture and eventually led to death (Figure 2A). In addition, the tails of bmp7a ^-/- mutants showed severe deformities at 16 hpf. Approximately, 34% of the bmp7a ^-/- mutants developed their tails to the ventral side, 66% to the dorsal side (Figure 2B). Moreover, we investigated the spatiotemporal expression pattern of bmp7a using whole-mount in situ hybridization. In 12 hpf wild-type embryos, bmp7a was expressed on the entire blastoderm, especially on the ventral side. In 14 hpf wild-type embryos, the head and tail ends were significantly expressed, but no expression was observed in the dorsal area. In 24 hpf wild-type embryos, bmp7a was significantly expressed in the head and abdomen, but not in the tail. However, in the mutant embryos, the expression of bmp7a disappeared at alive time points (Figure 2C).

We studied the effect of the bmp7b gene on bone development by staining adult fish with whole-mount alizarin red. The results showed that approximately 80% of adult bmp7b ^-/- (n = 20) had bend haemal arches and nodules in the intermuscular bones (Figure 3A). Compared to the wild type fish, the external morphology of bmp7b ^-/- mutants was significantly different. That is, the growth was significantly slower, and these mutants had a darker body surface color (Figure 3B). To compare the effects of bmp7b deletion on growth and development, we measured the growth data (body length and weight) of the bmp7b ^-/- and wild-type zebrafish at three time points (30, 60, and 90 dpf). The results showed that the body length and weight of the bmp7b ^-/- zebrafish were significantly lower than those of wild types (p < 0.05) (Figure 3C).

Fish growth is mainly regulated by the hypothalamus-pituitary-liver growth axis (GH/IGF axis). To explore the causes for the differences in growth between the wild-type and bmp7b ^-/- zebrafish, we compared the expression levels of genes related to growth (growth hormone 1, gh1; growth hormone receptor a/b, ghra/b; insulin-like growth factor-1/2a/2b, igf-1/2a/2b; insulin-like growth factor 1 receptor a/b; igf1ra/b; and insulin like growth factor 2 receptor; igf2r) at different points (30, 60, and 90 dpf). The results showed that in the bmp7b ^-/- mutants, the expression of gh1 decreased significantly at all three points, and ghra and igf2b decreased significantly at both 60 and 90 dpf (p < 0.05); however, igf-1 and ghrb were increased significantly at 30 dpf, and igf1rb was increased at both 60 and 90 dpf (Figure 3D).

We were intrigued by the darker body color of bmp7b ^-/- mutants; therefore, we stained the transverse section of the skin at 30 dpf with HE. We compared the structure of the stripes in bmp7b ^-/- mutants and wild types in four parts (Figure 4A), and found that the skin of the bmp7b ^-/- mutants had a higher number of melanocytes with a wider distribution (Figure 4B). In addition, the melanin content in the skin of adult bmp7b ^-/- mutants and wild types were measured. The melanin content in the skin of wild types was 32.987 ± 0.477 A₂₉₀/g, and that in the bmp7b ^-/- mutants was 42.800 ± 0.904 A₂₉₀/g, indicating a significant difference between the both (p < 0.05) (Figure 4C). These results confirmed that the increase in melanin content in the skin of bmp7b ^-/- mutants led to the darkening of the body color.

To further explore the causes for the difference in body color, we compared the transcriptomes of the skin of wild types and bmp7b ^-/- mutants at 30 dpf using RNA-seq. After quality control of the original data, reference genome mapping, and gene expression level quantification, a total of 1128 DEGs between the two groups, including 627 upregulated genes and 501 downregulated genes in mutants, were identified based on an adjusted p-value < 0.05 (Figure 4D, Supplementary Table S3). Subsequently, KEGG pathway analysis revealed important signalling pathways that were significantly enriched in DEGs, such as the mitogen-activated protein kinase (MAPK) (KO: dre 04010), Wnt (KO: dre 04310), melanogenesis (KO: dre 04916), and tyrosine metabolism (KO: dre 00350) signaling pathways, associated with melanin production (Figure 4E, Supplementary Table S4). The upregulation of gna14 and downregulation of wnt7ba, erbb3b, and bambia might promote an increase in melanin production in the skin (Figure 4F). The expression of candidate genes was verified by qPCR (Figure 4G), and the results were consistent with the transcriptome data.

The results of whole-mount in situ hybridization of wild-type embryos showed that bmp7b was significantly expressed in the eyes in 24 hpf embryos, confirming that bmp7b may have a role in eye development in fish. Additionally, we also observed that the expression of bmp7b in the mutants was significantly reduced (Figure 5A). The eyes of bmp7b ^-/- mutants showed a larger iris (Figure 5B), and the ratio of the iris diameter to the eye diameter also increased significantly (p < 0.05) and gradually increased with age (Figure 5C).

Furthermore, we observed the histological characteristics of the eyes in wild-type and bmp7b ^-/- zebrafish at 60 dpf by longitudinal tissue sections and found an increase in the thickness of the pigment cell layer of the bmp7b ^-/- mutants (Figures 5D,E). Therefore, we speculated that the changes in the retinal structure caused by bmp7b deletion may impact the vision of bmp7b ^-/- mutants. To verify this hypothesis, we analyzed the feeding behavior of bmp7b ^-/- mutants by observing their reaction after the addition of food (brine shrimp). The response of the wild types to food was instant while that of the bmp7b ^-/- mutants was more sluggish (Figure 5F), consistent with their trajectory behavior. After adding brine shrimp, the wild types quickly found it and swam towards it, while the bmp7b ^-/- mutants took a longer time to discover it. Additionally, we also analyzed the behavior 5s after adding brine shrimp and the sensitivity to food and the moving distance of bmp7b ^-/- mutants was significantly lower than that of the wild types (p < 0.05) (Figure 5G). In order to further verify whether the eye development of the mutant was affected, we examined the angular velocities of wild type and bmp7b ^-/- mutants within 5 s after the light intensity was reduced (from light to dark). The results showed that the angular velocity of the mutant was significantly reduced (p < 0.05) (Figure 5H).

The transcriptomes of the eye tissues between the wild-type and bmp7b ^-/- zebrafish at 30 dpf were compared using RNA-seq. Transcriptomic analysis indicated significant differences between the wild-type and bmp7b ^-/- zebrafish (Figure 6A). Over 3,552 DEGs were detected in the eyes of bmp7b ^-/- fish compared to wild types, with 1,090 DEGs were upregulated and 2,462 DEGs were downregulated (Figure 6B, Supplementary Table S5). KEGG analysis showed that 20 enriched pathways for the DEGs between wild-type and bmp7b ^-/- mutants were clustered. These included six pathways related to diseases or immune responses (Supplementary Table S6), including apoptosis (KO: dre00010), necroptosis (KO: dre04217), phagosome (KO: dre04145), cell cycle (KO: dre04110), cellular senescence (KO: dre04218), and herpes simplex infection (KO: dre05168), suggesting that inflammatory and eye or retina-related diseases appeared during the development of the eye. In addition, changes in the structure of the eye and disease pathways may lead to the changes in phototransduction (KO: dre04744), which has been associated with dramatic morphological, structural, and functional changes in eye or retinal development. This resulted in a decrease in the vision of the mutants, which was also consistent with the behavioral results. Wnt (KO: dre04310), TGF-β (KO: dre04350), tyrosine metabolism (KO: dre00350), and MAPK signaling pathways (KO: dre04010) were all involved in the production and development of melanin (Figure 6C). Rgs9a, rgs9b, rcvrn2, guca1d, grk1b, opn1mw4, and gc2 genes might be associated with the light-sensing ability of the eye, while wnt7ba, gna14, erbb3b, and bambia genes might promote an increase in melanin in the iris of the eyes. The expression of candidate genes was verified by qPCR (Figure 6D), and the results were consistent with the transcriptome data.