All experiments strictly abided by the experimental animal guidance policy of the experimental animal ethics committee of the Institute of Oceanology, Chinese Academy of Sciences (Permit Number: IOCAS2022012).

O. punctatus were obtained from the Fish Breeding Base of The Institute of Oceanology, Chinese Academy of Sciences (Haihe Aquatic Seedling Co., Ltd., Weihai Wendeng). To detect the differential expression of BMP genes among the tissues of O. punctatus, we anaesthetized three adult O. punctatus with MS-222, dissected and collected eight different tissues (gills, liver, spleen, brain, the surrounding tissue of the beak-like teeth, intestine, muscle, and heart), and preserved the tissue samples in liquid nitrogen for further analysis.

Three-year-old Oplegnathus samples were selected, their teeth were dissected, and the surface structure of the beak-like teeth was observed under a ZEISS Stemi 2000-c stereoscope. Otherwise, samples were imaged with high-resolution, high-contrast scans by a ZEISS Xradia 515 Versa 3D X-ray microscope at 20.0 and 3.0 µm voxel resolutions. Afterwards, ZEISS 3D Viewer software was used to observe the reconstructed 3D structure of the beak-like teeth of Oplegnathus, and the teeth were rendered in 3D using Dragonfly software. Samples of O. punctatus juveniles were collected for skeletal staining observation. Bone staining was performed according to the methods described by Dingerkus and Uhler (1977). Samples were stored in glycerol after processing, and a small amount of thymol was added to prevent bacterial contamination. Bone-stained O. punctatus samples were placed in a Petri dish with glycerol, observed and photographed with a ZEISS Stemi 2000-c stereoscope.

All available BMP gene sequences and BMP amino acid sequences of nine species (Cyprinus carpio, Salmo salar, Danio rerio, Oryzias latipes, Xenopus laevis, Takifugu rubripes, Gallus gallus, Mus musculus, and Homo sapiens) were downloaded from the NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi) and Ensembl (http://asia.ensembl.org/) public databases. A simple HMM search in TBtools software and NCBI Blast were used to search the reference genomes of O. fasciatus and O. punctatus and preliminarily obtain candidate BMP gene sequences of these two species (Chen et al., 2020). Meanwhile, the NCBI-CDD database (https://www.ncbi.nlm.nih.gov/cdd/) was used to verify the identified characteristic domains of BMP genes in O. fasciatus and O. punctatus. The SMART (http://smart.embl-heidelberg.de) website was used to predict the structural information of BMP protein domains, such as TGFb and TGF_β. The mass fractions and isoelectric points of the BMP proteins were estimated by the ExPASy online tool (https://web.expasy.org/compute_pi). The reference genome accession codes used for analysis in this study are as follows: O. fasciatus (NCBI accession code: PRJNA393383) and O. punctatus (CNGB accession code: CNP0001488) (Xiao et al., 2019; Li et al., 2021).

The BMP protein sequences of O. fasciatus, O. punctatus and the other nine species (C. carpio, S. salar, D. rerio, O. latipes, X. laevis, T. rubripes, G. gallus, M. musculus, and H. sapiens) were aligned using MUSCLE of TBtools (Edgar, 2004). Then, the created alignment file was used for phylogenetic analysis with MEGA11 software to construct a neighbour-joining (NJ) tree and a maximum likelihood (ML) tree (1,000 bootstraps) for comparison (Tamura et al., 2021). The online software EvolView (https://evolgenius.info//evolview-v2/) was used to beautify the evolutionary tree. Multiple sequence alignment of O. fasciatus was performed using Jalview software. MEME (http://meme-suite.org/tools/meme) online software was used to analyse BMP motifs. Then, the results were visualized with TBtools software. The conserved domains of BMPs in Oplegnathus were obtained using CDD of the NCBI and visualized by TBtools. Chromosome position and gene density information were obtained through the Oplegnathus genome database, the BMP gene positions were marked on the chromosomes, the collinear relationships between different chromosomes were obtained through the genome database, and the results were visualized using TBtools software.

Phylogenetic trees were constructed individually for BMP gene families from 11 species using the NJ method in MEGA11. Using the CODEML program of PAML4.9j (Yang, 2007), selection pressure analysis was performed on the main clades using different models (branch model, site model and branch-site model). First, branch models are used to determine whether the selection pressure on each branch of the evolutionary tree is significant (Yang, 1998). The main branch models are the one-ratio model (M0), in which the ω-values of all branches in the phylogenetic tree are equal; the free-ratio model (M1), in which the ω-values among branches in the phylogenetic tree are unequal; and the two-ratio model (M2), where the ω-values of foreground branches are different from those of background branches. In addition, the site model assumes that different branches of the phylogenetic tree are subject to the same selection pressure but that different amino acid sites experience different selection pressures. The following are the main site models: the one ratio model (M0), assuming that all sites have the same ω value; the near-neutral model (M1a), assuming that only conserved sites (0 < ω < 1) and neutral sites (ω = 1) are present; the positive selection model (M2a), in which sites with positive selection (ω > 1) are added to M1a; the discrete model (M3), assuming a simple discrete distribution trend of ω values for all sites; the beta model (M7), assuming that ω belongs to the matrix (0, 1) and follows a beta distribution for all sites; and the beta and ω model (M8), in which another class of ω values is added to M7, which can be obtained computationally and can be set to 1. Comparisons between M1 and M2, M0 and M3, and M7 and M8 are often used for site models, and M1 and M2 are more stable than M7 and M8. Moreover, in the analysis of the branch-site model, we tested for positively selected amino acid sites among the branches through the settings of different foreground branches and compared the null hypothesis of each Model A (MA) with the corresponding alternative hypothesis (Zhang et al., 2005). Then, the ratio ω of non-synonymous to synonymous substitutions (dN/dS) was calculated. When ω = 1, selection is neutral (that is, there is no selection); when ω > 1, there is positive selection; and when 0 < ω < 1, there is negative selection, also called purifying selection. In this paper, the likelihood ratio test (LRT) was used to compare model pairs, and the Bayesian method (BEB) was used to identify the sites with positive selection (Bielawski and Yang, 2003).

RT–PCR was performed with the CFX96 real-time detection system (Bio-Rad, United States) to determine BMP expression differences in different developmental stages and tissues of O. punctatus. Total RNA was isolated from samples of eight different tissues (gill, liver, spleen, brain, the surrounding tissue of the beak-like teeth, intestine, muscle, and heart) and five developmental stages (5 days past hatching (dph), 16, 30, 50, and 60 dph) in O. punctatus using TRIzol (Invitrogen) according to the manufacturer’s instructions. Considering that the fish were too small to precisely separate the surrounding tissue of the beak-like teeth, the sampling locations for the five developmental stages were on the head, except for days 5 and 16 dph, when whole fish were taken because they were too small to precisely separate the head. RNA concentration and purity were measured using an ultra-micro spectrophotometer (NanoDrop 2000, United States). In general, the OD260 nm ratio was OD280 nm ≥ 1.8, and the OD260 nm ratio was OD230 nm ≥ 1. RT-PCR was conducted to synthesize cDNA from 1 μg of total RNA by using an Evo M-MLV RT Kit (Accurate Biology, Changsha, China) following the manufacturer’s instructions. The internal positive control was the β-actin gene from O. punctatus, with the forward primer (5′-GCT​GTG​CTG​TCC​CTG​TA-3′) and reverse primer (5′-GAG​TAG​CCA​CGC​TCT​GTC-3′). Details of the BMP gene primer pairs are provided in Supplementary Table S1. qRT–PCR experiments were performed with a 20 μl mixture consisting of 0.4 μl of each primer, 10 μl of SYBR^® Premix Ex Taq™ II (Tli RNaseH Plus) (2X), 7.2 μl of ddH₂O, and 2 μl of template cDNA. The PCR program consisted of an initial denaturation step at 95°C for 30 s and 39 denaturation cycles: denaturation at 95°C for 5 s and annealing at 60°C for 30 s. The temperature depended on the primers. Melting curve analysis was performed at the end of the reaction to demonstrate its specificity. Each experiment was performed in triplicate. β-Actin was used for normalization of all gene expression data (Wang et al., 2020). According to the measured Ct value (the average value of three parallel samples for Ct), the relative expression levels of the BMP genes were calculated by Pfaffl’s method (Pfaffl, 2001). All data obtained in the experiment were analysed by one-way analysis of variance using GraphPad Prism 9 software, and p < 0.05 was indicative of a statistically significant difference. Quantitative results were plotted using the R environment and GraphPad Prism 9 software.

As shown in Figure 1, the Oplegnathus species have a unique healing beak-like teeth structure, with the upper and lower jaw teeth showing healing and a gap in the middle. In addition, there are obvious traces of the accumulation of replacement teeth on the outer edges of the beak-like teeth. Moreover, the outer layer of the beak-like teeth is filled with smooth calcium, while the inner layer is relatively rough. As shown in Figure 1, the 3D rendering analysis of the thickness of the teeth by micro-CT revealed a nested internal arrangement of the teeth of the Oplegnathus species, extending from far from the rostral end to near the rostral end. In addition, several irregularly distributed internal teeth could be observed on the inner edge of the beak-like teeth (Figure 1E). The micro-CT scan results at 30, 50 and 70 dph showed beak-like teeth distributed in a single row of separation at 30 dph, replacement teeth at 50 dph, beak-like teeth that generally tended to heal, and complete rostral tooth healing at 70 dph (Figures 1G–I). The results of bone staining are shown in Figures 1J–L. At 45 dph, the beak-like teeth are individually arranged conical teeth, and most of the teeth and bones are stained blue with bromophenol blue, showing a cartilaginous state; only the tips of the teeth, which are beginning to ossify, are red. After 50 dph, the discrete arrangement of conical teeth shows a tendency to heal, and all of the beak-like teeth are stained with alizarin red, indicating that ossification has been completed. By 90 dph, the beak-like teeth of O. punctatus had healed. External observation showed clear signs of replacement tooth accumulation on the outer edges of the upper and lower teeth, with three to four rows of replacement teeth.

The cDNA sequence characteristics and accession numbers of all cloned BMP genes of Oplegnathus are shown in Table 1. The accession numbers of all cloned BMP genes of the other 11 species are shown in Supplementary Table S2. 19 and 16 different BMP genes were found in the genomes of O. fasciatus and O. punctatus, respectively (Table 1). In Table 1, genes with the prefix “Of” before “BMP” are BMP genes in O. fasciatus, and those with the prefix “Op” are BMP genes in O. punctatus. The length of the coding sequence (CDS) is 966–1,971 bp, the length of the encoded protein is 321–656 amino acids, the molecular weight is between 36 and 73 kDa, and the isoelectric point is between 5.09 and 9.72. As shown in Figure 2, to deeply investigate the structural diversity of the BMP gene family, the conserved domains and motifs of the Oplegnathus BMP proteins were visualized by TBtools software. Except BMP1, the Oplegnathus BMPs have similar TGF_β domains, indicating that the BMP gene family belonging to the TGF-β superfamily is highly conserved. In addition, most BMP genes, except BMP3 and BMP15, also possess a TGFb propeptide domain, revealing their potentially differentiated functions in biological processes. Furthermore, a total of 10 distinct conserved motifs were identified. The conserved domains and motifs were compared according to the phylogenetic relationships between Oplegnathus and T. rubripes (please refer to Supplementary Figure S1 for a description of the names and lengths of the motifs). The highly consistent conserved domains and motifs in the same subfamilies indicated similar protein functions among members of the same subfamily. To explore the conserved domains and phylogenetic relationships of the BMP proteins of O. fasciatus and O. punctatus, multiple sequence alignment of their TGF_β domains was performed, revealing that all BMP family members contain a highly conserved TGF_β region at the C-terminus, which consists of approximately 100–110 amino acids (Figure 3). As shown in Figures 4A, 19 BMPs were distributed among 15 of the 24 chromosomes in O. fasciatus. Both BMP3b-1 and BMP9 were located on chromosome 10 and had tandem duplications. As shown in Figures 4B, 16 BMPs were distributed among 13 of the 24 chromosomes in O. fasciatus. Both BMP3b-1 and BMP9 were located on chromosome 22. As shown in Figure 5A, the homologous fragment of the OfBMP3b-1 gene located on chromosome 18 of O. fasciatus is the OfBMP3b-2 gene located on chromosome 12. Moreover, the homologous fragment of the OfBMP7a gene located on chromosome 8 of O. fasciatus is the OfBMP7b gene located on chromosome 4. Furthermore, the homologous fragment of OfBMP13 on chromosome 3 of O. fasciatus is located on chromosome 21. As observed in O. fasciatus, the homologous fragment of the OpBMP3b-1 gene located on chromosome 22 in O. punctatus is the OpBMP3b-2 gene located on chromosome 15. Additionally, the homologous fragment of the OpBMP7a gene located on chromosome 3 in O. punctatus is the OpBMP7b gene located on chromosome 7. Furthermore, the homologous fragment of OpBMP13 on chromosome 13 of O. punctatus is located on chromosome 9 (Figure 5B).

To demonstrate the phylogenetic relationships of BMP genes between Oplegnathus and other species, we selected a total of 198 BMP amino acid sequences from 11 species, namely, O. fasciatus, O. punctatus, C. carpio, S. salar, D. rerio, O. latipes, X. laevis, T. rubripes, G. gallus, M. musculus, and H. sapiens. The phylogenetic tree was constructed by the ML method using MEGA11 and TBtools. As shown in Figure 6, the high topological consistency of the phylogenetic relationships of the BMPs indicates that they are highly evolutionarily conserved. Based on structural homology, except for BMP1, the BMPs were divided into six subfamilies, namely, the BMP2/4/16, BMP5/6/7/8, BMP9/10, BMP12/13/14, BMP3/15, and BMP11 groups. Both O. fasciatus and O. punctatus contain three members of the BMP11 gene family: BMP11-1, BMP11-2, and BMP11-3. O. fasciatus has three members of the BMP3 gene family: OfBMP3a, OfBMP3b-1, and OfBMP3b-2. O. punctatus has two members of the BMP2 gene family: OpBMP3b-1 and OpBMP3b-2. The copy number of each BMP gene in nine vertebrates was summarized and analysed. The details of the numbers of BMP members among the 11 species are provided in Table 2. Except in C. carpio, the number of BMP family members in the vertebrates ranges from 14 to 22. C. carpio (Xu et al., 2014; Chen et al., 2017) experienced an additional whole-genome duplication event, giving it approximately double the number of paralogous genes in zebrafish. As shown in Figure 7, BMP1, BMP12, and BMP14 were not found in Oplegnathus or T. rubripes. Additionally, BMP16 has been found only in fish species. BMP13 and 14 are absent in X. laevis, while BMP8 and BMP12 are absent in G. gallus.

In the branch-site model, we set the branch on which the target BMP genes of O. fasciatus and O. punctatus were located as the foreground branch, and the branch on which the BMP genes of the other species were located was set as the background branch. According to the results in Table 3, we identified three genes (BMP3, BMP7, and BMP9) with positively selected sites under the branch-site model. In contrast, no positively selected sites were detected for BMPs in the branch model and site model, and the omega values were less than 1, indicating that the BMP genes are strongly negatively selected under this model. In the branch-site model, for the BMP3 gene, the following 11 amino acid sites were positively selected: 136H, 138K, 140V, 141F, 143F, 145L, 146S, 148I, 150E, 151S, and 153L. Similarly, 559T of BMP7 and 26F of BMP9 were identified as positively selected sites. Furthermore, the LRTs show that the above results are reliable, and the alternative hypothesis Model A is accepted. The specific parameters and results of the model are shown in Supplementary Tables S4–S9. Please refer to the Supplementary Material for details on the model parameters and LRTs for other BMP genes.

The qRT–PCR expression results showed that the BMP genes were expressed in various healthy tissues of O. punctatus, but the expression patterns of these genes varied among the tissues. The gill expression values were selected as a reference for comparison with the expression values of other tissues during the analysis of PCR gene expression data, so the relative expression value for gill tissue was set at 1. In general, most BMP genes were widely expressed, and there was a trend for higher relative expression of most BMP genes in the gill, as this tissue is more likely to be associated with BMP-related pathways. As shown in Figure 8, regarding our focus on the epithelial tissue and mesenchymal tissue of the beak-like teeth, the high relative expression of BMP2, BMP3b-1, BMP3b-2, BMP4, and BMP11-1 in the epithelial tissue and mesenchymal tissue of the beak-like teeth (all values above 1) is interesting. However, the relative expression of BMP3b-2 in the spleen and brain was higher than that in the epithelial tissue and mesenchymal tissue of the beak-like teeth. The relative expression of the remaining BMP genes in the epithelial tissue and mesenchymal tissue of the beak-like teeth was low (all below 0.6). Based on the tissue expression patterns of BMP genes, several BMP genes (BMP2, BMP3b-1, BMP3b-2, BMP4, and BMP11-1) that were highly expressed in the surrounding tissue of the beak-like teeth were screened. As shown in Figure 10, their developmental expression patterns at 30, 50, and 60 dph were also analyzed in correlation with the results of alizarin red staining at early developmental stages. As shown in Figures 8A–C, there are obvious differences in the expression levels of the BMP2/4/16 subfamily. The expression abundance of BMP2 was higher in the gills, liver, intestine and the epithelial tissue and mesenchymal tissue of the beak-like teeth, followed by the spleen and brain, with lower levels in the muscle and heart. BMP3b-1 and BMP4 exhibited similar expression patterns; however, their expression in the liver was significantly lower than that of BMP2. In contrast, BMP16 was highly expressed in the muscle, heart and gill, while the lowest expression was detected in the epithelial tissue and mesenchymal tissue of the beak-like teeth. Unlike BMP3b-1, BMP3b-2 is most highly expressed in the brain, followed by the spleen and gill, and weakly expressed in the muscle and heart. As shown in Figure 8B, the BMP5/6/7/8 subfamily exhibited various expression levels across the examined tissues. There was significant downregulation in the epithelial tissue and mesenchymal tissue of the beak-like teeth for the BMP6, BMP7a, BMP7b, and BMP8 genes. The expression of BMP6 and BMP8 was the highest in the gill. BMP7a was highly expressed in the brain, while BMP7b was highly expressed in the muscle. As shown in Figure 8C, the BMP9/10 subfamily showed higher expression in the gill, muscle and heart. For BMP9, the expression level was the lowest in the liver, followed by the brain, spleen, the surrounding tissues of beak-like teeth, and intestine. The expression of BMP10 was higher in the liver, followed by the spleen, brain and intestine, and the lowest expression was observed in the surrounding tissues of the beak-like teeth. As shown in Figure 8C, members of the BMP11 subfamily were mainly expressed in the gill. Additionally, BMP11-1 was also highly expressed in the intestine and the surrounding tissues of the beak-like teeth, while BMP11-2 and BMP11-3 were highly expressed in the muscle and heart. Furthermore, there was significant downregulation in the liver for BMP11-2.

The temporal expression patterns of BMP genes during different developmental stages of O. punctatus were analyzed by real-time qRT–PCR. The mRNA abundances of different BMP genes in five developmental stages were compared. BMP8 was the most abundantly expressed BMP gene during O. punctatus development. Notably, the relative expression levels of all BMP genes were significantly upregulated from 5 to 30 dph and peaked at 30 dph, with the exception of BMP9. Moreover, the lowest relative expression of all genes was observed at 5 dph. As shown in Figure 10, most of the BMP genes, especially BMP2, BMP3b-1, BMP3b-2, BMP4, and BMP11-1, were abundantly expressed at 30 dph, when the beak-like teeth were conical teeth that had not yet healed and only the tips of the teeth were ossified. The expression of BMPs significantly increased during the period from 16 to 30 dph, which indicates that osteoblasts are rapidly generated at this stage and that activation of ossification in O. punctatus occurs in preparation for the sclerosing of bones and beak-like teeth. The expression of BMP genes across development can be divided into three patterns. BMP2, BMP3b-1, BMP4, BMP11-1, BMP11-2, BMP11-3, and BMP16 exhibited the first expression pattern (Figure 9A): the expression levels increased from 5 to 30 ph, with the highest expression at 30 dph, and decreased sequentially from 30 to 60 dph, but the relative expression at 60 dph was still higher than the starting value (5 dph). In addition, the relative mRNA expression levels of BMP2, BMP3b-1, BMP11-1, BMP11-2, and BMP11-3 showed similar trends at the different developmental stages of O. punctatus. However, the relative mRNA expression levels of BMP4 and BMP16 were significantly higher than those of the above five genes. The second expression pattern (Figure 9B) was different from the first, with the expression levels of BMP3b-2, BMP6, BMP7a, BMP7b, BMP8, and BMP10 increasing from 50 to 60 dph. Compared with other BMP genes, BMP7a showed no significant change in mRNA abundance from 5 to 60 dph. In addition, BMP8 showed a significant increase in mRNA expression at 30–60 dph. For the third expression pattern (Figure 9C), there was a continuous and subtle trend of increasing BMP9 mRNA expression from 5 to 50 dph, with the expression decreasing to a minimum at 60 dph.