Transcriptome Sequencing Analysis of lncRNA and mRNA Expression Profiles in Bone Nonunion

Background Bone nonunion is a serious complication of fracture. This study explored the differentially expressed lncRNAs (DELs) and mRNAs (DEGs) and identified potential lncRNA-mRNA interactions in bone nonunion. Methods We extracted total RNA from three bone nonunion and three bone union patient tissue samples. RNA sequencing was performed to detect DELs and DEGs between bone nonunion and union tissue samples. The lncRNAs and genes with absolute log2-fold change (log2FC) > 1 and adjusted p value < 0.05 were further chosen for gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. lncRNA and targeted mRNA interaction networks were constructed. Results We observed 179 DELs and 415 DEGs between the bone nonunion and union tissue samples. GO analysis indicated that DELs and DEGs were mainly enriched in the chondroitin sulfate proteoglycan biosynthetic process. DELs and DEGs were enriched in “ECM-receptor interaction” and “Staphylococcus aureus infection” KEGG pathways. Several potential lncRNA-mRNA interactions were also predicted. Conclusions This study identified bone nonunion-associated lncRNAs and mRNAs using deep sequencing that may be useful as potential biomarkers for bone nonunion.


Introduction
Bone nonunion, a serious complication of fracture, occurs in approximately 5-10% of patients with bone fractures [1][2][3][4][5][6][7]. Infected bone nonunion is caused by many factors, including fractures and accidents. Bone nonunion may lead to delayed union and amputation, which further contributes to functional limitation, disability, and poor quality of life [1,2]. The most common causes of bone nonunion include infection, insufficient local blood supply, separation of fracture ends, and insufficient fracture stabilization [8,9]. The use of antibiotics has improved the treatment of bone infections, but bone nonunion remains an obstacle in the repair of damaged bone [3,4].
Currently, nonunion is a serious challenge in the treatment of bone loss associated with bone infections. The process of bone remodeling includes the breakdown and resorption of bone and the formation of new bone. Osteoclasts are responsible for bone breakdown and resorption, whereas osteoblasts are responsible for new bone formation. Osteoclasts attach to the older bone area, secrete acidic substances to dissolve minerals, secrete protease to digest the bone matrix, and form bone resorption lacuna. Subsequently, osteoblasts migrate to the resorbed site and secrete the bone matrix that is then mineralized to form new bone. The balance between osteoclastic and osteogenic processes is substantial in maintaining the normal bone mass. However, this balance is compromised because resorption replaces formation in case of bone infection, resulting in bone loss or bone nonunion [10]. A previous study indicated that differentially expressed miRNAs might be a potential diagnostic and therapeutic biomarker for infected tibial nonunion [11]. Additionally, the data from the GEO dataset indicated that ADAMTS18 and TGFBR3 genes were differentially expressed in nonunion skeletal fracture [12]. Moreover, the coinjection of BMP and DCN into the bone nonunion area could improve the induction of bone formation [13,14]. However, the exact molecular mechanisms underlying bone nonunion remain unclear at present. Therefore, it is critical to explore the etiological mechanism of bone nonunion and to develop new targets for the diagnosis and treatment of infected bone nonunion.
Long noncoding RNAs (lncRNAs) are a class of RNAs longer than 200 nucleotides. Previous studies have indicated that lncRNAs can act as master regulators, affecting target gene expression levels [15]. Growing evidence suggests that lncRNAs are involved in the epigenetic regulation of gene expression, transcription, cell death, and other important biological processes [16,17]. Previous reports showed that many lncRNAs are related to osteoclast and osteoblast cell functions. For instance, lncGHET1, lncRhno1, lncTUG1, and lncUCA1 are identified to be involved with osteoblast proliferation and differentiation [18][19][20][21], whereas lncXIST [22], lncNeat1 [23], and lncCRNDE [24] are reported to be associated with osteoclast differentiation. However, at present, there is insufficient information on specific lncRNAs involved in bone nonunion.
In this study, we obtained bone nonunion and union tissue samples from patient fracture sites. We performed transcriptome sequencing of these tissues to determine the differentially expressed lncRNAs (DELs) and differentially expressed mRNAs (DEGs) and identify potential lncRNA-mRNA interactions. Our findings may provide new insights to further elucidate the pathogenesis of, and develop biomarkers for, bone nonunion.

Sample Collection.
The samples used in this study were obtained from three patients with normal fracture healing and three patients with bone nonunion ( Table 1). The specimens were collected from the normal healing fracture site and scar tissue approximately 3 mm in size at the bone nonunion site. The diagnosis of bone nonunion was based on the definition given by the Food and Drug Administration. First, the fractures went unhealed for six months and there was no further healing trend within three months. Clinical X-ray examination was performed to confirm bone nonunion, and surgery further confirmed the formation of a small amount of scar tissue and callus at the fracture end or with only a small amount of scar tissue. None of the patients included in the study had infections, tumors, autoimmune diseases, bone nonunion caused by pathological fractures, history of hormone use, and history of smoking.
2.2. Total RNA Extraction. We extracted total RNA using TRIzol reagent following the manufacturer's protocol. The absorbance ratio at 260/280 nm (A260/A280) was measured using SmartSpec Plus to determine the concentration and purity of the isolated RNA. The integrity of the extracted RNA was confirmed using electrophoresis (1.5% agarose gel). The RNA was then transcribed into first-strand cDNA using the First-Strand cDNA Synthesis Kit (TaKaRa) for performing gene expression analysis.
2.3. lncRNA and mRNA Sequencing. We used 3 μg RNA per sample for sample preparation. Following ribosomal RNA (rRNA) depletion, the RNA was fragmented and a cDNA library was constructed using the VAHTS Total RNA-seq (HMR) Library Prep Kit. Libraries were sequenced on an Illumina HiSeq 2500 platform according to the manufacturer's instructions, and 125 bp paired-end reads were produced (Table S1 & S2). DELs and DEGs between samples were identified using the Cuffdiff program in the Cufflinks package. As cutoff criteria, p values < 0.05 and jlog 2FCj > 1 were used.

Analysis of DELs and DEGs. Gene Ontology (GO)
Enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for these DEGs and predicted target genes for DELs were conducted using the clusterProfiler R package. We obtained all the gene sets used in these functional annotations from the DAVID database. p values were adjusted by Benjamini & Hochberg methods, and FDR < 0:05 was defined as significantly enriched.

Prediction of Cis-and Trans-Regulated Target Genes of
DELs. lncRNAs directly regulate adjacent target genes in the genome and this is termed cis-acting regulation. According to the taxonomic annotation information of lncRNA, neighboring known genes are predicted to be potential cisregulated target genes. lncRNAs that are located far from their target genes play an indirect regulatory role through sequence complementarity, which is referred to as transacting regulation. We used RepeatMasker to search the Alu repeat structure of lncRNAs and 3 ′ -UTRs. We used BLASTN sequence alignment to search for complementary sequence regions of lncRNAs and 3 ′ -UTRs. The thermodynamic stability and binding ability of complexes formed by lncRNAs and 3 ′ -UTRs were predicted by RNAplex and RIsearch, with an aim to predict trans-regulated target genes of lncRNAs.

Boxplot and Principal Component Analysis (PCA)
Diagram. Boxplots describe data using five statistics, Oxidative Medicine and Cellular Longevity including the minimum, first quartile (25%), median (50%), third quartile (75%), and maximum. Through boxplots, we can gauge the symmetry of the data and the degree of dispersion of the distribution. As shown in Figure 1(a), we observed that the gene expression level of samples 1-5 was stable, while sample N6 was different, which may have been caused by a more serious fracture in patient 6. It is possible to observe the similarity between samples through PCA plots. The closer the distance between samples on the PCA diagram, the closer the expression trend of sample genes is. As shown in Figure 1(b), the PCA diagram revealed that the expression features of samples 1-5 were similar, while sample N6 was different, which may be caused by a more serious fracture in patient 6.  Tables 2 and 3, respectively. Additionally, we annotated and classified the studied lncRNAs and they were mainly divided into intergenic lncRNAs, sense lncRNAs, intronic lncRNAs, antisense lncRNAs, sRNA host lncRNAs, enhancer lncRNAs, and bidirectional lncRNAs, accounting for 75.1%, 15.3%, 3.4%, 2.8%, 1.7%, and 0.6%, respectively ( Figure S1).

Function and Pathway Predictive Analysis of DELs and
DEGs. The biological processes (BP), cellular components (CC), and molecular functions (MF) of DEGs and DELs were analyzed using the DAVID database. GO analysis of DELs in terms of MF showed carbohydrate derivative transporter activity and Se-containing compound transmembrane transporter activity as most enriched. GO analysis of DELs in terms of CC was enriched in the CD95 deathinducing signaling complex in cellular components and smooth muscle contractile fiber. GO analysis of DELs in terms of BP was predominantly enriched in the regulation of peptidase activity and dermatan sulfate proteoglycan metabolic process ( Figure 4).
GO analysis of DEGs for MF was primarily enriched in MHC class II receptor activity and chemokine activity, CC showed enrichment in cell surface and plasma membrane components, and BPs showed developmental processes and cell migration ( Figure 5). 3.4. lncRNA Cis-and Trans-Regulated Genes. Differential expression of lncRNA cis-and trans-regulated genes was predicted. As shown in Figure 6, lncRNA ENST0000

Discussion
With the development of sequencing technology, transcriptome sequencing has been used to understand a variety of diseases [25,26]. Wei et al. used an miRNA expression profile of bone nonunion and union tissues to find nine upregulated and nine downregulated miRNAs [27]. Long et al. reported 557 differentially expressed miRNAs in bone nonunion tissues and further explored that miR-381 can modulate human bone mesenchymal stromal cell osteogenesis [28]. The results obtained in different transcriptome sequencing studies may vary greatly, which may be related to different sequencing technologies, samples, and sequencing methods. Previous studies indicated that lncRNAs achieve their functions in tumors through a wide range of mechanisms [29][30][31]. However, lncRNAs have been rarely studied in orthopedics, especially with respect to bone nonunion [32], thus limiting the detection and treatment of bone nonunion to a certain extent.
In this study, transcriptome sequencing was performed on bone tissue samples collected from long bones (tibia, femur, and humerus) of patients with bone nonunion and normal bone union. We detected and analyzed 179 DELs and 415 DEGs. GO analysis showed that DELs were primarily enriched in carbohydrate derivative transporter activity in MF, CD95 death-inducing signaling complex in CC, and regulation of peptidase activity in BP. The DEGs were mainly involved in MHC class II receptor activity for MF, cell surface, and developmental processes for CC and BP. The KEGG pathway enrichment of the DELs showed the ECM-receptor interaction pathway and viral carcinogenesis pathway. The KEGG pathway enrichment in DEGs showed the S. aureus infection pathway and rheumatoid arthritis pathway. The ECM-receptor interaction pathway primarily functions through three ECM proteins, including collagen, fibronectin, and laminin. Laminin is involved in osteogenesis and promotion of bone defect repair [33,34]. Studies have suggested that collagen type XV may be involved in ECM organization early in the osteogenesis process, a prerequisite for promoting subsequent mineral matrix deposition [35]. Immunohistochemical and transcriptomic studies have shown the expression and dynamic regulation of fibronectin in several stages of      Oxidative Medicine and Cellular Longevity fracture healing [35,36]. Single-cell RNA sequencing of the injury site revealed an early increase in mesenchymal progenitor cell (MPC) genes associated with cell adhesion pathways and ECM receptor interactions. The ECM creates a microenvironment with a MPC differentiation bias closer to a specific stiffness role in tissues [37][38][39][40]. For example, a rigid environment that mimics the natural bone favors differentiation into osteoblasts [41], whereas a softer matrix promotes the development of adipocyte fate [42,43].
lncRNAs can cis-regulate the transcription of adjacent protein-coding genes, thereby regulating the expression of such genes and participating in developmental and other biological processes associated with them. Cis-regulation refers to the transcriptional activation and expression regulation of noncoding RNAs to adjacent mRNAs. In this study, we found that lncRNA ENST00000453060 may cisregulate COL6A1. Previous studies have indicated that genetic deletion of COL6A1 impairs osteoblast connections and communication [44]. COL6A1 plays a substantial role     Oxidative Medicine and Cellular Longevity in osteoblasts, and lncRNA ENST00000453060 may regulate osteoblasts via cis-regulation of COL6A1. Additionally, our previous study confirmed that lncRNA ENST00000563492 could promote the osteogenesisangiogenesis coupling process in bone marrow stromal cells [45]. lncRNA trans-regulation is the regulation of distal mRNA transcription. lncRNAs can regulate the expression of distant genes by binding to enhancers and promoters. lncRNAs regulate the activity of bound proteins or RNAs in the cytoplasm or nucleus in a dose-dependent manner. The lncRNA UCSC_TCONS_00017121 trans-regulates RBP1. Previous studies have shown that RBP1 promotes differentiation of osteoblasts [46]. The lncRNA UCSC_ TCONS_00017121 may also regulate osteoblasts via cisregulating COL6A1.
There were certain limitations in the present study. First, no validation assays, including qPCR and histological analysis, were performed to confirm the differential expression, thereby demanding the need for further experimental studies. Second, patient matching, including differences in ages of enrolled patients and differences involving sites of sample collection from bone nonunion tissues, was not well handled. Third, only 6 patients were enrolled in this study, the results from which need to be confirmed in a further study with greater numbers of samples. Fourth, the sample size of this study was relatively small and the results need to be interpreted with caution.

Conclusions
A total of 179 DELs and 415 DEGs were identified between bone nonunion and bone union tissue samples. All of these lncRNAs and mRNAs may be related to the occurrence and development of bone nonunion. GO, KEGG, and regulatory analysis for these lncRNAs and mRNAs were performed to detect their potential functions. This study identified potential biomarkers for bone nonunion, but a validation cohort is still essential to confirm the applicability of these biomarkers.

Data Availability
The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

Ethical Approval
This study was approved by the Second Xiangya Hospital of Central South University Committee for Clinical Research, and all of the methods were in accordance with the

Consent
All studies included in this study got informed consent from each study participant, and each study was approved by the ethics committee or institutional review board.

Disclosure
The study funders/sponsors had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication.

Conflicts of Interest
The authors declare that they do not have any competing interests.