Differential Expression Pro les and Function Prediction of tRNA-Derived Fragments in Fibrous Dysplasia


 Background: Fibrous dysplasia (FD) is a benign bone disease which normal bone matrix is replaced by fibrous tissue and immature bone tissue. Transfer RNA (tRNA)-derived RNA fragments (tRFs) and tRNA halves (tiRNAs) are a type of small non-coding RNA in the transcriptome of eukaryotes that produced by specific shearing of mature tRNA. Here, we conducted a comparative analysis of the expression of tRFs/tiRNAs in BMSCs and FD BMSCs using a high‐throughput sequencing technique. Quantitative real‐time polymerase chain reaction (qRT-PCR) was used to validate the differential expressed tRFs/tiRNAs between two samples. Results: The results showed that tRF-34-87R8WP9I1EWJIQ [tDR‐001276] was significantly upregulated in FD BMSCs, and 3 downregulated tDRs (tRNA-derived small RNAs) (tRF-22-8EKSP1852 [tDR‐006826], tRF-18-H9R8B7D2 [tDR‐006049] and tRF-33-86V8WPMN1E8Y0E [tDR‐001271]) were also detected. Prediction of target genes and gene ontology (GO) and KEGG pathway analysis indicated that the upregulated tRF was mainly involved in regulation of immune response and osteoclast differentiation, which may be the underlying mechanism of FD pathological features. The downregulated tRFs/tiRNAs were related to calcium ion transport (tDR-006826), apoptotic signaling pathway and cell proliferation (tDR-006049) and endocrine system development (tDR-001271). The upregulated tRFs/tiRNAs was related to immune response (tDR-001276). The protein-protein interaction network analysis for predicted target genes established by the STRING database showed that PPP2R5A, ADAMTS1, PPARA, and POLR2C were the most frequently interacted proteins in target genes of tDR-006826, tDR-006049, tDR-001271 and tDR-001276, respectively. Conclusions: Our findings provided a comprehensive analysis of the expression of tRFs/tiRNAs in FD BMSCs and BMSCs. These differential expression tRFs/tiRNAs may be novel regulatory factors involved in FD BMSCs, and they could serve as potential therapeutic targets.


Background
Fibrous Dysplasia (FD) is a skeletal disorder resulting in local bulging deformities, pain, and pathological fractures, which seriously affect the health and life quality of patients [1,2]. FD is recognized as a stem cell disease which normal bone tissue is replaced by excessively proliferated brous tissue caused by mosaic mutations in GNAS gene [3,4]. Mutations in the α subunit of activated G protein leading to excessive accumulation of cyclic adenosine monophosphate (cAMP) in FD bone marrow stromal cells (BMSCs), which is the key cause of pathogenesis [5]. The application of osteoclast inhibitory drugs is the main way to relieve symptoms, but there is still a lack of effective and safe treatment methods [6].
Non-coding RNAs (ncRNAs) play signi cant roles in different disease pathogenesis. Recent advances in high-throughput RNA sequencing technology have revealed a class of ncRNAs which are derived from various RNA species such as ribosomal RNA, messenger RNA (mRNA), and transfer RNA (tRNA) [7][8][9].
It has been found that tRFs/tiRNAs are closely related to many human diseases such as tumors, metabolic diseases, neuropsychiatric diseases, and infectious diseases [16]. There have been accumulating studies in tRFs/tiRNAs which are involved in the regulation of biogenesis processes in different types of mammalian stem cells [8]. Recent studies have identi ed that tRFs play important roles in regulating translation, as well as microRNAs [15,17,18]. Moreover, some tRFs are also reported mediating gene silencing by binding to proteins [12,19]. Study also uncovers a critical function of epigenetic modi cation of tRFs in directing regulation of translation in stem cells. An additional posttranscriptional layer of regulation has been discovered which may directly affect the spatiotemporal regulation of gene expression during development and diseases [20].
The functions of messenger RNAs in FD have been well studied, however, there is no research on the role of tRFs in FD. In the present study, small RNA (miRNA) sequencing analysis was performed to detect tRFs/tiRNAs in BMSCs and FD BMSCs, and the differential expression tRFs/tiRNAs were further validated by real-time quantitative polymerase chain reaction (qRT-PCR). Our ndings indicate that tDR-006826, tDR-006049 and tDR-001271 are downregulated in FD BMSCs, whereas tDR-001276 is upregulated.
Moreover, the biological functions of four tiRNAs were evaluated to reveal potential pathological mechanisms of FD. Taken together, our ndings might provide a theoretical basis for exploring novel therapeutic targets of FD.

Expression pro les of tRFs/tiRNAs in FD BMSCs and normal BMSCs
We used a high-throughput sequencing technique to detect expression pro les of tRFs/tiRNAs in FD BMSCs and normal BMSCs. For the six libraries, an average of 10 million raw reads per library was obtained, whereas adapter-trimmed reads (length >= 16nt) were about 8.8 million (88.79%) (Additional le 1). Among the mapped reads, miRNAs were the most common with an average of 3154113 reads (to adapter-trimmed reads, 35.54%), and tRFs had an average of 1117973 reads (to adapter-trimmed reads, 12.60%) (Additional le 1). The copy number and sequence were recorded for each unique read and showed in the sequence read length distribution (Fig. 1a, b). tRNAs with different sequences may have the same anticodon and transfer the same amino acid. Stacked plot showed the number of tRFs/tiRNAs derived from the same anticodon tRNA in FD BMSCs and normal BMSCs (Fig. 1c, d).

Differential expression of tRFs/tiRNAs in FD BMSCs and BMSCs
The differential expression tRFs/tiRNAs analysis was performed between FD BMSCs and BMSCs, and the results were shown in the heat map (Fig. 2a). The principal component analysis (PCA), which showed distinguishable tRFs/tiRNAs expression pro les among 6 samples, was used to reduce the dimensionality of the data sets (Fig. 2b). The scatter plot assessed tRFs/tiRNAs expression variation between FD BMSCs and BMSCs (Fig. 2c), and the upregulated and downregulated tRFs/tiRNAs are showed in Fig. 2d. The Venn diagram exhibited that there were 797 types of tRFs/tiRNAs speci cally in BMSCs, and 17 types of tRFs/tiRNAs speci cally in FD BMSCs (Fig. 2d). We analyzed the percentage of each subtype of the tRFs/tiRNAs expressed in FD BMSCs and BMSCs, and the pie chart revealed that FD BMSCs mainly increased the expression of tiRNA-5, whereas decreased the expression of tRF-1 and tRF-3 (Fig. 2e).

Validation of differential expression tRFs/tiRNAs in BMSCs and FD BMSCs
In order to validate the differential expression, 3 signi cantly downregulated tRFs

PPI network analysis
Protein-protein interaction (PPI) network for the target genes was established using the STRING database (https://www.string-db.org) (Fig. 6). We analyzed the PPI results and found that PPP2R5A, ADAMTS1, PPARA and POLR2C were the most frequently interacted proteins in target genes of tDR-006826, tDR-006049, tDR-001271 and tDR-001276, respectively.

Discussion
FD is a slowly progressing bone tumor-like lesion which main pathological manifestation is that normal bone tissue is replaced by hyperproliferative brous tissue and dysplastic immature trabecular bone [21]. Existing studies suggest that FD is a stem cell disease [22]. Mutation of activated guanine nucleotide binding protein G protein α subunit (Gsα) in local FD leads to over-accumulation of the second messenger cAMP and activates downstream signaling pathways, which further causes the formation of amounts of brous tissue and immature bone tissue [23]. The clinical symptoms of FD which occur in craniomaxillofacial bones often show local bulging deformities, occlusal disorders, and pathological fractures [24]. Although studies are gradually deepening, there is still a lack of effective and safe treatment methods for FD.
tRNA is considered to be an adaptor molecule that helps ribosomes to decode mRNA and synthesize protein [25]. Recent studies have revealed that tRNAs are a major source of small non-coding RNAs that possess varied functions [26]. Studies have shown that in many species of cells, mature tRNA or precursor tRNA is speci cally sheared to produce tRNA-derived fragments (tRFs) or tRNA half-molecules (tiRNAs) [27,28]. tRFs/tiRNAs are widely present in tissue cells of various organisms and have tissue speci city and disease relevance [29]. Although the speci c biological functions of tRFs/tiRNAs have not been fully elucidated, more and more evidence suggests that they have a variety of biological functions, such as binding protein to affect the stability of mRNA, interaction with cytochrome C regulating apoptosis, changing the gene transcription cascade process of offspring, exerting negative gene regulation in the form of miRNA [13]. It has been found that tRFs/tiRNAs are closely related to many human diseases [30]. In our present study, high-throughput RNA sequencing detected four differential expression tRFs/tiRNAs in FD BMSCs and BMSCs. Meanwhile, one high-expression tiRNAs and three lowexpression tRFs/tiRNAs were veri ed via qRT-PCR, which were consistent with the high-throughput RNA sequencing data. Furthermore, bioinformatics analysis was performed to demonstrate the biological functions of the differential expression tRFs/tiRNAs.
In order to investigate the functions of tRFs/tiRNAs targets, we employed GO and KEGG analysis. Most of the results were enriched in processes related to calcium ion transport (tDR-006826), apoptotic signaling pathway and cell proliferation (tDR-006049), endocrine system development (tDR-001271), and immune response (tDR-001276). In our previous study, we have found that FD is relating to the proliferation of BMSCs which contributes to abnormal brous tissue proliferation [31]. Meanwhile, FD BMSCs was also exhibited weak apoptotic and strong proliferation ability [31]. When FD presents with skin pigmentation and precocious puberty, it is known as McCune-Albright, which inferring FD is an endocrine relating disease [24]. Pain is a common clinical manifestation in FD and in uences the health-related quality of life of FD patients [32]. Local in ammation may be the main cause of pain. Therefore, the four tRFs/tiRNAs may play important roles in FD.
Although the speci c biological functions of tRFs/tiRNAs have not been fully elucidated, more and more evidences show that they have a variety of biological functions. Some studies have shown that tRFs/tiRNAs affect the stability of mRNA via binding protein [13,33]. Other studies have observed that tRFs/tiRNAs perform a miRNA-like negative gene-regulation mechanism [34,35]. It is also report that tRFs could form a double-stranded molecule with complementary RNA, acting as a primer binding site or base pairing with the target RNA in Argonaute [12,35]. Because of the incomplete of tRNA database, the target prediction of tRFs/tiRNAs-related mRNAs is still based on miRNA-like methods. In the present study, we used two common miRNA databases to predict the targets of four differential expression tRFs/tiRNAs, and the roles of each tRFs/tiRNAs were inferred by the associated biological functions of the corresponding target mRNAs. We further established a PPI network for predicted target genes of the four tRFs/tiRNAs, and found that PPP2R5A, ADAMTS1, PPARA and POLR2C were the most frequently interacted proteins in target genes of these tRFs/tiRNAs.

Conclusion
In conclusion, our study provided a comprehensive analysis of tRFs/tiRNAs in FD, and the results indicated that tDR-006826, tDR-006049, tDR-001271 and tDR-001276 might play key regulatory roles in FD, which was expected to provide a theoretical basis and data support for clinical treatment. Culture medium was changed every 3 days until 70-80% con uence was achieved.

Methods
RNA extraction and quality control Total RNA from cells was extracted using TRIzol reagent (Takara, Japan) according to the manufacturer's instruction. The integrity and quantity of RNA samples were checked using agarose gel electrophoresis and Nanodrop™ instrument. The optical density 260/280 absorbance ratio of total RNAs was between 1.8 and 2.0.

High-throughput sequencing
Total RNAs were pretreated to remove RNA modi cations. Puri ed RNA was sent to KangChen Bio-tech (Shanghai, China) for small RNA-seq library construction. The libraries are quali ed and absolutely quanti ed using Agilent BioAnalyzer 2100 (Agilent, California, CA). Small RNA sequencing analysis was performed using an Illumina NextSeq 500 system (Illumina, California, USA) according to the instructions of manufacturer.
Analysis of high-throughput sequencing data The differential expression of tRFs/tiRNAs was analyzed using R package edgeR based on the highthroughput sequencing data. PCA, pie plots, Venn plots, hierarchical clustering, scatter plots, and volcano plots were analyzed in the R environment.

Real-time PCR and statistics
Total RNAs were reverse transcribed into cDNA using Bulge-loop TM microRNA (miRNA) qRT-PCR Primer Sets (Ribobio, Guangzhou, China). Quantitative real-time PCR analyses were performed with U6 small nuclear RNA (snRNA) as an internal control, and reactions were detected using an Applied Biosystems 7900HT Fast Real-time PCR system (Applied Biosystems, Gaithersburg, CA, USA). The primers of tRFs were designed by Ribobio. The relative expression fold of tRFs in different samples was calculated using the 2 −ΔΔCt method. All reactions were performed in triplicate.

Target gene prediction of tRFs/tiRNAs
The tRNA database applied to map the reads was based on the Genomic tRNA Database (http://gtrnadb.ucsc.edu/). Because of the miRNA-like functions of tRFs/tiRNAs, TargetScan and miRanda databases were used to predict tRFs/tiRNAs target genes.

GO, KEGG, and PPI network analyses
To determine the biological function of the target genes, GO was performed to analyze the target genes. The KEGG database (http://www.kegg.jp/) was used to identify signi cantly enriched pathways of differential expression tRFs/tiRNAs for pathway analyses. STRING database (https://string-db.org/cgi) was used to construct a PPI network among the predicted target genes with con dence ≥ 0.4 as the threshold condition. The network was constructed by cytoscape and cytoHubba plug-in was used to obtain core genes. Informed consent for all patients was obtained before specimen collection in this study.

Consent for publication
Not applicable Availability of data and materials The datasets used and analysed during the current study are available from the corresponding author on reasonable request.

Competing interests
The authors declare that they have no competing interests.      Prediction of tRFs/tiRNAs target genes. a The potential targets components and prediction target genes of tDR-006826. b The potential targets components and prediction target genes of tDR-006049. c The potential targets components and prediction target genes of tDR-001271. d The potential targets components and prediction target genes of tDR-006826.