Expression profiles of tRFs/tiRNAs in FD BMSCs and normal BMSCs
We used a high‐throughput sequencing technique to detect expression profiles 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 file 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 file 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 profiles 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 specifically in BMSCs, and 17 types of tRFs/tiRNAs specifically 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 significantly downregulated tRFs (tRF-22-8EKSP1852 [tDR‐006826], tRF-18-H9R8B7D2 [tDR‐006049] and tRF-33-86V8WPMN1E8Y0E [tDR‐001271]) and 1 upregulated tRF (tRF-34-87R8WP9I1EWJIQ [tDR‐001276]) in BMSCs and FD BMSCs were measured by qRT-PCR (Fig. 3a). The four tRFs positions on the cloverleaf secondary structure of derived tRNA were shown in Fig. 3b.
Prediction of target genes with bioinformatics tool
To further investigate the function of the differential expression tRFs/tiRNAs in FD BMSCs and BMSCs, we first performed target gene prediction of these tRFs (tDR‐006826, tDR‐006049, tDR‐001271, and tDR‐001276) via TargetScan and miRanda databases and showed the target components in Fig. 4a, b, c, d. tDR‐006826, tDR‐006049, tDR‐001271 and tDR‐001276 were predicted to 107, 78,148, and 275 target transcripts, respectively (Fig. 4 a, b, c, d).
GO enrichment and KEGG pathway analysis
The Gene Ontology (GO) analysis, which including biological process (BP), cellular components (CC), and molecular function (MF), was performed for the functional analysis for target genes. The target genes were most involved in phosphate ion binding (GO: 0042301, tDR‐006826) (Fig. 5a), transcription corepressor activity (GO: 0003714, tDR‐006049) (Fig. 5c), Transcriptional activator activity, RNA polymerase II proximal promoter sequence-specific DNA binding (GO: 0001077, tDR‐001271) (Fig. 5e), and four-way junction DNA binding (GO: 0000400, tDR‐001276) (Fig. 5g). By BP assay, it is found that tDR‐006826 was involved in Chromosome, centromeric region (GO: 0000775) (Fig. 5a), tDR‐006049 in intrinsic component of plasma membrane (GO: 0031226) (Fig. 5c), tDR‐001271 in intracellular membrane-bounded organelle (GO: 0043231) and tDR‐001276 in nucleoplasm part (GO: 0044451) (Fig. 5e, g), and their molecular functions included calcium ion transmembrane import into cytosol (GO: 0097553), regulation of plasma membrane organization (GO: 1903729), endocrine system development (GO: 0035270) and isoprenoid biosynthetic process (GO: 0008299), respectively (Fig. 5a, c, e, g). After mapping the targeted genes to terms in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database, we found that the target genes participated in Th17 cell differentiation (tDR‐006826), FoxO signaling pathway (tDR‐006049), cellular senescence (tDR‐001271), and osteoclast differentiation (tDR‐001276) (Fig. 5b, d, f, h).
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.