Continuous dynamic identification of key genes and molecular signaling pathways of periosteum in guided bone self-generation in swine model

Background Guided bone self-generation with periosteum-preserved has successfully regenerated mandibular, temporomandibular and interphalangeal joint. The aim of this study was to investigate the dynamic changes of gene expression of periosteum which was involved in the guided bone self-generation. Methods Rib defects of critical size were created in mature swine with periosteum-preserved. The periosteum was sutured into a sealed sheath that closed the bone defect. The periosteum of trauma and control sites were harvested at postoperative 9 time points, and total RNA was extracted. Microarray analysis was conducted to identify the differences in the transcriptome of different time points between two groups. Results The differentially expressed genes (DEGs) between control and trauma group were different at postoperative different time points. The dynamic changes of the number of DEGs fluctuated a lot. There were 3 volatility peaks, and we chose 3 time points of DEG number peak (1 week, 5 weeks and 6 months) to study the functions of DEGs. Oxidoreductase activity, oxidation–reduction process and mitochondrion are the most enriched terms of Go analysis. The major signaling pathways of DEGs enrichment include oxidative phosphorylation, PI3K-Akt signaling pathway, osteoclast differentiation pathway and Wnt signaling. Conclusions The oxidoreductase reaction was activated during this bone regeneration process. The oxidative phosphorylation, PI3K-Akt signaling pathway, osteoclast differentiation pathway and Wnt signaling may play important roles in the guided bone self-generation with periosteum-preserved. This study can provide a reference for how to improve the application of this concept of bone regeneration. Supplementary Information The online version contains supplementary material available at 10.1186/s13018-023-03524-y.


Introduction
Periosteum, a highly vascularized connective tissue, plays a key role in the growth, development and regeneration of bone [1]. It covers the external surface of bones except for the sites of articulation surface and muscle attachment. The inner cambium layer of periosteum contains skeletal stem cell population [2]. The skeletal stem cell population roles importantly in the fracture healing and bone development, which could differentiate into chondrogenic and osteogenic lineages [3,4].
Bone formations include intramembranous bone formation and endochondral bone formation [5]. During the process, the periosteum lays down intramembranous woven bone. The tissue damage closer to the fracture site could block the blood supply and result in the formation of cartilaginous masses [6]. Then the endochondral bone formation was observed. The chondrocytes were reported to derive predominantly from the periosteum. Therefore, periosteum is key factor in the two kinds of bone formation.
With high regenerative capacity, periosteum has high potential of the utilization for bone regeneration [7]. In clinical studies, it has been reported that spontaneous bone regeneration could be observed to re-establish the bony continuity for bone tumor segmentectomy with periosteum preservation [8,9]. This phenomenon provides enlightening resolution for endogenous bone tissue engineering for bone defects. We have created animal models of bone defects to study the periosteum and bone regeneration, and reconstructed new bone tissues which were histologically, anatomically and functionally similar to normal bone tissue [10][11][12]. The guided production of autologous bone without exogenous cells, scaffolds and growth factors is a novel concept deserving exploring in further clinical application like self-generation joint and cartilage. Therefore, a better understanding of molecular mechanisms is important for accelerating the bone regeneration process.
In this study, we used microarray to investigate the dynamic changes of gene expression of the periosteum involved in the bone regeneration for bone defect at different time points. The study may help us to have a better understanding of the molecular regulation mechanisms that periosteum is involved in during the bone regeneration of periosteum-preserved bone defect.

Materials and methods
This study was approved by the Animal Research Ethics Committee of the authors' institution and carried out in accordance with the approved guidelines (SH9H-2019-A645-1).

Subjects and surgical procedure
Thirty-six healthy skeletally mature female swine were used in this study. The age was ranging from 4 to 6 months. Under general anesthesia, the left side of ribs was chosen as experimental group and the right side as control group. Figure 1a shows the schematic diagram of this procedure. At the surface of the fourth rib, the chondro-osseous transition point was touched and an incision beginning from the transition point and extending for 3 cm toward the spine was created. The third, fourth and fifth ribs were exposed. The periosteum was incised at the middle of the ribs, and a segment of rib was removed, resulting in a 3 cm defect (Fig. 1b). The periosteum was sutured with 6-0 PDS sutures, creating a closed space and forming a sealed sheath (Fig. 1c). Then the wound was closed.

Next-generation sequencing (NGS)
On postoperative 1 day, 3 days, 1 week, 2 weeks, 1 month, 5 weeks, 3 months, 6 months and 7 months, the swine were killed, and the periostea were harvested (Fig. 1d). The periostea were placed into cryotube for RNA-seq. The periostea and formed bone were stained to observe the osteogenesis. The process includes total RNA extraction, quality detection, mRNA purification, mRNA fragmentation, cDNA synthesis, PCR enrichment of library fragments, library quality inspection and on-machine sequencing.
The raw data were filtered, and the filtered high-quality sequence was compared to the reference genome. The expression level of each gene was calculated. Further analysis of expression differences, enrichment analysis and cluster analysis were performed. We perform FPKM density distribution analysis, saturation analysis, sample correlation test and principal components analysis (PCA) analysis. DESeq software package in the R language was used to analyze the PCA principal components based on the expression level.
The R language ggplots2 software package was used to draw volcano maps of DEGs. R language Pheatmap software package was used to perform two-way clustering analysis on the union of different genes and samples of all comparison groups, clustering according to the expression level of the same gene.
TopGO was used for GO enrichment analysis. GO classification was carried out according to molecular function (MF), biological process (BP) and cell component (CC). The enrichment degree was measured by Rich factor and FDR value. The top 20 GO term entries with the smallest FDR value were displayed. According to the KEGG enrichment analysis results of DEGs, the degree of enrichment was measured by Rich factor and FDR value. The top 20 KEGG pathways with the smallest FDR value were also displayed.

Statistical analysis
Data are expressed as the mean ± S.D. Fold change and P value were calculated with t-test to identify DEGs. The threshold set for differentially expressed genes was the |log2FoldChange|> 1 and the P-value < 0.05. Functional analysis of these DEGs was conducted with associated database mentioned in the "Materials and methods" section.

Results
All the surgeries were completed successfully. The RNA quality inspection results suggested the RNA quality of all samples was good. Therefore, the RNA samples were reliable for further transcriptomic analysis. In Additional file 1: Fig. S1, the volcano maps show the DEGs between two groups were different at each time points. The dynamic changes of the number of DEGs are shown in Fig. 2, in which we can intuitively see that the number of up-regulated and down-regulated genes fluctuated a lot. There were 3 volatility peaks, and we chose 3 time points of the volatility peaks (1 week, 5 weeks and 6 months) to study the functions of DEGs.
Cluster analysis was used to determine the expression patterns of differentially expressed genes under different experimental conditions. Genes with high expression correlation between samples are grouped into a class; usually these genes have actual correlation in some biological process, or a metabolic or signaling pathway. Therefore, through expression clustering, we can find the unknown biological relationship between genes. The results of cluster analysis are shown in Additional file 2: Fig. S2 and show biological connections between genes. Genes are represented horizontally, with each column representing one sample, with red representing highexpressed genes and green representing low-expressed genes. From the figure, we can see that the expression levels of these genes in different samples and the expression patterns of different genes in the same sample are clustered, indicating that there are actual links in the biological processes of these genes, or in a certain metabolic or signaling pathway.
The top 10 GO term entries classified into MF, BP and CC are shown in Additional file 3: Fig. S3. At postoperative 1 week, the most top terms of three categories were cytoplasm, oxidoreductase activity and oxidation-reduction process (Additional file 3: Fig. S3a). At postoperative 5 weeks, the most top terms were mitochondrion, oxidoreductase activity and oxidation-reduction process  (Additional file 3: Fig. S3b). At postoperative 6 months, the most top terms were extracellular region, glycosaminoglycan binding and extracellular matrix organization (Additional file 3: Fig. S3c). The top 20 GO term entries with the smallest FDR value are shown in Fig. 3. The most significant enrichment of the top 3 GO terms according to FDR value at postoperative 1 week include oxidationreduction process, cytoplasm and myofibril ( Table 1). The most significant enrichment of the top 3 GO term entries according to FDR value at postoperative 5 weeks include mitochondrion, cytoplasm and oxidation-reduction process ( Table 2). The most significant enrichment of the top 3 GO term entries according to FDR value at postoperative 6 months include extracellular region, extracellular space and extracellular matrix ( Table 3).
As to the results of KEGG analysis, the top 20 pathways with the smallest FDR value are shown in Fig. 4. The top 3 pathways of KEGG enrichment analysis of postoperative 1 week include Parkinson's disease, oxidative phosphorylation and thermogenesis. At 5 weeks after the surgery, the top 3 pathways include Parkinson's disease, citrate cycle (TCA cycle) and non-alcoholic fatty liver     The PI3K-Akt signaling pathway, which was ranked the third top KEGG enrichment of 6 months, was also ranked 91th at postoperative 1 week and 123th at postoperative 5 weeks. The PI3K-Akt signaling pathway is shown in Additional file 4: Fig. S4. The associated DEGs at postoperative 6 months include EGF, FN1, ANGPT4, SPP1, COMP, TLR4, COL1A2, CSF1R, ITGA11, CREB3L1, NGFR, COL6A1, TLR2, PCK1, THBS3, COL6A2, FGF10 and IGF1.
The osteoclast differentiation pathway, which was ranked the 14th top KEGG enrichment of 6 months, was also ranked 212th and 149th at postoperative 1 week and 5 weeks. The osteoclast differentiation pathway is shown in Additional file 5: Fig. S5. The associated DEGs at postoperative 6 months include FOS, CSF1R, NFATC2, PPARG, TYROBP, NCF2 and NCF4.
The Wnt signaling pathway, which was ranked the 15th top KEGG enrichment of 6 months, was also ranked 129th and 128th at postoperative 1 week and 5 weeks. The Wnt signaling pathway is shown in Additional file 6: Fig. S6. The associated DEGs at postoperative 6 months include WNT2B, FZD7, NFATC4, WNT16, NFATC2, SFRP2, ROR2, BAMBI and SFRP1.
The results of periosteal and bone staining show that at postoperative 3 days, it is mainly a period of hematoma formation and a small amount of cartilage formation (Fig. 6A). At postoperative 1 week, the hematoma gradually disappeared and more cartilage was formed (Fig. 6B). At postoperative 5 weeks, the periosteal sheath was filled with cartilage tissue and new trabeculae gradually replaced the cartilage tissue (Fig. 6C). At postoperative 6 months, the periosteal sheath was already filled with new bone tissue that had been shaped (Fig. 6D).

Discussion
Current treatments for large bone defect mainly include autologous, allogeneic and alloplastic grafts, and are very limited due to a lack of donor sites, potential immunogenic response [13,14]. Guided self-generation has promising potential of clinical translation, since it need not exogenous cells, growth factors and scaffolds, which are necessary in traditional tissue engineered techniques [11,12]. Periosteum is considered to be crucial in the procedure of guided self-generation, but little is known about the molecular regulation mechanisms. This study focuses on the molecular regulation mechanisms that periosteum is involved in during the bone regeneration. The current study outlines the dynamic changes of transcription associated with periosteum in guided bone self-generation. The dynamic changes of gene expression were firstly reported in so detailed time points. The dynamic changes of the number of DEGs fluctuated very much. From the dynamically changing volatility graph of DEGs, we could find 3 volatility peaks and we chose 3 time points of DEG number peaks (1 week, 5 weeks and 6 months) to study the function of DEGs.
The cytoplasm was one of the top terms. Special chemical and physical characteristics of the extracellular matrix (ECM) could pass on the stimulus from outside the cell to the cytoplasm and then results in osteocyte cytoskeletal changes, ECM remodeling and altering bone [15]. Oxidoreductase activity was the top term both at postoperative 1 week and 5 weeks. It was reported oxidoreductase activity was inversely correlated with collagen synthesis in osteoblasts [16]. Oxidation-reduction process was also the top term. Oxidation-reduction process was significantly dysregulated following the formation of heterotopic ossification [17]. Genes in the oxidation-reduction process category were down-regulated in the healing bones between titanium orthodontic miniimplants site bone and sites of surgical defects [18]. The top terms suggest oxidoreductase reaction has an important role in bone regeneration. Mitochondrion was the top term of postoperative 5 weeks, and mitochondrion plays important role in generating reactive oxygen species (ROS), thus affecting bone remodeling [19,20]. Extracellular region and extracellular matrix organization were the top terms of postoperative 6 months. Extracellular matrix such as enzymes and ions is often involved in bone metabolism. Extracellular Zn 2+ and its transporters participate in regulating physiological or pathological processes [20,21]. Nucleotide-metabolizing ectoenzymes has been reported to be associated with bone modeling [22][23][24]. However, the important role of these enzymes in bone metabolism has not been fully understood yet.
Oxidative phosphorylation is one of the top pathways both at postoperative 1 week and 5 weeks. It is consistent with the results of GO analysis. Studies have found there existed a metabolic shift for stem cell differentiation [25][26][27]. Glycolytic metabolism is the main metabolism in stem cells for their self-renewal, and it has a benefit for cell proliferation [28]. There is an increased cellular energy demand when stem cells commit to differentiate. During this procedure, glycolytic metabolism could be shifted to mitochondrial oxidative phosphorylation (OXPHOS). Ram et al. investigated mitochondrial biogenesis during differentiation of periosteum-derived mesenchymal stem cells (POMSC), which express the typical mesenchymal stem cells (MSC) surface markers and have multipotency [28,29]. Mitochondrial proteins including OXPHOS and VDAC, and mitochondrial DNA (mtDNA) contents were all increased during the differentiation of POMSC. However, glycolytic metabolism is decreased during osteogenic differentiation. Furthermore, osteogenic differentiation of POMSC could be prevented by reducing mtDNA content with ethidium bromide treatments. Therefore, OXPHOS may play a key role in the differentiation of POMSC. POMSC differentiation and regenerative medicine could be regulated by modulation of OXPHOS with pharmaceutical [30,31].
The PI3K-Akt signaling pathway was ranked the third top KEGG enrichment of 6 months. It was reported that FN1 could activate PI3K-Akt signaling pathway to modulate the biology of cells [32]. The up-regulation of PI3K-Akt signaling pathway has been reported to stimulate fracture healing [33,34]. Heli et al. [35] found that overexpression of FN1 could make contributions to fracture healing through the activation of TGF-β/PI3K-Akt signaling pathway. In this current study, FN1 was one of the up-regulated DEGs. We could infer that the FN1 may activate the PI3K-Akt signaling pathway to affect bone self-generation. Concetta et al. reported that periosteumderived progenitor cells could be induced toward osteoblastic differentiation under IGF1 or IGF1/VEGF stimuli with the activation of PI3K-Akt signaling pathway, which plays important role in osteogenesis processes [36]. The IGF1 was also differentially up-regulated in our study, which suggests that IGF1 may affect the activation of PI3K-Akt signaling pathway, thus regulating the bone regeneration.
The osteoclast differentiation pathway was ranked the 14th top KEGG enrichment of 6 months. The function coordination of skeletal cells which include bone-forming osteoblasts and bone-resorbing osteoclasts maintains the homeostasis of bone [37]. Many factors affect the differentiation of osteoclasts, including receptor activator of nuclear factor (NF)-κB ligand (RANKL), tumor necrosis factor (TNF) family cytokine and macrophage colonystimulating factor (M-CSF) [38]. The c-Fos is an essential molecule for osteoclastogenesis and was up-regulated. We speculate that the activation of this pathway is mainly associated with callus shaping at different time points. The Wnt signaling pathway was top pathway of 3 time points. It was reported that Wnt signaling was involved in different cellular activities during cell differentiation, organogenesis and early embryonic development [39]. Fracture repair was considered to be a recapitulation of embryonic development, and members of the Wnt signaling pathway were activated [40]. Nan et al. found that Wnt signaling pathway was activated during bone regeneration [40]. The Wnt signaling pathway may play an important role in the regenerative process.
Several strategies, including gene therapy and tissue engineering together with mesenchymal stem cells (MSC), have been proposed to promote the healing of the musculoskeletal tissue. Moreover, a recent technology has revolutionized gene editing: Clustering regulatory interval short palindromic repeats (CRISPR) features simple target design, affordable, versatile and efficient, but requires more research to be the preferred platform for genome editing. Predictive genomics DNA analysis can understand which genetic advantages (if any) can be exploited and why specific rehabilitation programs are more effective in some people than others [41,42]. Therefore, a better understanding of the genetic impact on musculoskeletal system function and disease healing is needed to plan and develop patient-specific management strategies. Currently, while some results are promising, all biological interventions are experimental and the cost/effectiveness has not been proven. In addition, the short follow-up time of most studies questioned the durability of treatment [42]. In this study, autologous periosteum was used to guide bone regeneration in vivo, and the regenerated bone was used for precise repair of the body. This technology has a promising clinical application prospect. The related differential genes and signaling pathways identified in this study can provide a rich theoretical basis for later gene and molecular intervention.
Though this in vivo study provides a better understanding in the molecular mechanisms involved in guided self-generation, it also has limitations. First, micro-CT can be conducted to evaluate the conditions of bone regeneration at different time points. The association of changes of molecules and bone regeneration can be analyzed to provide more precise explanations of the mechanisms. Besides, different parts of the regeneration may be at different stage of healing process.