Osteocytes as main responders to low-intensity pulsed ultrasound treatment during fracture healing

Ultrasound stimulation is a type of mechanical stress, and low-intensity pulsed ultrasound (LIPUS) devices have been used clinically to promote fracture healing. However, it remains unclear which skeletal cells, in particular osteocytes or osteoblasts, primarily respond to LIPUS stimulation and how they contribute to fracture healing. To examine this, we utilized medaka, whose bone lacks osteocytes, and zebrafish, whose bone has osteocytes, as in vivo models. Fracture healing was accelerated by ultrasound stimulation in zebrafish, but not in medaka. To examine the molecular events induced by LIPUS stimulation in osteocytes, we performed RNA sequencing of a murine osteocytic cell line exposed to LIPUS. 179 genes reacted to LIPUS stimulation, and functional cluster analysis identified among them several molecular signatures related to immunity, secretion, and transcription. Notably, most of the isolated transcription-related genes were also modulated by LIPUS in vivo in zebrafish. However, expression levels of early growth response protein 1 and 2 (Egr1, 2), JunB, forkhead box Q1 (FoxQ1), and nuclear factor of activated T cells c1 (NFATc1) were not altered by LIPUS in medaka, suggesting that these genes are key transcriptional regulators of LIPUS-dependent fracture healing via osteocytes. We therefore show that bone-embedded osteocytes are necessary for LIPUS-induced promotion of fracture healing via transcriptional control of target genes, which presumably activates neighboring cells involved in fracture healing processes.

Fracture healing is a complex, well-orchestrated, regenerative process involving various tissues and cell types. The repair process is divided into four main stages: inflammation, soft callus formation, hard callus formation, and remodeling 1, 2 . The injury initiates an inflammatory response with white blood cells accompanied by the secretion of inflammatory cytokines. This response induces coagulation of the hematoma around the fracture site and is a template for callus formation. A cartilaginous callus is formed by a collagen matrix produced by specific mesenchymal stem cells derived from the surrounding soft tissues and bone marrow. Within the callus, endochondral formation and endochondral ossification occur. The primary soft callus is resorbed by osteoclasts and replaced with a hard callus by osteoblasts. The hard callus is subsequently resorbed by osteoclasts for remodeling into the bone's original cortical structure, characterized by coupled cycles of osteoblast and osteoclast activity. During this period, marrow space is re-established, vascular remodeling takes place, and, finally, new bone is generated. Thus, fracture healing requires the appropriate interaction of various tissues and cell types at each step 3 .
Mechanical stimulation by treatment with low-intensity pulsed ultrasound (LIPUS) is an established therapy for bone fracture treatment, and the US Food and Drug Administration approved EXOGEN (a LIPUS system) in the 1990s for the accelerated healing of certain fresh fractures 4 . Because mechanical loading, generally supplied by gravity and exercise, is effective at inducing bone remodeling and maintaining bone mass due to the regulating biological functions of bone cells 5,6 , it seems reasonable to assume that LIPUS stimulation promotes fracture healing through control of bone cells. Both bone-forming osteoblasts and mechanosensitive osteocytes are capable of responding to LIPUS, and the synergistic action of these cells is likely important for the efficacy of LIPUS treatment 7 . It has been reported that osteoblasts are sensitive to LIPUS and contribute to fracture healing through bone formation and inflammatory regulation 8 . On the other hand, although they are less well studied than osteoblasts, osteocytes are also known to respond to LIPUS and have an important role in the biological function of fracture healing 9,10 .

Results
LIPUS promotes fracture healing in zebrafish but not in medaka. To examine whether osteocytes are necessary for LIPUS to affect fracture healing, we compared fracture healing rates in medaka and zebrafish. Medaka have acellular bone lacking osteocytes, whereas zebrafish have osteocyte-containing bone, but apart from this both types of bone are extremely similar in structure, mechanical properties, and mineral density 15 . The fin rays contain segmented dermal bones and we fractured some of these bones using a scalpel in three places in each fish (Fig. 1a).
Alizarin red and Alcian blue staining was performed to observe the fracture healing processes in both species (Fig. 1b). After the fracture, blood clot forms hematoma which induce inflammatory response and be a template for callus formation. About 1 week, the fractured region formed an Alcian blue-positive callus and was bridged. About 2 weeks, the callus was gradually replaced by an Alizarin red-positive hard callus. About 3 weeks, the Alizarin red-positive calluses were then remodeled into a smooth surface. Since the surface was still not back to normal after a few months, at this point we considered the bone fracture have completely healed. The time frame of each phase of fracture healing, inflammation, chondrogenesis, ossification and remodeling, in medaka The body weight (upper panel) and total length (lower panel) of experimental group of medaka and zebrafish before the LIPUS application (n = 3-5). (e) Upper: Days required for complete fracture healing in LIPUS-treated and untreated medaka and zebrafish. Lower: Promotion rate of fracture healing induced by LIPUS stimulation in medaka and zebrafish (n = 3-5 fishes × 3 fractures, representative data). The time point that the fracture was completely healed was measured in a blinded manner. This experiment was repeated three times with similar results. (f) Upper: The body weight and length of the fish were measured every week, and the condition factor (K-factor; g/cm 2 × 10 3 ) was calculated at the indicated points. Lower: Percent change in condition factor (n = 3-5). Data are presented as mean ± SEM. p Values were determined using Student's t-test. * p < 0.05, ** p < 0.01, *** p < 0.001.The photographs in (a) and (b) were taken by T. S. ◂ and zebrafish are depicted in Fig. 1c. The medaka bone took approximately seven days longer to heal completely than the zebrafish bone (Fig. 1b,c). There was no difference in the process and the period of fracture healing between male and female in both species.
Next, we examined the effect of LIPUS on accelerating bone fracture healing in anosteocytic and osteocytic bone with medaka and zebrafish. In fish fin rays, bone segments in the distal fin were Alcian blue positive unmineralized cartilage (Fig. 1a) because fish fins grow by sequentially adding new segments of bone to the distal end of fin ray 22 . We selected Alizarin red positive mineralized bone located in proximal part of fin ray to induce bone fracture and three fractures per fish were induced (Fig. 1a). Before the LIPUS treatment, medaka and zebrafish were divided into two groups of each and we confirmed that there was no difference in the total length and body weight between two groups (Fig. 1d). Then, the fish divided into two groups were stimulated or unstimulated with LIPUS for 20 min at 30 mW/cm 2 once a day until the fracture was completely healed. Three fractures per fish scored separately because there are about 1-3 days differences in healing time within a fish. The LIPUS treatment accelerated fracture healing in the zebrafish but not in the medaka (Fig. 1e). The condition factor (K-factor: weight [g]/(length [cm] 3 ) × 1000) that is used as an indicator of the physiological state of fish dependent on adequate management (feeding density and climate) was not affected by the LIPUS treatment in either species (Fig. 1f). The LIPUS treatment provided an advantage during fracture healing, probably through direct effects on the skeletal tissue rather than causing systemic changes in physiological condition. Because the effect of LIPUS on fracture healing is apparent in zebrafish bone, which contains osteocytes, but not in medaka bone, which does not, osteocytes appear to be important for LIPUS-mediated fracture healing.
Global transcriptome analysis in LIPUS-stimulated osteocytes. Osteocytes therefore responded to LIPUS stimulation and contributed to the accelerated fracture healing. In order to understand the molecular events in osteocytes triggered by LIPUS, we performed RNA sequencing (RNA-seq) analysis of MLO-Y4 oateocyte-like cells exposed to LIPUS and of controls. Cells were stimulated for 20 min with LIPUS and subsequently placed in a CO 2 incubator for 30 min and RNA was then extracted. The experiment was repeated independently three times and the transcriptomes were determined by using three independent RNA samples. Principal component analysis (PCA) was performed to investigate the variability between the biological triplicates of LIPUS treated and non-treated samples (Fig. 2a). The first component (PC1) that separates two of control samples from LIPUS-treated group represents more than 75% of the variability. However, one of the control biological replicates differ from the other two samples. Since LIPUS treated and non-treated (control) cells were generated in one cell culture plate, leakage of LIPUS stimulation into control cells may lead to the variation among control group. Of 16,476 genes identified, 179 were significantly affected by LIPUS stimulation (Fig. 2b). Of these 179 genes, 93 were upregulated and 86 were downregulated by LIPUS stimulation (Fig. 2b). Since the expression levels of most of the bone-related genes did not change (Fig. 2c), LIPUS-stimulated osteocytes may participate in fracture healing not through their osteogenic function.
To clarify the biological meaning of these 179 genes, functional annotation analysis was performed using the functional annotation tool DAVID and Uniprot keywords were selected as the database. Top-ranked major annotated term categories included Immunity, Antiviral defense, Cytoplasm, Secreted, Extracellular matrix, and Transcription (Fig. 2d). Further, these genes were clustered into nine groups, and enrichment scores greater than 1.3 (which corresponded to the negative log of q-values < 0.05) were considered significant, resulting in seven clusters that were significantly enriched (Table 1). Similarly to the annotation categories, Immunity, Transcription, and Secreted (bolded in Table 1) genes were in the top four result clusters. It seems reasonable that LIPUS modulates these functions in osteocytes, because in fracture healing an inflammatory immune response is necessary and osteocyte-derived secreted factors are likely to play a role in the recruitment of mesenchymal stem cells (MSCs) and blood cells to form a callus and induce vascularization.
These transcription factors may also be important in regulating related gene expression in all fracture healing processes, including the inflammatory and bone remodeling phase. We investigated individual gene expression in three clusters: Immunity, Secreted, and Transcription. 23 of 33 genes genes belonging to the Secreted category were upregulated by LIPUS, whereas 17 of 20 Immunity genes tended to be downregulated (Fig. 2e). In the Transcription category, there were almost equal numbers of up-(12) and downregulated (18) genes after LIPUS stimulation (Fig. 2e). These results show that osteocytes respond to LIPUS via 179 genes that are involved in the inflammatory response, transcriptional regulation, and protein secretion, all of which affect the fracture healing rate.
Gene expression altered by LIPUS in zebrafish fin rays. Next, to examine whether genes in the categories Immunity, Transcription, and Secreted are also modulated by LIPUS in vivo gene expression levels of each cluster category were measured in zebrafish fins. Bone fractures were induced in caudal fins and stimulated with LIPUS for 20 min once a day. On days 1 and 7, RNA was isolated from the fin rays 30 min after stimulation. To perform qPCR with these samples, we isolated the appropriate genes using the steps described below for the Immunity category in cluster A, which included 20 genes (Table 1). From among these genes, we first selected the 14 that had the smallest q-value (0.00388). Next, we eliminated genes that do not have zebrafish orthologs or for which it was not possible to design primer sets for qPCR. Finally, six genes were left (labeled with asterisks in Table 1) and we measured the changes in their expression using qPCR. Similarly, 17 genes were selected from the Transcription category in cluster B and 16 from the Secreted category in cluster D for qPCR (labeled with asterisks in Table 1).
Most of the genes in the Immunity and Secreted categories were not affected by LIPUS stimulation of the zebrafish fins. In the Immunity category, only one out of six genes (16.6%), C1s2 (Complement component 1, s subcomponent 2), was upregulated, and in the Secreted category the expression of only one out of 16 genes www.nature.com/scientificreports/  www.nature.com/scientificreports/ (1.5%), Saa3 (Serum amyloid A-3 protein), was affected (Fig. 3a). However, the expression of 10 out of the 17 genes (58.8%) in the Transcription category was significantly altered (Fig. 3b). Reactivity to LIPUS of the Transcription-classified genes was approximately parallel in MLO-Y4 cells and zebrafish fin rays. Therefore, we focused on Transcription-categorized genes for target gene analysis.

Target gene analysis of transcription factors modulated by LIPUS.
To identify the target genes regulated by the 10 LIPUS-responding transcription factor genes, we utilized the online ChIP-seq database ChIP-Atlas and obtained the binding scores for Egr1, Egr2, nuclear receptor 4a1 (Nr4a1), JunB, and NFATc1 in the promoter regions of target genes. Unfortunately, information regarding Btg2, Snai2, Cited2, FoxQ1, and Helz2 was not available in this database (labeled N/A in Table 2). The top 10 potential target genes with the highest average scores were selected and their expression levels in osteocytic MLO-Y4 cells were evaluated using RNA-seq data (column 4 in Table 2). RNA-seq data showed that all target genes were not changed by LIPUS (column 4-6 in Table 2), probably because there is a time-lag for target genes to be initiated transcription by these transcription genes. Eleven genes which abound with high expression levels in MLO-Y4 cells (FPKM value > 10; bolded in Table 2) were selected and their functions were investigated. These genes participate in energy metabolism (Adssl1), cell survival (Ncapg2, Mcl1, and Brox), and cell-cell contact (Arhgap17), and involve several signaling pathways (Cbr1, Dynlt1b, Aida, and Htra3) ( Table 2). The TGF-β and Cox2/PGE2 pathways, which are modulated by Htra3 and Cbr1 factors 23,24 , are especially interesting because they are induced and secreted in the early stages of fracture healing (inflammation phase) as proinflammatory molecules 3 . Cst3, which is the most abundantly expressed Maff Vdr Cluster www.nature.com/scientificreports/ gene in MLO-Y4 cells, encodes a secreted factor and promotes the differentiation of osteoblasts associated with BMP signaling 25 . Previous reports were referred to for other transcription factor genes (Btg2, Snai2, Cited2, FoxQ1, and Helz2) because their target genes were not available in the ChIP-Atlas database. Btg2 can code for transcriptional coactivators and allows interaction with and modulation of various nuclear receptors, such as all-trans retinoic acid receptors, estrogen receptors, and androgen receptors 26 . It has been reported that estrogen receptors in osteocytes are important for bone formation 27 . Snai2 factors can act as transcriptional repressors and regulate cell apoptosis, migration, and detachment from and attachment to the extracellular matrix 28 . Cited2 encodes a CBP/p300-binding transcription coactivator and enhances TGF-β-mediated transcription 29 . Interestingly, a previous study using a rat fracture model reported that a Cited2 factor was a negative regulator of fracture healing www.nature.com/scientificreports/ www.nature.com/scientificreports/ associated with an increase in matrix metalloprotease 30 . LIPUS stimulation reduced Cited2 expression in the fractured fin rays (Fig. 3b), which is consistent with the positive effect of LIPUS on fracture healing. FoxQ1 encodes a member of the forkhead transcription factor family and regulates cell proliferation via TGF-β and Wnt signaling 31 . Helz2 encodes transcriptional coactivator helicase with zinc finger 2, a coregulator of peroxisome proliferator-activated receptor gamma (PPARγ) 32 . Since PPARγ in osteocytes is important for the coupling of bone formation and resorption 33 , the upregulation of Helz2 by LIPUS treatment (Fig. 3b) may control fracture healing processes such as bone remodeling from osteoblasts and osteoclasts. Taken together, these results show that the 10 transcription genes shown in Table 3 were modulated by LIPUS stimulation not only in vitro but also in vivo. Further, their target genes may be involved in the promotion of fracture healing, probably through activation of the inflammatory response and osteogenic cell differentiation.

Measurement of candidate transcription genes affected by LIPUS in medaka fin rays.
In zebrafish bone, which contains osteocytes, several transcription genes reacted to LIPUS under fracture-healing conditions. To assess whether osteocytes are important for LIPUS efficacy via these transcriptional factor genes, we performed another in vivo experiment with medaka, whose bone lacks osteocytes. Medaka fin rays with bone fractures were isolated after LIPUS treatment and RNA was extracted, as with the zebrafish samples. Of the ten transcription genes shown in Table 2, eight genes could measure by qPCR. Primer sets of other two genes were not be able to make because medaka reference sequences were not found in the public database.
Interestingly, Egr1, Egr2, FoxQ1, Helz, JunB, and NFATc1 did not respond to LIPUS in the medaka (Fig. 3b); these are considered osteocyte-dependent genes that are induced by LIPUS treatment. The expression of other genes, however, such as Snai2 and Nr4a1, was altered by LIPUS in the medaka, similarly to zebrafish (Fig. 3b). Since MLO-Y4 cells are osteocyte-like, in other words, late mature osteoblasts, MLO-Y4 cells seem to partially retain an osteoblastic signature 34 . Therefore, Snai2 and Nr4a1 may react to LIPUS not only via osteocytes but also via osteoblasts.
These in vivo results show that osteocytes contribute to LIPUS-promoted fracture repair via transcriptional regulation of Egr1, Egr2, FoxQ1, Helz, JunB, and NFATc1. This can induce cellular metabolism and survival mediated by the target genes or trigger various signaling pathways, such as the TGF-β-and Wnt-dependent pathways that are required for fracture healing.

Discussion
In this study, we demonstrated that zebrafish which have osteocytic bone, is more receptive to gain LIPUSinduced promotion of fracture healing than in medaka, which have anosteocytic bone. LIPUS stimulation altered the gene expression profile in osteocytes, and some of the transcription genes whose expression was affected appear to be important for exerting the role of osteocytes in modifying the function of osteogenic and inflammatory cells that are involved in the fracture healing process (Fig. 4).
Ample in vitro studies have shown that osteoblasts and chondrocytes rather than osteocytes, respond to LIPUS, and directly activating their functions and differentiation is important for LIPUS-dependent acceleration of fracture healing 35,36 . Indeed, osteoblasts and chondrocytes directly generate new bone, and it has been proposed that they are specific key players in the endochondral ossification that is part of the fracture healing process. Whereas, osteocytes are terminally differentiated cells and do not exert a direct effect on osteogenesis, so how osteocytes contribute to the effect of LIPUS on fracture healing is unclear, despite their mechano-sensitivity. In vivo studies have demonstrated that LIPUS accelerates all stages of the fracture repair process (inflammation, To clarify the roles of osteoblasts and osteocytes in the response to LIPUS, in vivoexperiments with genetic animal models lacking each of these cell types would be useful. However, since osteocytes are derived from osteoblasts, it is difficult to conduct long-term targeted ablation of one of these cell types but not the other in animal models. Thus, in vivo studies still have not defined the main target cells contributing to the efficacy of LIPUS. In this report, we utilized medaka as a natural osteocyte knock-out model to elucidate the effect of osteocytes in LIPUS-induced fracture repair. Zebrafish, which have osteocyte-rich bone, were utilized as a comparative model. Because the size, ecology, and swimming mode of both fish species are similar, they have previously been used in a comparative study to examine the role of osteocytes in bone modeling 2 . Comparisons of these species are partially suitable for bone research focusing on osteocytes, but the differences between them due to evolutionary distance should not be ignored. It is estimated from genomic comparisons that medaka and zebrafish separated from their last common ancestor approximately 110 million years ago 40,41 . This evolutionary distance is reflected in their biological functioning. For example, regeneration of the heart following injury is observed in zebrafish and not in medaka 42 , and the response to the estrogen 17β-estradiol occurs faster in medaka than in zebrafish 43 . Although fracture healing processes seem to be similar in the two species, their biological differences could affect experiments on fracture healing with LIPUS stimulation. To obtaining the more precise findings by in vivo using fish, the present study, utilizing only medaka and zebrafish, is not sufficient. Other than zebrafish, goldfish (Carassius auratus), carp (Cyprinus carpio), catfish (Mystus macropterus), and salmon (Salmo salar and performed a comparative study with the common carp (with cellular bone) and tilapia (with acellular bone) to investigate the mechanical properties of the two types of bone tissue 46 . Further studies using other fish species are needed, to reveal more clearly the role of osteocytes in LIPUS efficacy. We focused on Egr1, Egr2, FoxQ1, Helz, JunB, and NFATc1as osteocyte-specific LIPUS-sensitive genes. Egr1 and 2, and especially Egr1, have been reported to encode transcription factors that are induced by mechanical stimulation in cells such as endothelial cells 47 , tendon cells 48 , and myocytes 49 . In osteoblastic cells, 15 min of gravity loading induced an increase in the expression of Egr1 50 , and substrate stretching (Flexcell) also upregulated Egr1 expression 51 . In osteocytic cells, Egr1 expression was upregulated by a lack of the phosphate-regulating gene PHEX 52 and activation of parathyroid hormone (PTH) signaling 53 . However, we found no studies indicating an interaction between mechanical stress, such as LIPUS treatment, and Egr1 expression in osteocytes. Also, we could not find any research articles related to the role of Helz in bone cells including osteocytes and in mechanical stimuli. It has been reported that JunB expression is altered by mechanical stretching or loading in murine osteoblastic cells 54 and in rat chondrocytes 55 , but not established in osteocytes. FoxQ1 expression regulates the osteogenic differentiation of mouse bone MSCs 56 , but its role in osteocytes and osteoblasts is still unknown. NFATc1 was highly expressed in the mouse parietal bones and MLO-Y4 cells and nuclear translocation of NFATc1 was induced by the compressive force in osteocytes embedded in murine parietal bones 57 . In the present study, we found several novel candidate transcription genes that are mechano-sensitive and contribute to the promotion of fracture healing by LIPUS, most likely in an osteocyte-dependent manner.
To predict the functional meaning of the LIPUS-regulated expression of Egrs, FoxQ1, Helz, JunB and NFATc1in osteocytes, we performed target gene analysis and focused on types of molecular signaling that may increase the potency of osteocytes in repairing bone fractures in cooperation with neighboring cells. The target genes of these transcription factors trigger TGF-β signaling, and PGE2 synthesis. A previous global gene expression analysis using exon arrays in rat bone following mechanical loading identified loading-induced genes and performed clustering analysis 58 . In that study, the forelimbs of rats were loaded for three minutes per day, and global gene expression changes were evaluated over a time course of four hours to 32 days. Early-response genes that were upregulated within 12 h included JunB but not Egrs or FoxQ1. Interestingly, many gene groups known to be essential for bone formation were identified within clusters related to matrix formation, Wnt/β-catenin www.nature.com/scientificreports/ signaling, and TGF-β signaling. This result was partially consistent with our RNA-seq analysis of MLO-Y4 cells treated with LIPUS. Sclerostin (SOST) expressed in osteocytes is a key factor in regulating mechanical stimulation-induced bone change through TGF-β signaling and PGE2 synthesis in mammalian. Osteocytes secrete the protein SOST, encoded by the SOST gene, which negatively regulates bone mass and osteoblast differentiation by inhibiting Wnt/β-catenin signaling 59 . It has been reported that SOST is induced by mechanical loading via a TGF-βdependent mechanism 60 , and osteocyte-intrinsic TGF-β elevates the expression of SOST genes 61 . Indeed, in osteocyte-specific TGF-β receptor II-deficient mice, osteocyte-intrinsic TGF-β signaling maintained bone quality and fracture resistance through perilacunar-canalicular re-modeling, which involved the activity of osteocytes, osteoblasts, and osteoclasts 61 . Osteocytes also release PGE2, which may stimulate osteoblastic activity via control of the Wnt/β-catenin pathway 62 , probably because PGE2 represses SOST expression 63 . Interestingly, SOST and PGE2, which control osteoblastic activity, are modulated by mechanical loading in osteocytes and enhance bone formation and remodeling as a consequence of the mechanical response of osteocytes 64,65 . In addition, TGF-β and Cox2/PGE2 are important molecules for causing inflammatory cells to commit to the early phases of fracture repair 3 . However, previous reports showed that SOST expression was not observed in osteocytes in both zebrafish and medaka 15 and we confirmed that by qPCR with fin ray sample. MLO-Y4 cells is known to express a very low level of SOST and our RNA-seq data defined that. Therefore, although it is difficult that our results are extrapolated to mammalian system simply, our findings could be cue for investigation of SOST-independent alternative system for LIPUS-mediated fracture healing through TGF-β signaling and PGE2 synthesis.
Despite the fact that the effect of "active" bone-forming osteoblasts and chondrocytes on LIPUS-induced fracture repair is well reported, the role of "static" osteocytes, which control bone generation indirectly, is still not fully understood. In addition, the synergistic activity of the many cell types triggered by LIPUS during fracture healing in situations mimicking physiological conditions remains to be elucidated. The present in vivo and in vitro studies show that osteocytes are LIPUS-sensitive cells, and the reactions they assist as a result of early-response transcription genes that affect bone-forming and inflammatory cells seem to be necessary for the promotion of fracture healing. However, further experiments are required to investigate whether the LIPUS-targeted signals and molecules in osteocytes do actually affect the functions of other collaborating cells, and to determine what kind of molecular and functional changes occur in other cells to influence biological or physiological activity during fracture repair. Because the molecular response of osteocytes to LIPUS can change in the long term as fracture healing progresses, it is also important to examine changes in the effect of LIPUS on the control of osteocytes over other cells and on the harmonization of the functions of various cells in each phase of fracture healing.

Methods
Ethics statement. There was no need of obtaining permissions for conducting experiments using fish in this study. In the current laws and guidelines of Japan relating to animal experiments of fish, experiments using fish are allowed without any ethical approvals from any authorities. All the experiments presented here were conducted in accordance with the guidelines established by the Hokkaido University Animal Experiment Committees and the ARRIVE guidelines (http:// www. nc3rs. org. uk/ page. asp? id= 1357).
Fish and creating bone fracture model. Sexually mature adult (5 to 12 months of age) male and female adult zebrafish (Danio rerio) and medaka (Oryzias latipes) were utilized. In vivo experiments shown in Fig. 1, wild type zebrafish were obtained from RIKEN Brain Science Institute (Saitama, Japan) and wild type medaka were obtained from the Japan National BioResource Project (NBRP) Medaka (Okazaki, Japan). For RNA isolation, both zebrafish and medaka were obtained from a local aquarium store (Homac, Sapporo, Japan). The fish were kept at 28 °C under a 14 h light/10 h dark cycle and automatically fed TetraMin Tropical Flakes (Tetra, Melle, Germany) twice a day. To create the bone fracture model, the fish were anesthetized with tricaine (Sigma-Aldrich, St. Louis, MO, USA) and their tail bones were fractured, using a scalpel, under a microscope. To evaluate the time point that the fracture was completely healed, analysts who did not know the data source were blindly measured it.
Bone staining. Fishes were sacrificed with an overdose of Tricaine solution (4 g/L in water). The fishes' caudal fin rays were fixed overnight at 4 °C with 4% paraformaldehyde in PBS for Alizarin red and Alcian blue staining. After fixation, the fin rays were stained with Alcian blue solution (70% ethanol and 30% acetic acid containing 0.1% Alcian blue) at room temperature overnight. The fins were treated with ethanol and washed in water. Alizarin red solution (4% Alizarin red, 0.5% potassium hydroxide) was added and the fins were incubated overnight at room temperature. Stained samples were stored in 80% glycerol for imaging..
In vivo LIPUS stimulation. One day after the bone fracture, the fish were placed in six-well plates with water (one fish per well). LIPUS (Teijin Pharma, Tokyo, Japan) was generated using an array of six PZT-4 (leadzirconate titanate) transducers (2.5 cm diameter) fixed with a locking device. The plate containing the fish was placed on the array and the locking device was immersed in a water tank. The LIPUS signal was set to 1.5 MHz, 200 ms burst width sine waves at 1.0 kHz, which were delivered at an intensity of 30 mW/cm 2 . The fish were exposed to LIPUS for 20 min daily for up to four weeks or until the fracture was completely healed. The bone fracture sites were observed using an optical microscope (Nikon, Tokyo, Japan) every day.
Total RNA extraction from fish fin ray. The  www.nature.com/scientificreports/ USA) and total RNA was extracted following the instructions for the RNeasy Mini RNA isolation kit (Qiagen, Hilden, Germany).
Cell culture and LIPUS stimulation. MLO-Y4 cells were kindly provided by Dr. Lynda F. Bonewald from the Department of Oral Biology at the Kansas City School of Dentistry, University of Missouri, Kansas City, MO, USA. Cells were maintained in α-modified essential medium (α-MEM) containing 10% fetal bovine serum (FBS) at 37 °C in a humidified atmosphere of 95% air and 5% CO 2 . Prior to LIPUS stimulation, the cells were seeded on collagen-coated six-well dishes and precultured overnight. The cell culture plates were loaded onto the LIPUS system array and fixed with the locking device. The locking device was immersed in a water tank and LIPUS stimulation was provided for 20 min (1.5 MHz, pulsed-wave mode intensity of 30 mW/cm 2 ). LIPUS system has six ultrasound transducers for use in a six-well cell culture plate. To generate LIPUS treated and non-treated (control) cells in one plate, three ultrasound transducers were plugged and the other three were unplugged to shut down the LIPUS stimulation. The cells were harvested after the LIPUS treatment and total RNA was extracted using a spin column kit (Zymo Research, Orange, CA, USA).
Quantitative RT-PCR. Total RNA was extracted from the fish fin rays and MLO-Y4 cells. RNA (0.5-1 µg) was reverse transcribed using a high-capacity complementary DNA (cDNA) reverse transcription kit (Applied Biosystems, Forster City, CA, USA). qRT-PCR assays were run and quantified in the ABI STEP one real-time PCR system using SYBR green PCR Master Mix (Qiagen). Relative mRNA expression was determined using the ΔCt method, and the values were normalized to the expression of β-actin. The primer sets used for qPCR are shown in Table 3.

RNA-sequencing and analysis.
Total RNA was isolated from the MLO-Y4 cells 30 min after LIPUS stimulation, using a spin column kit (Zymo Research). Total RNA integrity and purity were assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA). The Illumina TruSeq Stranded mRNA protocol was used for the preparation of RNA-seq libraries and sequenced on a HiSeq 2500 machine (Illumina, San Diego, CA, USA) as paired-end, 100-base pair reads and 20 million reads per sample were obtained. Pairedend sequencing reads from the HiSeq 2500 were trimmed using Trimmomatic v0.36 and trimmed reads were aligned to the GRCm38 genome assembly using TopHat v2.1.0 and bowtie v2.2.6.0. Gene expression intensity was normalized to Fragments per kilobase of exon per million reads mapped (FPKM) which was calculated using Cufflinks (version 2.2.1) with a transcriptome reference (Ensembl Mouse Transcript). Signature genes for each group were identified using an adjusted q-value of < 0.05. Significant difference was defined by q-value. q is an false discovery rate (FDR)-adjusted enrichment p value, and q < 0.05 (i.e., − log 10(q) > 1.3) was defined as significant. A principal component analysis (PCA) was performed using PCAGO, an interactive web service (https:// pcago. bioinf. uni-jena. de) to explore the variation between samples. We used a subset of the top 500 genes with most variable expression. The Functional Annotation Tool of the Database for Annotation, Visualization, and Integrated Discovery (DAVID) v6.8 was used to characterize the gene annotation enrichment analysis.
Target gene analysis of transcription factor. To identify the putative target genes of the transcription factors, the mouse ChIP-seq database in ChIP-Atlas (http:// chip-atlas. org/) was used. If a transcription factor's ChIP-seq peaks settled within the gene's promoter region (5 kb around the transcription start site), the gene was identified as a target gene 66 . Statistics. All statistics were calculated using Microsoft Excel and GraphPad Prism. All data are presented as a mean ± SEM and comparisons were made using Student's t-test. The SEM was chosen to compare the population mean of the two groups. Significance levels of results were defined as follows: * p < 0.05, ** p < 0.01, and *** p < 0.001. All experiments presented were performed in two to three independent experiments, except for the RNA-seq study. The experiments were not randomized and sample size was not predetermined.

Data availability
Data supporting the conclusions are available from the corresponding author upon reasonable request. RNA-seq data have been deposited to Gene Expression Omnibus (GEO) repository with accession number GSE162674. www.nature.com/scientificreports/