Differential expression and effect analysis of lncRNA-mRNA in congenital pseudarthrosis of the tibia

Background: To analyze the lncRNA-mRNA differential expression and co-expression network of periosteal stem cells (PSCs) from congenital pseudarthrosis of the tibia (CPT) and normal patients, and to explore the role of key lncRNAs.

Methods: Differentially expressed lncRNAs and mRNAs in PSCs were obtained by sequencing, and biological functions of differentially expressed mRNAs were detected by gene ontology (GO), Kyoto encyclopedia of genes and genomes (KEGG) pathway and protein -protein interaction (PPI) analysis. The co-expression network of lncRNA-mRNA was constructed by correlation analysis of differentially expressed lncRNAs and mRNAs, and the key lncRNAs were screened according to the connectivity degree. After that, the cis-regulated target genes of differential expressed lncRNAs and mRNAs were predicted.

Results: A total of 194 differentially expressed lncRNAs were identified, including 73 upregulated and 121 downregulated genes. A total of 822 differentially expressed mRNAs were identified, including 311 upregulated and 511 downregulated genes. GO, KEGG and PPI enrichment analysis showed that the regulatory function of differentially expressed mRNAs were mainly gathered in skeletal system development and tissue morphogenesis. The co-expression network with 226 nodes and 3,390 edges was constructed based on correlation analysis. A total of 10 key lncRNAs, including FAM227B, POM121L9P, AF165147 and AC103702, were screened according to connectivity degree. Prediction of target genes indicated that FAM227B-FGF7 and AC103702-HOXB4/5/6 may play an important role in the pathogenesis of CPT.

Conclusion: A total of 10 key lncRNAs, including FAM227B, POM121L9P, AF165147, and AC103702, occupy the core position in the co-expression network, suggesting that these lncRNAs and their target genes may play an important role in the pathogenesis of CPT.

1 Introduction

Congenital pseudarthrosis of the tibia (CPT) is regarded as a rare, refractory and long-term congenital disease in children. It is reported that the incidence rate of CPT ranges from about 1: 140,000 to 1: 250,000 (Hefti et al., 2000). Since Paget reported the disease in 1891, its pathogenesis remains to be unclear (BOYDHarold, 1982; Pannier, 2011). Currently, surgical treatment is the main treatment for CPT patients in clinic. However, according to the literatures, the risk of complications such as bone non-union, re-fracture, limb length difference and ankle valgus after CPT can reach 50% (Li et al., 2022), which results in reoperation, as well as even severe consequences such as disability and amputation (Zhou et al., 2022) and brings huge psychological and economic burden to patients and their families. Therefore, it is the key to seek more effective treatment through deeply studying the pathogenesis of CPT.

Previous studies have found that tibial periosteum lesion is an important pathogenic factor causing the decline of local bone osteogenesis, which results in the formation of CPT. Therefore, it has gradually become a hot spot to explore the biological characteristics of periosteum lesion in CPT. Periosteum is a connective tissue membrane which covers the bone surface and possesses rich micro-vascular tissue. It has good osteogenic ability and plays an important role in bone growth, development and injury repair (Filion et al., 2017). Periosteum is mainly composed of three layers: fibrous layer, undifferentiated layer and germinal layer. Germinal layer is composed of periosteal stem cells (PSCs), mesenchymal stem cells, fibroblasts and other cells (Augustin et al., 2007). Recent studies have confirmed that periosteal stem cells in germinal layer have excellent osteogenic potential, and they can differentiate into osteoblasts under specific conditions, which are important functional cells involved in the process of bone remodeling and repair (Doherty et al., 2020; Matthews Brya et al., 2021). The osteogenic differentiation ability of PSCs is affected by many factors, such as local micro-environment (various hormones, cytokines, transcription factors) (Granchi et al., 2010), but the specific ways to affect the osteogenic differentiation ability of PSCs still remain to be unknown.

In recent years, long-chain non-coding RNA (lncRNA) has been the hot spot of gene expression regulation. LncRNA is a non-coding RNA composed of more than 200 nucleotides (Miao and Lee, 2015). Initially lncRNA was considered as noise in the transcription process, without biological function (Iyer Matthew et al., 2015). However, more and more studies have shown that lncRNA is closely related to many physiological activities and pathological processes. Some studies reported that lncRNA participates in the process of regulating many common diseases, such as various cancers, cardiovascular diseases, metabolic diseases (Pearson et al., 2016; Mila and Alessia, 2019; Macdonald and Mann, 2020). With the development of technology, more and more scholars found that lncRNA is also related to some rare diseases, such as Prader-Willi syndrome (Wu et al., 2016), but there are still few studies on CPT.

In order to better understand the pathogenesis of CPT and explore the clues of the key pathogenic genes, we used the lncRNA microarray technology to detect the differences in lncRNA and messenger RNA (mRNA) expression profiles of PSCs between normal people and CPT patients. In addition, the function and interaction of lncRNAs and mRNAs were comprehensively analyzed, aiming to provide experimental basis for finding new therapeutic targets for CPT.

2 Materials and methods

2.1 Sample acquisition

The samples of this study came from patients who were admitted into the department of Orthopaedics of The First Affiliated Hospital of Zhejiang University, School of Medicine and Hunan Children’s Hospital from 2018 to 2020. Normal periosteal tissues were obtained from the operation of healthy children with tibial fracture caused by trauma in The First Affiliated Hospital of Zhejiang University, School of Medicine, and the diseased periosteum tissues were obtained from CPT patients during surgery in Hunan Children’s Hospital. A total of 10 periosteum samples were obtained and divided into two groups: control group and CPT group, with five from CPT patients (3–8 years old) and five from traumatic children (3–6 years old) separately. This study had been approved by the Ethics Committee of the First Affiliated Hospital of Zhejiang University, School of Medicine and Hunan Children’s Hospital, and informed consent had been obtained from all patients.

2.2 Isolation, purification and identification of PSCs

After the periosteal tissue sample was cut to 1 mm³, it was digested at 37°C for 12 min with collagenase A (2 mg/mL, Roche) and trypsin (2.5 mg/mL, Roche). After DMEM culture solution containing 10% FBS was added, the supernatant was removed by centrifugation. DMEM culture solution containing 10% FBS was added, and the tissue was resuspended gently for 5–10 times. The suspension was filtered by a 70-um cell sieve, and the filtrate was inoculated into a culture dish. When the cell confluence reached 90%–100%, the cells were digested by trypsin for subculture. STRO-1 was used as marker, and PSCs were purified by indirect immunomagnetic bead positive selection. DSB-X Botin STRO-1 labeled antibodies were co-incubated with PSCs at 4°C for 30 min, and the supernatant was removed after centrifugation at 1000 rpm for 8 min 80 ul of magnetic beads and PSCs were mixed for 30 min and put them in the magnetic field for 3 min to remove the supernatant. Flowcomp release buffer was added. After repeated mixing, magnetic field was placed to demagnetize the beads so as to obtain purified PSCs. Fluorescently labeled monoclonal antibodies CD31, CD34, CD44, and CD90 were added to the purified PSCs, and the mixture was mixed and left for 15 min at 4°C. After washing 3 times with PBS, 500 ul of PBS containing 0.5% FBS was added and transferred to a flow cytometer for detection.

2.3 RNA extraction library construction and sequencing

LC bio-technology Co., Ltd. (Hangzhou, China) provided microarray technical support for analyzing and evaluating lncRNA and mRNA expression profiles. The microarray sequencing process is shown in Figure 1. Total RNA was extracted using Trizol reagent (hermofisher, 15596018) following the manufacturer’s procedure. The total RNA quantity and purity were analyzed of Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, United States, 5067-1511) with RIN number >7.0. Approximately 5 ug of total RNA was used to deplete ribosomal RNA according to the manuscript of the Ribo-Zero Gold rRNA Removal Kit (Illumina, cat. MRZG12324, San Diego, United States). After removing ribosomal RNAs, the left RNAs were fragmented (NEBNext^® Magnesium RNA Fragmentation Module, cat. E6150S, United States) into short fragments using divalent cations under high temperature. Then the cleaved RNA fragments were reverse transcribed to create the cDNA by SuperScript™ II Reverse Transcriptase (Invitrogen, cat. 1896649, United States), which were next used to synthesise U-labeled second-stranded DNAs with E. coli DNA polymerase I (NEB, cat. m0209, United States), RNase H (NEB, cat. m0297, United States) and dUTP Solution (Thermo Fisher, cat. R0133, United States). An A-base was then added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contains a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single- or dual-index adapters are ligated to the fragments, and size selection (300–600 bp) was performed with AMPureXP beads. After the heat-labile UDG enzyme (NEB, cat. m0280, United States) treatment of the U-labeled second-stranded DNAs, the ligated products were amplified with PCR by the following conditions: initial denaturation at 95°C for 3 min; 8 cycles of denaturation at 98°C for 15 s, annealing at 60°C for 15 s, and extension at 72°C for 30 s; and then final extension at 72°C for 5 min. The average insert size for the final cDNA library was 300 ± 50bp. At last, we performed the 2 × 150 bp paired-end sequencing (PE150) on an Illumina Novaseq™ 6000 (LC Bio-Technology CO. Ltd. Hangzhou, China) following the vendor’s recommended protocol.

FIGURE 1

FIGURE 1. The schematic diagram of sequencing.

2.4 Bioinformatics analysis

2.4.1 Different expression analysis of mRNAs and lncRNAs

Genes/RNAs differential expression analysis was performed by DESeq2software (Smyth and McCarthy, 2010; Pertea et al., 2015) between two different groups (and by edgeR between two samples). The genes/mRNA/lncRNA with the parameter of false discovery rate (FDR) below 0.05 and absolute fold change ≥2 were considered differentially expressed genes/mRNA/lncRNA. Differentially expressed coding RNAs were then subjected to enrichment analysis of gene ontology (GO) functions and Kyoto encyclopedia of genes and genomes (KEGG) pathways. Cluster analysis was carried out and plotted by heat map, volcanic map.

2.4.2 GO enrichment analysis

GO (http://www.geneontology.org/) is an international standardized gene functional classification system which offers a dynamic-updated controlled vocabulary and a strictly defined concept to comprehensively describe properties of genes and their products in any organism. GO has three ontologies: molecular function, cellular component and biological process. The basic unit of GO is GO-term. Each GO-term belongs to a type of ontology. GO enrichment analysis provides all GO terms that significantly enriched in DEGs comparing to the genome background. Firstly, all DEGs were mapped to GO terms in the Gene Ontology database, gene numbers were calculated for every term, significantly enriched GO terms in DEGs comparing to the genome background were defined by hypergeometric test and then ranked according to the -log10 (p). The calculating formula of p-value is:

$P=1-{\displaystyle \sum }_{i=0}^{m-1}\frac{\left(\begin{array}{l}M\\ i\end{array}\right)\left(\begin{array}{l}N-M\\ n-i\end{array}\right)}{\left(\begin{array}{l}N\\ n\end{array}\right)}$

Here N is the number of all genes with GO annotation; n is the number of DEGs in N; M is the number of all genes that are annotated to the certain GO terms; m is the number of DEGs in M. GO terms meeting this condition with p < 0.05 were defined as significantly enriched GO terms in DEGs. Significant enriched GO terms were then hierarchically clustered into a tree based on Kappa-statistical similarities among their gene memberships. 0.3 Kappa score was applied as the threshold to cast the tree into term clusters. This analysis was able to recognize the main biological functions that DEGs exercise.

2.4.3 Pathway enrichment analysis

Genes usually interact with each other to play roles in certain biological functions. Pathway-based analysis helps to further understand genes biological functions. KEGG (www.genome.jp/kegg) is a public database on genome deciphering. The databases such as GENES and Pathway contain information of cellular biochemical processes, including genome sequence, metabolism, membrane transport, chemical substances and enzyme reactions. In organisms, different genes coordinate their biological functions, and pathway enrichment analysis identified significantly enriched metabolic pathways or signal transduction pathways in DEGs comparing with the whole genome background. The calculating formula is the same as that in GO analysis. Here N is the number of all genes that with KEGG annotation, n is the number of DEGs in N, M is the number of all genes annotated to specific pathways, and m is number of DEGs in M. Pathways meeting this condition with p < 0.05 were defined as significantly enriched pathways in DEGs.

2.4.4 Protein -protein interaction (PPI) network analysis

String software (http://string-db.org) was used to analyze the interaction between target genes encoded by differentially expressed lncRNA-mRNA and corresponding proteins, and the key proteins were screened out. Minimum required interaction score >0.7 was set to obtain the interactive relationship between proteins and construct the PPI network. MCODE (Molecular Complex Detection) algorithm was then applied to this network to identify neighborhoods where proteins are densely connected. GO enrichment analysis was applied to each MCODE network to extract “biological meanings” from the network component, where top three best p-value terms were retained. Cytoscape software was imported for visualization.

2.4.5 LncRNA-mRNA co-expression network

The differentially expressed lncRNA-mRNA co-expression network was constructed to explore their interaction in the of normal and CPT patients. By calculating Pearson’s correlation coefficient (r) of lncRNAs and mRNAs, the correlation of their expression levels was obtained. The significantly correlated lncRNA-mRNA pairs were screened by setting |r| > 0.9 and p < 0.01 for constructing lncRNA-mRNA co-expression network, and then Cytoscape software was imported for visualization. Connectivity degree is an index used to evaluate the importance of lncRNA-mRNA in the network. The greater degree is, the greater its regulation effect on the network is. By calculation, the top 10 differentially expressed lncRNAs were listed.

2.4.6 Prediction of lncRNA target gene

The main regulation modes of lncRNA include cis-regulation and trans-regulation (Tachiwana et al., 2020). The lncRNA cis-regulation target genes were mainly predicted according to the positional relationship. The differentially expressed lncRNA and mRNA within 100 kbp upstream and downstream of chromosome were defined as cis regulation. Differential lncRNAs and mRNAs in cis-regulation were predicted and analyzed.

2.5 Quantitative reverse transcription-polymerase chain reaction (qRT-PCR)

We selected the top five differentially expressed lncRNAs mentioned above which were verified by qRT-qPCR. The analysis was performed according to the instructions: the cells were lysed, and RNA was extracted with Trizol reagent (Invitrogen, 15596018). The EntiLink™first strand cDNA Synthesis Kit (ELK Biotechnology, EQ003) was used to synthesize the first strand cDNA. qRT-PCR was performed using an ABI Quant Studio six Real-Time PCR instrument (Life technologies). Three double Wells were made for each sample and SYBR Green PCR Super Mix kit (VAZYME, Q111-02) was used. The reaction conditions were 95°C (15 s), 60°C (60 s), 95°C (15 s), and 40 cycles. The quantitative method was 2^^−ΔΔ CT method. The primer sequences are shown in Table1.

TABLE 1

TABLE 1. The primer sequences.

3 Results

3.1 Differential lncRNA and mRNA screening

The normal and CPT periosteal cells were isolated and purified through the above steps, and the purified cell surface markers were identified by flow cytometry. The results showed that the average positive rate of the identified cell lines CD31 and CD34 was lower than 1%, and the average positive rate of CD44 and CD90 was higher than 98%, which was identified as PSCs (Figure 2).

FIGURE 2

FIGURE 2. Flow cytometry was used to detect the positive ratio of molecular markers CD31, CD34, CD44 and CD90 of purified PSCs derived from different periosteal tissues. PSCs, periosteal stem cells.

Further, we used the purified PSCs for sequencing analysis. According to the screening criteria (absolute fold change ≥2), 194 differentially expressed lncRNAs were identified by sequencing, including 73 upregulated genes and 121 downregulated genes, with 116 newly discovered lncRNAs and 78 known lncRNAs. There were 822 differentially expressed mRNAs, including 311 upregulated genes and 511 downregulated genes.

At the same time, cluster analysis was carried out, and the first 100 lncRNAs and mRNAs with the smallest p-value were plotted by heat map (Figures 3A, B). Moreover, the expression difference of lncRNA and mRNA in normal and CPT PSCs were analyzed by volcanic map (Figures 3C, D).

FIGURE 3

FIGURE 3. (A), the first 100 lncRNAs with the smallest p-value were plotted by heat map. (B), the first 100 mRNAs with the smallest p-value were plotted by heat map. (C), the expression difference of lncRNA in normal and CPT PSCs was analyzed by volcanic map. (D), the expression difference of mRNA in normal and CPT PSCs was analyzed by volcanic map. CPT, congenital pseudarthrosis of the tibia; PSCs, periosteal stem cells.

3.2 Gene function annotation and pathway enrichment analysis

In order to further analyze the function of differentially expressed mRNAs in normal and CPT PSCs, GO and KEGG pathway analysis were carried out. GO analysis showed that the biological processes of differentially expressed mRNAs was mainly focused on signal transduction, multicellular organism development and positive regulation by transcription RNA polymerase II; secondly the cell components were concentrated in cell membrane, plasma membrane and cytoplasm; thirdly, molecular biological functions were focused on protein binding, metal ion binding and DNA binding, etc (Figure 4A). Significant enriched GO terms were then clustered based on Kappa-statistical similarities among their gene memberships. The results showed that the biological functions with the most significant differences were mainly gathered in skeletal system development, tissue morphogenesis and NABA core matrisome (shown in Figures 4B, C).

FIGURE 4

FIGURE 4. (A), GO analysis revealed gene enrichment in different cellular biological processes. (B), the top 20 significant GO terms were hierarchically clustered into a tree based on Kappa-statistical similarities among their gene memberships. (C), the representative terms from the full cluster were selected and converted into a network layout. Each term was represented by a circle node, where its size is proportional to the number of genes fall under that term, and its color represent its cluster identity (nodes of the same color belong to the same cluster). Terms with a similarity score >0.3 were linked by an edge (the thickness of the edge represents the similarity score). GO, Gene ontology.

The KEGG pathways enrichment analysis showed that 59 pathways had significant difference (p < 0.05), with 30 upregulated pathways and 29 downregulated pathways. A total of 20 pathways with the most significant regulation were selected. KEGG analysis showed that the pathways with the most significant differences were mainly concentrated in pathways in cancer, PI3K-Akt signaling pathway and cytokine-cytokine receptor interaction (Figure 5A). The regulatory changes of related mRNAs involved in the pathways of cancer and PI3K-Akt signaling pathway were shown in Figures 5B, C.

FIGURE 5

FIGURE 5. (A), the top 20 significant signaling pathways were hierarchically clustered into a tree, which showed that the pathways with the most significant differences were mainly concentrated in pathways in cancer and PI3K-Akt signaling pathway. (B), the regulatory changes of related mRNAs involved in the pathways of cancer. Red represented upregulation and blue represented downregulation. (C), the regulatory changes of related mRNAs involved in PI3K-Akt signaling pathway. Red represented upregulation and blue represented downregulation.

3.3 PPI network analysis

For differentially expressed mRNAs, protein-protein interaction enrichment analysis was carried out through String (Damian et al., 2019). MCODE (Bader and Hogue, 2003) was also used to identify closely connected network components, and five functional components with the highest p-value were obtained. The PPI included 118 nodes and 144 interactions (Figure 6A). According to the annotation of MCODE, PPI enrichment analysis showed that five MCODE modules were enriched (Figure 6B) and the regulatory function of differentially expressed mRNAs were mainly concentrated in skeletal system development, tissue morphogenesis, MAPK signaling pathway regulation, etc (shown in Supplementary Table S1).

FIGURE 6

FIGURE 6. (A), all protein-protein interactions among differentially expressed mRNAs were extracted from PPI data source and formed a PPI network. (B), MCODE algorithm was applied to identify neighborhoods where proteins were densely connected. Five MCODE modules were enriched. PPI, protein-protein interaction; MCODE, Molecular Complex Detection.

3.4 Co-expression network and target prediction of lncRNA-mRNA

In order to explore the interaction between differentially expressed lncRNAs and mRNAs, we constructed a lncRNA-mRNA co-expression network with 226 nodes and 3,390 edges, according to |r| > 0.9 and p < 0.01 (Figure 7A). Moreover, we sorted the differentially expressed lncRNAs by connectivity degree, and the top 10 lncRNAs were listed in Table 2, such as FAM227B, POM121L9P, AF165147, AC103702, for further study.

FIGURE 7

FIGURE 7. (A), a lncRNA-mRNA co-expression network with 226 nodes and 3,390 edges was constructed, and the top 10 lncRNAs ranked by connectivity degree in the network were highlighted as squares. The larger the square, the higher connectivity degree. Blue squares represented lncRNAs and yellow circles represented mRNAs. (B), cis-regulation predictive analysis on the differentially expressed lncRNAs and mRNAs was performed. The top five lncRNAs with significant differences were listed separately, with green diamonds represented lncRNAs and yellow circles represented mRNAs.

TABLE 2

TABLE 2. The top 10 lncRNAs were sorted by connectivity degree.

3.5 Target gene prediction of lncRNA

Furthermore, cis-regulation predictive analysis on the differentially expressed lncRNAs and mRNAs was performed (Figure 7B). At the same time, according to the correlation coefficient, the ranking was carried out (Table 3). Comparing the results of prediction with the above 10 differentially expressed lncRNAs, it was found that FGF7 and HOXB4/5/6 might play an important role in the occurrence and development of CPT as targets of FAM227B and AC103702. These interactions between lncRNAs and mRNAs, as the possible mechanism of CPT, might become potential molecular targets for CPT therapy.

TABLE 3

TABLE 3. Cis-regulation predictive analysis on the differentially expressed lncRNAs and mRNAs was performed according to the correlation coefficient.

3.6 qRT-PCR

To verify the reliability of microarray, five differentially expressed lncRNAs were selected for qRT-PCR. The results showed that the expression of FAM227B, POM121L9P, TARID, SRP9, and AF165147 were downregulated (p < 0.05), which were consistent with microarray (Figure 8).

FIGURE 8

FIGURE 8. Five differentially expressed lncRNAs were selected for qRT-PCR. qRT-PCR, quantitative reverse transcription-polymerase chain reaction.

4 Discussion

In this study, high-throughput microarray sequencing technology was used to identify the genome of periosteal stem cells from normal and CPT patients, detect the differentially expressed lncRNAs and mRNAs, and statistically analyze their biological functions.

Since CPT has been reported, its pathogenesis remains to be unclear (BOYDHarold, 1982; Pannier, 2011). Through previous epidemiological and pathological studies, some scholars have found that 40%–80% of CPT patients are NF-1 gene mutation carriers, and there are neurofibroma-like tissues growing and invading in the diseased periosteum (Carlier et al., 2016; Van Royen et al., 2016). Histological studies have also found that the periosteum of the lesion site of CPT has abnormal thickening (Hermanns-Sachweh et al., 2005; Pannier, 2011), accompanied by fibrous tissue hamartoma-like hyperplasia, which was similar to the pathology of NF-1 patients (Van Royen et al., 2016). We also found in the previous study that the exosomes of CPT patients with NF-1 can significantly inhibit bone formation and stimulate bone absorption, thereby reducing the bone mass of trabecular and cortical bone in patients (Yang et al., 2021). However, some scholars suggested that whether CPT is combined with NF-1 has no significant difference in clinical characteristics, histopathology and postoperative tibial healing (Guanghui et al., 2019). At the same time, some scholars believed that the local bone vasculature disorder caused by the abnormal periosteal blood vessels in CPT is also one of the important factors for the formation of pseudarthrosis. The critical role of the vasculature in bone health has been repeatedly demonstrated in previous studies (Chen et al., 2021; Charlotte and Anjali, 2022). Benita et al. found that the periosteal arteries of the lesion were occluded, resulting in decreased osteogenic capacity in CPT patients (Hermanns-Sachweh et al., 2005). Tommi et al. confirmed that the vessel wall at the pseudarthrosis was thickened, the lumen was narrowed, and the local blood supply was blocked (Tommi et al., 2006).

With the continuous progress of molecular biology technology in recent years, our understanding of the pathogenesis of CPT is also deepening. Lei et al, (2001) found that there is excessive proliferation of fibroblasts and abnormal expression of related cytokines in periosteum of CPT lesions, which can result in abnormal thickening of periosteum, suggesting that CPT may be a disease originated from periosteum. Periosteum is composed of three layers: fibrous layer, undifferentiated layer and germinal layer. The fibrous layer is mainly composed of collagen fibers. The germinal layer is rich in PSCs, osteoblasts, mononuclear macrophages and other cells, with good osteogenic ability (Lin et al., 2014). The balance of bone metabolism is maintained by osteoblasts and osteoclasts, which are differentiated from PSCs under specific conditions. Studies have found that compared with normal children, diseased PSCs from CPT patients show lower osteogenic differentiation ability and faster aging speed (Cho et al., 2008; Bajada et al., 2009; Granchi et al., 2012). Some scholars have used stem cells in the treatment of local CPT and achieved good curative effect (Tikkanen et al., 2010). These results suggest that the periosteal lesions of CPT are caused by abnormal PSCs to some extent.

In recent years, high-throughput transcriptome sequencing technology has been increasingly recognized by scholars due to its fast speed, high throughput and good accuracy. By comparing and analyzing non-coding RNA such as lncRNA and miRNA between diseased and normal cells, key genes and regulatory pathways related to diseases can be screened out, which is one of the mainstream methods to explore the pathogenesis of diseases at present. Many studies have shown that non-coding RNA plays an important role in the differentiation and regulation of stem cells such as osteogenesis and adipogenesis. Zhuang et al. (Zhuang et al., 2015) sequenced the stem cells from ankylosing spondylitis and normal patients and found that Lnc ZNF354A-1, Lnc LIN54-1, Lnc FRG2C-3 were involved in abnormal osteogenic differentiation of MSCs in ankylosing spondylitis patients. Song et al. (Song et al., 2015) induced the osteogenic differentiation of BMSCs and found that there was a complex interaction between 217 lncRNAs and miRNAs, which affected the osteogenic differentiation of BMSC by regulating ECM receptor, adhesion spot and other pathways. Huang et al. (Huang et al., 2022) found that differentially expressed mRNA was mainly downregulated in CPT MSCs, which was related to the decline of osteogenic differentiation ability of MSCs. However, there are few studies on lncRNA in the pathogenesis of CPT. Therefore, by exploring the regulatory mechanism of lncRNA-mRNA in the occurrence and development of CPT, it may provide potential targets for the prevention and treatment of CPT. In this study, lncRNA-mRNA were identified in five pairs of normal and CPT PSCs by high-throughput microarray sequencing technology, and differentially expressed lncRNAs and mRNAs were screened out. The results showed that there were 194 differentially expressed lncRNAs (73 upregulated genes and 121 downregulated genes), and 822 differentially expressed mRNAs (311 upregulated genes and 511 downregulated genes). Both lncRNAs and mRNAs were mainly downregulated in CPT PSCs, which may be closely related to the decline of osteogenic differentiation ability and abnormal bone metabolism.

The enrichment analysis of differentially expressed mRNAs showed that the biological functions with the most significant differences were mainly concentrated in skeletal system development, tissue morphogenesis and NABA core matrisome. Signaling pathway analysis showed that the pathways with the most significant differences were mainly concentrated in pathways in cancer, PI3K-Akt signaling pathway and cytokine-cytokine receptor interaction, which was similar to the previous literature (Song et al., 2015). Moreover, through PPI enrichment analysis, we found that the function of differentially expressed mRNA regulatory proteins was related to skeletal system development, tissue morphogenesis, MAPK signaling pathway regulation, etc. The results further showed that the development and construction of skeletal system plays an important role in the pathogenesis of CPT.

LncRNA-mRNA co-expression network can directly show the interaction of the whole transcription level. According to the position of genes in the co-expression network, the key genes can be efficiently inferred and screened out, which may play an important role in the pathogenesis. In this study, a series of key genes, such as FAM227B, POM121L9P, AF165147, and AC103702, were listed by analyzing the differences of co-expression network and the changes of connectivity degree abundance. These 10 differentially expressed lncRNAs had the highest connectivity degree and were located in the center of the co-expression network, which means some of them might be closely related to the pathogenesis of CPT. Interestingly, these key lncRNAs were all downregulated in CPT PSCs. This was also consistent with the results of enrichment analysis mentioned above. At the same time, we predicted the targets of differentially expressed lncRNAs and ranked them according to the correlation coefficient. The results showed that FAM227B-FGF7 and AC103702-HOXB4/5/6 were potential target pathways with clinical significance. Homologous (HOX) gene is a relatively evolutionarily conserved transcription factor, which plays an important regulatory role in embryonic development (Deschamps and and van Nes, 2005). Studies have shown that HOX is a potential mechanism to maintain the tripotency of PSCs, and inhibition of HOX expression will lead to the loss of tripotency in stem cells (Jacopo et al., 2012). Bradaschia et al. found that HOX-positive periosteum contains more PSCs and shows better differentiation ability (Bradaschia-Correa et al., 2019). Liu et al. found that HOXB7, as the downstream target of miR-24, can block the inhibition of miR-24 on osteogenic differentiation of MSCs (Huang et al., 2022). FGF7 is a regulator of interstitial cell secretion, and its combination with FGFR plays an important role in regulating cell growth, migration and differentiation (Pastar et al., 2014; Liu et al., 2020). Previous studies have found that FGF7 can improve bone formation and promote osteogenic differentiation of stem cells (Poudel et al., 2017; Liu et al., 2018). Jeon et al. (2013) found that FGF7 can induce osteogenic differentiation of embryonic stem cells through ERK/Runx2 signaling pathway. Qiu et al. (2022) reported that over-expression of miR-381-3p inhibits osteogenic differentiation of BMSCs by targeting negative regulation of FGF7. All the above studies indicated that HOX and FGF7 genes play a positive role in the process of cell osteogenesis. In this study, it was found that both FAM227B and AC103702 showed a downregulation trend in CPT PSCs, while their downstream targets FGF7 and HOXB4/5/6 were also at low expression levels. It is reasonable to speculate that these two signaling pathways are closely related to the decline of osteogenic differentiation ability of CPT PSCs. Of course, the function of lncRNA/mRNA in CPT PSCs still remains largely unknown and needs further study.

There were still some limitations in this study. First of all, the sample size included in this study was small, so there might be some certain sample heterogeneity. The results could not fully represent the characteristics of all patients and need to be verified by larger samples in the future. Secondly, the basic data of patients could not be homogenized, so patients of different ages and types might have different gene expressions. Finally, on the basis of bio-informatics analysis, we need to further verify the function in animal experiments in the future, and further study the role of lncRNA-mRNA in the occurrence and development of CPT, so as to provide new ideas and scientific basis for the prevention and treatment of CPT.

In this study, the differential expression of lncRNAs and mRNAs in normal and CPT PSCs was detected by high-throughput sequencing analysis. The biological processes and signaling pathways of some differentially expressed lncRNAs and mRNAs were related to skeletal development and construction. In addition, by exploring lncRNA-mRNA co-expression and lncRNA target prediction, two significantly correlated lncRNA-mRNA signaling pathways were screened out. These results were helpful to deepen our understanding of the pathogenesis of CPT and provide potential therapeutic targets.

Data availability statement

The original contributions presented in the study are publicly available. This data can be found here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE200069

Ethics statemant

This study had been approved by the Ethics Committee of the First Affiliated Hospital of Zhejiang University, School of Medicine and Hunan Children’s Hospital.

Author contributions

ZL and GY designed and directed the study, ZL. and KL wrote the manuscript. ZL and HM were in charge of the experiment and collected results. HM and KL performed Statistical analysis. KL and GY critically reviewed the manuscript. All authors read and approved the final manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (82101818). Natural Science Foundation of Hunan Province (2022JJ70007). Zhejiang Medical and Health Science and Technology Project (2023RC157).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1094298/full#supplementary-material

Abbreviations

PSCs, periosteal stem cells; CPT, congenital pseudarthrosis of the tibia; lncRNA, long-chain non-coding RNA; mRNA, messenger RNA; GO, Gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; PPI, Protein -protein interaction; MCODE, Molecular Complex Detection; HOX, Homologous; NF-1, Type I neurofibromatosis.

References

Augustin, G., Antabak, A., and Davila, S. (2007). The periosteum. Part 1: Anatomy, histology and molecular biology. Injury-international J. Care Inj. 38 (10), 1115–1130. doi:10.1016/j.injury.2007.05.017

CrossRef Full Text | Google Scholar

Bader, G. D., and Hogue, C. W. V. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinforma. 4, 2. doi:10.1186/1471-2105-4-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Bajada, S., Marshall, M. J., Wright, K. T., Richardson, J. B., and Johnson, W. E. B. (2009). Decreased osteogenesis, increased cell senescence and elevated Dickkopf-1 secretion in human fracture non union stromal cells. Bone 45 (4), 726–735. doi:10.1016/j.bone.2009.06.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyd, , and Harold, B. (1982). Pathology and natural history of congenital pseudarthrosis of the tibia. Clin. Orthop. Relat. Res. 166 (166), 5–13. doi:10.1097/00003086-198206000-00003

PubMed Abstract | CrossRef Full Text | Google Scholar

Bradaschia-Correa, V., Kevin, L., Josephson Anne, M., Lee, S., Palma, L., Litwa, H. P., et al. (2019). Hox gene expression determines cell fate of adult periosteal stem/progenitor cells. Sci. Rep. 9, 5043. doi:10.1038/s41598-019-41639-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Carlier, A., Ashbourn, J. M. A., Nica, I., Legius, E., and Geris, L. (2016). Capturing the wide variety of impaired fracture healing phenotypes in neurofibromatosis type 1 with eight key factors: A computational study. Sci. Rep. 7, 20010. doi:10.1038/srep20010

PubMed Abstract | CrossRef Full Text | Google Scholar

Charlotte, O-W., and Anjali, K. (2022). Fundamentals of bone vasculature: Specialization, interactions and functions. Semin. Cell Dev. Biol. 123, 36–47. doi:10.1016/j.semcdb.2021.06.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, J., Unnikrishnan, S., Lih, T. S., Lippo, L., De Angelis, J., Labella, R., et al. (2021). High-resolution 3D imaging uncovers organ-specific vascular control of tissue aging. Sci. Adv. 7, eabd7819. doi:10.1126/sciadv.abd7819

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, T. J., Seo, J. B., Lee, H. R., Yoo, W. J., Chung, C. Y., and Choi, I. H. (2008). Biologic characteristics of fibrous hamartoma from congenital pseudarthrosis of the tibia associated with neurofibromatosis type 1. J. Bone & Jt. Surg. 90 (12), 2735–2744. doi:10.2106/JBJS.H.00014

PubMed Abstract | CrossRef Full Text | Google Scholar

Damian, S., Gable Annika, L., David, L., Junge, A., Wyder, S., Huerta-Cepas, J., et al. (2019). STRING v11: Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 47, D607–D613. doi:10.1093/nar/gky1131

PubMed Abstract | CrossRef Full Text | Google Scholar

Deschamps and, J., and van Nes, J. (2005). Developmental regulation of the Hox genes during axial morphogenesis in the mouse. Development 132 (13), 2931–2942. doi:10.1242/dev.01897

PubMed Abstract | CrossRef Full Text | Google Scholar

Doherty, L., Wan, M., Kalajzic, I., and Sanjay, A. (2020). Diabetes impairs periosteal progenitor regenerative potential. Bone 143 (1), 115764. doi:10.1016/j.bone.2020.115764

PubMed Abstract | CrossRef Full Text | Google Scholar

Filion, T., Skelly, J. D., Huang, H., Greiner, D. L., Ayers, D. C., and Song, J. (2017). Impaired osteogenesis of T1DM bone marrow-derived stromal cells and periosteum-derived cells and their differential in-vitro responses to growth factor rescue. Stem Cell Res. Ther. 8 (1), 65. doi:10.1186/s13287-017-0521-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Granchi, D., Devescovi, V., Baglìo, S., Leonardi, E., Donzelli, O., Magnani, M., et al. (2010). Biological basis for the use of autologous bone marrow stromal cells in the treatment of congenital pseudarthrosis of the tibia. Bone 46 (3), 780–788. doi:10.1016/j.bone.2009.10.044

PubMed Abstract | CrossRef Full Text | Google Scholar

Granchi, D., Devescovi, V., Baglio, S. R., Magnani, M., Donzelli, O., and Baldini, N. (2012). A regenerative approach for bone repair in congenital pseudarthrosis of the tibia associated or not associated with type 1 neurofibromatosis: Correlation between laboratory findings and clinical outcome. Cytotherapy 14 (3), 306–314. doi:10.3109/14653249.2011.627916

PubMed Abstract | CrossRef Full Text | Google Scholar

Guanghui, Z., Zheng, Y., Yaoxi, L., Yan, A., Zhengmao, H., Yongjia, Y., et al. (2019). Identification and characterization of NF1 and non-NF1 congenital pseudarthrosis of the tibia based on germline NF1 variants: Genetic and clinical analysis of 75 patients. Orphanet J. Rare Dis. 14 (1), 221. doi:10.1186/s13023-019-1196-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Hefti, F., Bollini, G., Dungl, P., Fixsen, J., Grill, F., Ippolito, E., et al. (2000). Congenital pseudarthrosis of the tibia: History, etiology, classification, and epidemiologic data. J. Pediatr. Orthopaedics-part B 9 (1), 11–15. doi:10.1097/01202412-200001000-00003

PubMed Abstract | CrossRef Full Text | Google Scholar

Hermanns-Sachweh, B., Senderek, J., Alfer, J., Klosterhalfen, B., Buttner, R., Fuzesi, L., et al. (2005). Vascular changes in the periosteum of congenital pseudarthrosis of the tibia. Pathology - Res. Pract. 201 (4), 305–312. doi:10.1016/j.prp.2004.09.013

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Y., Yang, G., and Liu, Y. (2022). Proteomic analysis of serum-derived exosomes from congenital pseudarthrosis of tibia. J. Clin. Pediatr. Surg. 21 (06), 551–557. doi:10.3760/cma.j.cn101785-202102043-010

CrossRef Full Text | Google Scholar

Iyer Matthew, K., Niknafs Yashar, S., Rohit, M., Singhal, U., Sahu, A., Hosono, Y., et al. (2015). The landscape of long noncoding RNAs in the human transcriptome. Nat. Genet. 47, 199–208. doi:10.1038/ng.3192

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacopo, P., Trombi, L., Spugnesi, L., Barachini, S., Maroni, G., Brodano, G. B., et al. (2012). HOXandTALEsignatures specify human stromal stem cell populations from different sources. J. Cell. Physiology 228 (4), 879–889. doi:10.1002/jcp.24239

CrossRef Full Text | Google Scholar

Jeon, Y. M., Kook, S. H., Rho, S.-J., Lim, S. S., Choi, K. C., Kim, H. S., et al. (2013). Fibroblast growth factor-7 facilitates osteogenic differentiation of embryonic stem cells through the activation of ERK/Runx2 signaling. Mol. Cell. Biochem. 382 (1-2), 37–45. doi:10.1007/s11010-013-1716-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Lei, W. (2001). Histopathological study of congenital pseudarthrosis of tibia. J. Fourth Mil. Med. Univ. 22 (20), 3.

Google Scholar

Li, Z., Yu, H., Huang, Y., Liu, Y., Zhu, G., Tan, Q., et al. (2022). Analysis of risk factors affecting union and refracture after combined surgery for congenital pseudarthrosis of the tibia: A retrospective study of 255 cases. Orphanet J. rare Dis. 17 (1), 245. doi:10.1186/s13023-022-02375-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, Z., Fateh, A., Salem, D. M., and Intini, G. (2014). Periosteum: Biology and applications in craniofacial bone regeneration. J. Dent. Res. 93 (2), 109–116. doi:10.1177/0022034513506445

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X-Y., Xin, L., Ming-Ru, B., Chen, X., Wang, C. L., Xie, J., et al. (2018). FGF-7 dictates osteocyte cell processes through beta-catenin transduction. Sci. Rep. 8, 14792. doi:10.1038/s41598-018-33247-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Z., Deng, Y., Li, T., Zhu, F., and Zhou, X. (2020). The opposite functions of miR-24 in the osteogenesis and adipogenesis of adipose-derived mesenchymal stem cells are mediated by the HOXB7/β-catenin complex. FASEB J. 34 (7), 9034–9050. doi:10.1096/fj.202000006RR

PubMed Abstract | CrossRef Full Text | Google Scholar

Macdonald, W. A., and Mann, M. (2020). Long noncoding RNA functionality in imprinted domain regulation. PLoS Genet. 16 (8), e1008930. doi:10.1371/journal.pgen.1008930

PubMed Abstract | CrossRef Full Text | Google Scholar

Matthews Brya, G., Novak, S., Sbrana Francesca, V., Funnell, J. L., Cao, Y., Buckels, E. J., et al. (2021). Heterogeneity of murine periosteum progenitors involved in fracture healing. Elife 10, e58534. doi:10.7554/eLife.58534

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, S., and Lee, K. W. (2015). From discovery to function: The expanding roles of long non-coding RNAs in physiology and disease. Endocr. Rev. 36 (1), 25–64. doi:10.1210/er.2014-1034

PubMed Abstract | CrossRef Full Text | Google Scholar

Mila, G., and Alessia, C. (2019). Long noncoding RNA and epithelial mesenchymal transition in cancer. Int. J. Mol. Sci. 20, 1924. doi:10.3390/ijms20081924

PubMed Abstract | CrossRef Full Text | Google Scholar

Pannier, S. (2011). Congenital pseudarthrosis of the tibia. Orthop. Traumatology Surg. Res. 97 (7), 750–761. doi:10.1016/j.otsr.2011.09.001

CrossRef Full Text | Google Scholar

Pastar, I., Stojadinovic, O., Yin, N. C., Ramirez, H., Nusbaum, A. G., Sawaya, A., et al. (2014). Epithelialization in wound healing: A comprehensive review. Adv. wound care 3 (7), 445–464. doi:10.1089/wound.2013.0473

PubMed Abstract | CrossRef Full Text | Google Scholar

Pearson, M. J., Philp, A. M., Heward, J. A., Roux, B. T., Walsh, D. A., Davis, E. T., et al. (2016). Long intergenic noncoding RNAs mediate the human chondrocyte inflammatory response and are differentially expressed in osteoarthritis cartilage. Arthritis & Rheumatology 68 (4), 845–856. doi:10.1002/art.39520

CrossRef Full Text | Google Scholar

Pertea, M., Pertea, G. M., Antonescu, C. M., Chang, T. C., Mendell, J. T., and Salzberg, S. L. (2015). StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 33 (3), 290–295. doi:10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

Poudel, S., Bhattarai, G., Kim, J. H., Kook, S. H., Seo, Y. K., Jeon, Y. M., et al. (2017). Local delivery of recombinant human FGF7 enhances bone formation in rat mandible defects. J. bone mineral metabolism 35 (5), 485–496. doi:10.1007/s00774-016-0784-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, L., Cai, J., Zhang, N., Fan, F. Y., and Li, X. M. (2022). Effect of miR-381-3p/FGF7 axis on the osteogenic differentiation of bone marrow mesenchymal stem cells through MEK/ERK signaling pathway. Tissue & Cell 76, 101791. doi:10.1016/j.tice.2022.101791

CrossRef Full Text | Google Scholar

Smyth, G. K., and McCarthy, D. J. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26 (1), 139–140. doi:10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, W. Q., Gu, W. Q., Qian, Y. B., Ma, X., Mao, Y. J., and Liu, W. J. (2015). Identification of long non-coding RNA involved in osteogenic differentiation from mesenchymal stem cells using RNA-Seq data. Genet. Mol. Res. Gmr 14 (4), 18268–18279. doi:10.4238/2015.December.23.14

PubMed Abstract | CrossRef Full Text | Google Scholar

Tachiwana, H., Yamamoto, T., and Saitoh, N. (2020). Gene regulation by non-coding RNAs in the 3D genome architecture. Curr. Opin. Genet. Dev. 61, 69–74. doi:10.1016/j.gde.2020.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Tikkanen, J., Leskela, H. V., Lehtonen, S. T., Vahasarja, V., Melkko, J., Ahvenjarvi, L., et al. (2010). Attempt to treat congenital pseudarthrosis of the tibia with mesenchymal stromal cell transplantation. Cytotherapy 12 (5), 593–604. doi:10.3109/14653249.2010.487898

PubMed Abstract | CrossRef Full Text | Google Scholar

Tommi, K., Pentti, K., Marja, N., Alanne, M., Leskela, H. V., Lehenkari, P., et al. (2006). Vasculopathy in two cases of NF1-related congenital pseudarthrosis. Pathol. Res. Pract. 202, 687–690. doi:10.1016/j.prp.2006.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Van Royen, K., Brems, H., Legius, E., Lammens, J., and Laumen, A. (2016). Prevalence of neurofibromatosis type 1 in congenital pseudarthrosis of the tibia. Eur. J. Pediatr. 175 (9), 1193–1198. doi:10.1007/s00431-016-2757-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H., Yin, Q. F., Luo, Z., Yao, R. W., Zheng, C. C., Zhang, J., et al. (2016). Unusual processing generates SPA LncRNAs that sequester multiple RNA binding proteins. Mol. Cell 64, 534–548. doi:10.1016/j.molcel.2016.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, G., Yu, H., Liu, Y., Ye, W., Zhu, G., Yan, A., et al. (2021). Serum-derived exosomes from neurofibromatosis type 1 congenital tibial pseudarthrosis impaired bone by promoting osteoclastogenesis and inhibiting osteogenesis. Exp. Biol. Med. 246 (2), 130–141. doi:10.1177/1535370220962737

CrossRef Full Text | Google Scholar

Zhou, Y., Tan, Q., Liu, K., Liu, Y., Zhu, G., Mei, H., et al. (2022). Epidemiological and clinical characteristics of congenital pseudarthrosis of the tibia in China. Front. Pediatr. 10, 943917. doi:10.3389/fped.2022.943917

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhuang, W., Ge, X., Yang, S., Huang, M., Chen, P., et al. (2015). Upregulation of lncRNA MEG3 promotes osteogenic differentiation of mesenchymal stem cells from multiple myeloma patients by targeting BMP4 transcription. Stem Cells 33 (6), 1985–1997. doi:10.1002/stem.1989

PubMed Abstract | CrossRef Full Text | Google Scholar