tRNA‑derived small RNAs: A novel class of small RNAs in human hypertrophic scar fibroblasts.

tRNA‑derived small RNAs (tsRNAs) have been shown to play regulatory roles in many physiological and pathological processes. However, their roles in hypertrophic scars remain unclear. The present study investigated differentially expressed tsRNAs in human hypertrophic scar fibroblasts and normal skin fibroblasts via high‑throughput sequencing. Several dysregulated tsRNAs were validated by reverse transcription‑quantitative polymerase chain reaction (RT‑qPCR). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment, target prediction, coexpression networks and competing endogenous RNA (ceRNA) networks were evaluated to discover the principal functions of significantly differentially expressed tsRNAs. In total, 67 differentially expressed tsRNAs were detected, of which 27 were upregulated and 40 downregulated in hypertrophic scar fibroblasts. The GO analysis indicated that the dysregulated tsRNAs are associated with numerous biological functions, including 'nervous system development', 'cell adhesion', 'focal adhesion', 'protein binding', 'angiogenesis' and 'actin binding'. KEGG pathway analysis revealed that the most altered pathways include 'Ras signaling pathway', 'Rap1 signaling pathway' and 'cGMP‑PKG signaling pathway'. The target genes of the differentially expressed tsRNAs participate in several signaling pathways important for scar formation. The results of RT‑qPCR were consistent with those of sequencing. MicroRNA (miR)‑29b‑1‑5p was identified as a target of tsRNA‑23678 and was downregulated in hypertrophic scar fibroblasts, constituting a negative regulatory factor for scar formation. Furthermore, tsRNA‑23761 acted as a ceRNA and bound to miR‑3135b to regulate the expression of miR‑3135b targets, including angiotensin‑converting enzyme. Collectively, these findings reveal that tsRNAs are differentially expressed in human hypertrophic scar fibroblasts, and may contribute to the molecular mechanism and treatment of hypertrophic scars.

Introduction

Scar is a generic term that refers to the histopathological changes in normal skin tissue caused by various types of trauma, and scars are an inevitable product in the process of human trauma repair (1). After skin damage, normal scars are flat, thin and almost invisible. However, factors such as excessive deposition of extracellular matrix (ECM) and excessive proliferation of fibroblasts can lead to the formation of proliferative scars (2-4). A hypertrophic scar has a protruding surface, and is irregularly shaped with an uneven appearance, red color and solid texture, with burning and itching sensations on the skin surface (5). Hypertrophic scars can damage the function and appearance of the skin and impact the patient's psychology and physiology (6-9). Although proliferative scars can be treated with surgery, radiotherapy, steroid injection and other means, these treatments have not been optimized, and the molecular mechanisms of proliferative scars remain unknown (10).

Recent studies have discovered a new class of small RNAs known as tRNA-derived small RNAs (tsRNAs) that are different from typical noncoding RNA. tsRNAs are fragments derived from transfer RNA (tRNA). tsRNAs can be divided into two main types according to the length of the tRNA and the cleavage site: The first is the stress-induced tRNA fragment, which is produced by a specific cut in the 28-36-nucleotide (nt) anticodon ring and is a mature tRNA, known as tRNA-derived stress-induced RNA or tiRNA; the other tsRNA is 14-30 nt in length and known as a tRNA-derived fragment (tRF) (11). tRFs/tRNA-derived RNAs (tDRs) are derived from mature tRNA or precursor tRNA, which can be subdivided as follows: i) tRF-5′, corresponding to the 5′ end of the mature tRNA according to its corresponding position on the tRNA, obtained by D-loop cleavage; ii) tRF-3′, corresponding to the 3′ end of the mature tRNA and containing a CCA tail sequence, obtained by T-loop cleavage; and iii) tRF-1, which is the 3′-tail sequence of the precursor tRNA and contains a poly T sequence at the 3′ end (12-14).

tsRNA is not only a by-product of random tRNA cleavage but also a small noncoding RNA that has crucial roles in numerous specific physiological and pathological processes (11). Additionally, tsRNA has an impact on the function of various types of cells. Studies have demonstrated that tsRNAs are present in a variety of organisms and affect different biological processes, such as the growth and metastasis of cancer cells, tumor inhibition, regulation of gene and protein expression, initiation of viral reverse transcriptase, RNA processing, the DNA damage response and neurodegeneration (15-18). However, the role of tsRNAs in hypertrophic scar fibroblasts has not yet been reported.

Therefore, the present study investigated differentially expressed tsRNAs in human hypertrophic scar fibroblasts and normal skin fibroblasts via high-throughput sequencing. tsRNAs with partial differential expression were confirmed by reverse transcription-quantitative polymerase chain reaction (RT-qPCR). Bioinformatics techniques such as target gene prediction, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment were then used to investigate the role of tsRNA differential expression. Some of the differentially expressed tsRNAs were selected to generate coexpression and competing endogenous RNA (ceRNA) networks. The main purpose of the study was to elucidate the potential molecular mechanism of tsRNAs in proliferative scars and to lay a foundation for the treatment of such scars.

Materials and methods

Preparation of fibroblasts

The present study was approved by the Ethics Committee of the First Affiliated Hospital of Nanchang University (Nanchang, China). Written informed consent was obtained from the adult participants and from the parent of the child participant. Hypertrophic scars and adjacent normal skin tissues were harvested from three female patients (aged 6, 20 and 26 years) who were admitted to the First Affiliated Hospital of Nanchang University from April to October 2018. A typical manifestation of hypertrophic scar is a bulging, erythematous, itching and thickened scar restricted to the site of injury (19,20). Samples were obtained from patients with hypertrophic scars following perineal burns that had not been treated with drugs or radiotherapy prior to sample collection. The samples were washed three times with PBS, and the subcutaneous tissue was removed. The tissue was placed in a Petri dish and digested with 0.25% trypsin-EDTA (Gibco; Thermo Fisher Scientific, Inc.) at 4°C for 10-12 h. The dermal layer of the tissue was retained; the tissue epidermis was completely removed to avoid contamination of the epidermal cells, and the dermis was cut into 1×1-mm pieces with ophthalmic scissors (21). The cut dermal tissue pieces were plated in a Petri dish, and DMEM containing 10% fetal calf serum (Gibco; Thermo Fisher Scientific, Inc.) was added. The cells were incubated at 37°C under 5% CO2, and the medium was changed every 3 days. The cells disintegrated, expanded and fused into sheets covering 80% of the bottom of the culture dish; they were then digested and passaged. Fibroblasts of the third generation were used in the study. The fibroblasts derived from hypertrophic scar tissue were designated as the experimental group, and those derived from normal skin were designated as the control group.

RNA isolation and quality control

Total RNA was extracted using TRIzol reagent according to the manufacturer's protocol (Invitrogen; Thermo Fisher Scientific, Inc.). The quantity of extracted RNA was measured using a NanoDrop ONE spectrophotometer (Thermo Fisher Scientific, Inc.), and the quality of the extracted RNA, specifically, its integrity and purity, was assessed by standard denaturing agarose gel electrophoresis. The optical density ratio at wavelengths of 260 and 280 nm (OD260/280) was measured as an indicator of RNA concentration and was used to determine RNA purity; an OD260/280 ratio of between 1.8 and 2.1 is considered to indicate acceptable RNA quality (22,23).

Library preparation for RNA sequencing (RNA-seq)

The library was created by first adding a 3′ linker sequence and sequentially adding 1 µl SR RT Primer, 4.5 µl nuclease-free H2O for reverse transcription primer hybridization, and the 5′ adaptor sequence. RNA fragments were used to synthesize cDNA strands as follows: Degeneration at 94°C for 15 sec, annealing at 62°C for 30 sec, and extension at 70°C for 15 sec. The library was then examined on an 8% SDS-PAGE Tris-borate-EDTA gel after 1 h at 120 V. The library was constructed using a HiSeq X Reagent kit v2.5 (FC-501-2501; Illumina, Inc.). Library quality was assessed using an Agilent 2200 TapeStation system (Agilent Technologies, Inc.) prior to sequencing. Total concentrations were determined using a Qubit 3.0 Fluorometer (Thermo Fisher Scientific, Inc.). Libraries were sequenced using a HiSeq X Ten high-throughput sequencer (Illumina, Inc.) with 2×150-bp paired-end reads. Purified cDNA libraries were used for cluster generation (cBot Cluster Generator; Illumina, Inc.). The library was sequenced and constructed by Shanghai Personal Biotechnology Co., Ltd. (www.personalbio.cn/).

Data quality control

First, tsRNA sequences that did not match sequences in miRBase (http://www.mirbase.org) were matched to sequences in the GtRNA database version 18.1 (24). The tsRNA sequences that matched to the GtRNA database were then searched against the tRFdb 2012 (25) and tRF MINTbase v2.0 (26) databases. The expressed tsRNAs in the tissue were identified. In the miRNA analysis, the filtered clean reads were compared with the human miRNA database 22.1 (http://www.mirbase.org/) to obtain miRNA expression data. Levels of tsRNA and miRNA gene expression were obtained by RNA-seq. The level of gene expression is closely associated with the following two aspects: i) The level of gene expression and the number of reads in the mapping reference sequence, that is, the effective sequencing amount; ii) the sum of the length of all exons of the gene, as increased gene transcription generates more sequencing fragment results (27). The number of reads of each gene was then obtained. According to the known small RNA data, the differences in small RNAs between the hypertrophic scar fibroblasts and normal skin fibroblasts were determined. Differentially expressed genes were obtained using DEGseq Software (version 1.18.0) (28). The selection criteria were fold change >2.0, P<0.05 and false discovery rate (FDR) <0.05.

Target prediction

TargetScan software (www.targetscan.org/) was used to predict the potential binding targets of mRNAs, and target genes were predicted by searching for conserved 8- and 7-mer sites that matched each miRNA seed region. The miRanda (http://www.microrna.org) and RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/) algorithms were employed to predict mRNAs targeted by differentially expressed tsRNAs (screening criteria: miRanda, score ≥150 and energy <−20; RNAhybrid, energy <−25). The intersection of the two algorithms was taken as the final result of the target gene prediction. The miRanda score is a binding score of the tsRNA and its target gene, with a higher score indicating a more accurate target gene prediction.

GO analysis

The predicted target genes were subjected to GO analysis (www.geneontology.org/), and all involved GO genes were evaluated. The significance level (P-value) of each GO gene was calculated by Fisher's test, the results were corrected using the FDR method of multiple hypothesis testing to identify significant GO genes from differential gene enrichment at P<0.05. The P-values for the enrichment of GO terms were calculated by the default statistical algorithm of the GO analysis database. A smaller P-value indicates a more significant role of a GO term, and terms with P<0.05 were considered statistically significant. Under the same P-value, a larger enrichment score indicates a more specific effect of a gene.

Pathway analysis

KEGG pathway (www.kegg.jp/) enrichment analysis of the differentially expressed tsRNAs was performed. All pathways in which the genes are involved were determined based on the hypergeometric distribution of Fisher's test calculation of each pathway significance level (P-value), the results of multiple hypothesis correction tests and the FDR. Thus, significant pathways embodied by differential genes were identified. The standard for significance screening was P<0.05; biological pathways with P<0.05 were considered significant.

Validation by RT-qPCR

To confirm the RNA-seq data, three tsRNAs were selected for verification by RT-qPCR. First, cDNA was synthesized using an ABI Q6 RT-qPCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.). RT-qPCR was performed in a 12-µl reaction mixture containing 1 ng total RNA, 1 µl RT primer and 2 µl dNTP mix (HyTest Ltd.), 4 µl 5X buffer, 1 µl Protector RNase Inhibitor and 1 µl transcriptase (all from Epicentre; Illumina Inc.). The stem-loop RT reaction was performed at 65°C for 5 min, followed by 42°C for 60 min and 70°C for 5 min. RT-qPCR was performed with the 2X FastStart Universal SYBR-Green Master mix kit (Roche) with the following amplification procedure: 95°C for 10 min, followed by 45 cycles of 95°C for 15 sec and 60°C for 60 sec for fluorescence collection. After the amplification reaction was complete, a melting curve of the RT-qPCR product was established during fluorescence collection (95°C for 10 sec; 60°C for 60 sec; and 95°C for 15 sec) by slowly heating from 60 to 99°C (+0.05°C/sec). The target gene and internal reference gene of each sample were subjected to RT-qPCR. Primers were synthesized by Shanghai Bioligo Biotechnology Co., Ltd., as listed in Table I. The U6 gene was used as an internal reference for primer design. RT-qPCR quantification was performed using the 2−ΔΔCq method (29).

Table I

Primer sequences of tsRNAs used for reverse transcription-quantitative PCR.

Primer ID	Primer sequence (5′-3′)	Length (bp)	Annealing temperature (°C)
U6	F: CGATACAGAGAAGATTAGCATGGC	61	60
	R: AACGCTTCACGAATTTGCGT
tsRNA-23678	F: GGGGGTATAGCTCAGGGGTAGA	72	60
	R: GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACAGTCAAA
tsRNA-23727	F: GGGGGTATAGCTCAGTGGTAGAG	71	60
	R: GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACAGTCAAA
tsRNA-23761	F: GGGGGTGTAGCTCAGTGGTAGAG	71	60
	R: GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGACAGTCAAA

tsRNA, tRNA-derived small RNA; F, forward; R, reverse.

Expression of several tsRNAs and miRNAs in hypertrophic scar tissue and normal skin tissue

Detection of tsRNA-23678, tsRNA-23761, tsRNA-23727, miR 29b-1-5p and miR-3135b gene expression was performed using RT-qPCR. The hypertrophic scar and normal skin tissue were cut into small pieces and ground into fine powder in liquid nitrogen, using a mortar and pestle. Total RNA was extracted using TRIzol reagent. The Rayscript cDNA Synthesis kit (cat. no. GK8030; Shanghai Generay Biotech Co., Ltd.) was used to synthesize cDNA via reverse transcription using the following temperature protocol: 95°C for 10 min, followed by 45 cycles of 95°C for 15 sec and 60°C for 60 sec for fluorescence collection. The specific cDNA fragments were amplified using a RT-qPCR detection system and primers from Shanghai Bioligo Biotechnology Co., Ltd., as described above. Data are expressed as expression relative to the U6 gene. All reactions were analyzed using the 2−ΔΔCq method.

Statistical analysis

Differences in tsRNA levels between hypertrophic scars and adjacent normal skin samples were analysed using a paired Student's t-test The FDR was calculated to correct the P-values. Differential expression of tsRNAs was defined by a fold-change threshold of >2.0 and P<0.05. P<0.05 was considered to indicate statistical significance. Statistical analyses were performed using SPSS version 19 software (IBM Corp).

tsRNA-miRNA coexpression network and ceRNA network analyses

The differentially expressed tsRNAs contain miRNA binding sites. Three differentially expressed tsRNAs were selected and networks containing them were established. High-throughput sequencing results indicated that tsRNA-23678, tsRNA-23727 and tsRNA-23761 were upregulated in the hypertrophic scar fibroblasts. The tsRNAs and miRNAs associated with the same target gene were obtained. The miRNAs with the same expression trends as the tsRNAs were identified, and significant target genes were screened on the basis of their P-values (P<0.05). First, Cytoscape 3.6.1 (www.cytoscape.org/) was used to construct a tsRNA-miRNA coexpression network and a ceRNA network. The purpose of constructing these networks was to discover potential miRNA reaction elements. The overlap of a binding site for both the tsRNA and miRNA on the same mRNA may aid in the prediction of tsRNA-miRNA-mRNA interactions. The intersections of results from the commonly used software miRanda, Targetscan and psRobot (omicslab.genetics.ac.cn/psRobot/) were used to predict tsRNA-miRNA-mRNA interactions.

Results

Extraction of total RNA

Denaturing gel electrophoresis revealed intense and distinct 28S and 18S ribosomal RNA bands. A smaller, slightly diffuse band consisting of low molecular weight RNA (tRNA and 5S ribosomal RNA) was also observed (Fig. S1). As shown in Table II, the extracted RNAs were acceptable for subsequent tsRNA experiments.

Table II

RNA quantification and quality assurance by spectrophotometry.

Group	OD260/280 C ratio	oncentration (ng/µl)	Quantity (ng)	Result
HA	2.10	653.4	13.07	Qualified
HB	2.05	681.6	13.63	Qualified

OD260/280, ratio of optical density at 260 and 280 nm; HA, hyper-trophic scar; HB, normal skin.

Differential tsRNA and miRNA expression profiles

High-throughput sequencing is an efficient method for studying the biological function of RNAs. tsRNA expression profiles in fibroblasts were obtained by RNA-seq. Thousands of transcripts were detected in the hypertrophic scars and normal tissues by RNA-seq. Among them, 67 differentially expressed tsRNAs were detected (Tables III and IV), of which 27 were upregulated and 40 downregulated in hypertrophic scar fibroblasts. The heatmap of intersample association (Fig. 1) showed an evident difference in transcript levels between hypertrophic scar and normal skin fibroblasts. A total of 149 miRNAs were significantly differentially expressed (fold change >2.0, P<0.05). Of these, 120 miRNAs were upregulated and 29 downregulated. The top 15 upregulated and top 15 downregulated significantly differentially expressed miRNAs are shown Table V. These results reveal that the expression profiles of tsRNAs and miRNAs in hypertrophic scar fibroblasts differed from those in normal skin fibroblasts.

Table III

Data for significant differentially expressed upregulated tsRNAs of hypertrophic scar compared with normal skin.

Yingbio ID	Fold change	Database ID	Log2FC	Anticodons	Type
tsRNA-04990	4.2733	tRF-33-LQR47673FEWSD3	2.0953	GlyGCC	5′-half
tsRNA-04998	2.0229	tRF-33-LQ947673FEWSD3	1.0164	GlyCCC	5′-half
tsRNA-06150	2.0243	tRF-33-FQ947673FEWSD3	1.0174	GlyCCC	5′-half
tsRNA-14697	2.0476	tRF-31-P4R8YP9LON4VD	1.0339	GlyCCC	5′-half
tsRNA-14698	4.1770	tRF-32-P4R8YP9LON4V3	2.0624	GlyCCC	5′-half
tsRNA-14699	4.1187	tRF-33-P4R8YP9LON4VDP	2.0421	GlyCCC	5′-half
tsRNA-14700	4.0846	tRF-34-P4R8YP9LON4VHM	2.0302	GlyCCC	5′-half
tsRNA-14701	4.1227	tRF-35-P4R8YP9LON4VN1	2.0435	GlyCCC	tRF-5
tsRNA-14702	4.0472	tRF-36-P4R8YP9LON4VN1B	2.0169	GlyCCC	tRF-5
tsRNA-14703	3.9967	tRF-38-P4R8YP9LON4VN18	1.9988	GlyCCC	tRF-5
tsRNA-14704	3.9166	tRF-40-P4R8YP9LON4VN1EH	1.9696	GlyCCC	tRF-5
tsRNA-14705	3.9591	tRF-43-P4R8YP9LON4VN1EH	1.9851	GlyCCC	tRF-5
tsRNA-14771	2.0016	tRF-37-PNR8YP9LON4V4R1	1.0012	GlyCCC	tRF-5
tsRNA-23678	5.2361	tRF-32-RKVP4PRL5FZU3	2.3885	CysGCA	5′-half
tsRNA-23679	4.3826	tRF-33-RKVP4PRL5FZUDP	2.1317	CysGCA	5′-half
tsRNA-23680	3.9089	tRF-34-RKVP4PRL5FZUHM	1.9667	CysGCA	5′-half
tsRNA-23727	8.8052	tRF-32-RKVP4P9L5FZU3	3.1383	CysGCA	5′-half
tsRNA-23728	7.1609	tRF-33-RKVP4P9L5FZUDP	2.8401	CysGCA	5′-half
tsRNA-23729	7.3594	tRF-33-RKVP4P9L5FZUDP	2.8795	CysGCA	5′-half
tsRNA-23730	7.1105	tRF-34-RKVP4P9L5FZUHM	2.8299	CysGCA	tRF-5
tsRNA-23761	5.9145	tRF-35-RKVP4P9L5FZUNF	2.5642	CysGCA	5′-half
tsRNA-23762	4.3390	tRF-32-RK9P4P9L5FZU3	2.1173	CysGCA	5′-half
tsRNA-23772	3.4717	tRF-33-RK9P4P9L5FZUDP	1.7956	AlaCGC	5′-half
tsRNA-23774	3.5778	tRF-31-RK9P4P9L5HMVE	1.8391	AlaCGC	5′-half
tsRNA-23775	3.9072	tRF-32-RK9P4P9L5HMVQ	1.9661	AlaCGC	5′-half
tsRNA-26551	3.3115	tRF-33-RK9P4P9L5HMV05	1.7274	SerGCT	tRF-1
tsRNA-26578	2.1116	tsRNA-1005	1.0783	GlyTCC	tRF-1

tsRNA, tRNA-derived small RNA; database ID, from the tRFdb database 2012; tRF, tRNA-derived fragment.

Table IV

Data for significant differentially expressed downregulated tsRNAs of hypertrophic scar compared with normal skin.

Yingbio_ID	Fold change	Database_ID	Log2FC	Anticodons
tsRNA-05068	0.3915	tRF-34-LSXMSL73VL4YHE	−1.3527	HisGTG
tsRNA-05069	0.2660	tRF-34-LSXMSL73VL4YHE	−1.9100	HisGTG
tsRNA-05070	0.2673	tRF-35-LSXMSL73VL4YMY	−1.9033	HisGTG
tsRNA-05071	0.2711	tRF-36-LSXMSL73VL4YMYE	−1.8828	HisGTG
tsRNA-06211	0.3743	tRF-33-FSXMSL73VL4YDN	−1.4174	HisGTG
tsRNA-06212	0.2695	tRF-34-FSXMSL73VL4YHE	−1.8916	HisGTG
tsRNA-06213	0.2575	tRF-35-FSXMSL73VL4YMY	−1.9569	HisGTG
tsRNA-06214	0.2562	tRF-36-FSXMSL73VL4YMYE	−1.9642	HisGTG
tsRNA-10278	0.3460	tRF-32-6SXMSL73VL4YK	−1.5309	HisGTG
tsRNA-10279	0.3335	tRF-33-6SXMSL73VL4YDN	−1.5840	HisGTG
tsRNA-10280	0.2544	tRF-34-6SXMSL73VL4YHE	−1.9747	HisGTG
tsRNA-10281	0.2471	tRF-35-6SXMSL73VL4YMY	−2.0162	HisGTG
tsRNA-10282	0.2668	tRF-36-6SXMSL73VL4YMYE	−1.9057	HisGTG
tsRNA-10283	0.2766	tRF-37-6SXMSL73VL4YMYP	−1.8538	HisGTG
tsRNA-15103	0.3519	tRF-31-PW5SVP9N15WV0	−1.5064	HisGTG
tsRNA-15104	0.3098	tRF-32-PW5SVP9N15WVN	−1.6902	HisGTG
tsRNA-15105	0.2369	tRF-33-PW5SVP9N15WV0E	−2.0770	HisGTG
tsRNA-15106	0.2555	tRF-34-PW5SVP9N15WV2P	−1.9682	HisGTG
tsRNA-15107	0.2545	tRF-35-PW5SVP9N15WV7W	−1.9739	HisGTG
tsRNA-15108	0.2731	tRF-36-PW5SVP9N15WV7W0	−1.8720	HisGTG
tsRNA-15109	0.2728	tRF-39-PW5SVP9N15WV7WIV	−1.8737	HisGTG
tsRNA-15110	0.2806	tRF-40-PW5SVP9N15WV7W6S	−1.8331	HisGTG
tsRNA-15111	0.2966	tRF-45-PW5SVP9N15WV7W6SJF	−1.7533	HisGTG
tsRNA-17436	0.3188	tRF-34-94U47P299DZUFJ	−1.6489	GlnTTG
tsRNA-18787	0.2153	tRF-35-VHJYOZ4E8BRY93	−2.2149	ArgTCT
tsRNA-19074	0.4575	tRF-37-V4QP596VZ631QJJ	−1.1279	GluTTC
tsRNA-21604	0.3724	tRF-32-XSXMSL73VL4YK	−1.4247	HisGTG
tsRNA-21605	0.2994	tRF-33-XSXMSL73VL4YDN	−1.7395	HisGTG
tsRNA-21606	0.2578	tRF-34-XSXMSL73VL4YHE	−1.9551	HisGTG
tsRNA-21607	0.2407	tRF-35-XSXMSL73VL4YMY	−2.0546	HisGTG
tsRNA-21608	0.2578	tRF-36-XSXMSL73VL4YMYE	−1.9551	HisGTG
tsRNA-21609	0.2653	tRF-37-XSXMSL73VL4YMYP	−1.9141	HisGTG
tsRNA-21610	0.2803	tRF-40-XSXMSL73VL4YMY91	−1.8348	HisGTG
tsRNA-21611	0.2802	tRF-42-XSXMSL73VL4YMY91M -	1.8350	HisGTG
tsRNA-24180	0.3091	tRF-32-ROD8N0X0JYOYO	−1.6937	TyrGTA
tsRNA-24181	0.2991	tRF-33-ROD8N0X0JYOY0F	−1.7411	TyrGTA
tsRNA-24182	0.2995	tRF-34-ROD8N0X0JYOY26	−1.7389	TyrGTA
tsRNA-24183	0.3123	tRF-36-ROD8N0X0JYOYUED	−1.6787	TyrGTA
tsRNA-24184	0.3030	tRF-37-ROD8N0X0JYOYUE3	−1.7222	TyrGTA
tsRNA-26579	0.4958	tsRNA-1042	−1.0121	ThrAGT

tsRNA, tRNA-derived small RNA; database ID, from the tRFdb database 2012; tRF, tRNA-derived fragment.

Figure 1

Heatmap showing expression profiles of differentially expressed tsRNAs. The map is based on the expression values of all expressed tsRNAs detected by RNA sequencing with fold change >2.0, P<0.05 and a false discovery rate <0.05. The expression values are indicated by the color scale; red indicates high relative expression and green indicates low relative expression. tsRNA, tRNA-derived small RNA.

Table V

Data of significant differentially expressed miRNAs of hypertrophic scar compared with normal skin (log2FC >1.0, P<0.05).

AccID	Fold change	Log2FC	Regulation
miR-29b-1-5p	0.3214	−1.6373	Down
miR-1275	0.4445	−1.7057	Down
miR-148b-5p	0.3617	−1.4671	Down
miR-151b	0.3072	−1.7025	Down
miR-18a-5p	0.3044	−1.7156	Down
miR-200a-3p	0.4043	−1.3064	Down
miR-215-5p	0.2924	−1.7736	Down
hsa-miR-222-5p	0.4006	−1.3196	Down
miR-25-5p	0.4175	−1.2601	Down
miR-324-5p	0.4280	−1.2242	Down
miR-335-5p	0.2566	−1.9623	Down
miR-491-5p	0.3522	−1.5054	Down
miR-590-3p	0.3644	−1.4563	Down
miR-664a-3p	0.4117	−1.2802	Down
miR-7-1-3p	0.3386	−1.5621	Down
miR-3135b	2.5475	1.3491	Up
miR-10a-5p	39.3914	5.2998	Up
miR-10b-5p	57.0105	5.8331	Up
miR-129-5p	2.4527	1.2944	Up
miR-133a-3p	4.0421	2.0151	Up
miR-137	2.8405	1.5061	Up
miR-1-3p	4.1297	2.0460	Up
miR-140-3p	3.5733	1.8372	Up
miR-140-5p	2.3318	1.2215	Up
miR-142-3p	10.0859	3.3342	Up
miR-142-5p	7.4310	2.8935	Up
miR-143-3p	4.9665	2.3122	Up
miR-143-5p	8.9891	3.1681	Up
miR-145-3p	4.5983	2.2011	Up
miR-200b-3p	3.2765	1.7121	Up

miR, microRNA; AccID, gene name according to miRbase.

To explore the potential functions of tsRNAs in hypertrophic scars, GO analysis of differentially expressed tsRNAs was performed. The important GO terms and the genes involved were obtained via the gene function analysis of target genes. Fig. 2 shows the degree of enrichment of the top 15 genes identified, with the enrichment degree listed in descending order. In the −log10 calculation of the P-value, a smaller P-value corresponds to a greater −log10 P-value. The GO project covers three main domains, namely biological processes (BP), molecular function (MF) and cellular component (CC). In total, 5,891 genes showed differences in the GO analysis. The data indicate that genes in the BP category are mainly associated with the following: 'Homophilic cell adhesion' via plasma membrane adhesion molecules; 'cell adhesion'; 'positive regulation of transcription,' DNA-templated; 'epidermal growth factor receptor'; and 'negative regulation of cell migration'. Genes in the MF category are mainly associated with 'metal ion binding', 'protein binding' and 'calcium ion binding'; genes in the CC category are mainly associated with 'cytoplasm', 'cell junction', 'plasma membrane' and 'cytosol'. In summary, these data indicate that differentially expressed target gene tsRNAs may be involved in cell division and proliferation.

Figure 2

GO analysis of differentially expressed tsRNAs. The top 15 enriched terms covering (A) BP, (B) MF and (C) CC are presented. GO, Gene Ontology; tsRNAs, tRNA-derived small RNA. BP, biological processes; MF, molecular functions; CC, cellular components.

KEGG pathway analysis

KEGG pathway analysis indicated that the dysregulated tsRNAs are involved in numerous signaling pathways, with 283 KEGG pathways associated with the differentially expressed tsRNAs. A histogram showing the -log10 P-values of KEGG pathways is shown in Fig. 3. The −log10 (P-value) denotes an enrichment score for the significance of pathway correlations. The histogram shows the top 15 enriched pathways, and the most significantly altered pathways included the 'Ras signaling pathway', 'Rap1 signaling pathway' and the 'cGMP-PKG signaling pathway'. These genes may serve important roles in cell proliferation, differentiation, metastasis and apoptosis.

Figure 3

KEGG pathway analysis of differentially expressed tsRNAs. The top 15 enriched pathways of tsRNAs in the KEGG analysis are presented. -Log10 represents the -log10 of the P-value, and a smaller P-value is indicated by a larger -log10 P-value. KEGG, Kyoto Encyclopedia of Genes and Genomes; tsRNAs, tRNA-derived small RNAs.

mRNA-target prediction of tsRNA

The results showed that one tsRNA can target multiple mRNAs and that one mRNA can also be a target gene for multiple tsRNAs. As shown in Table VI, the target genes of tsRNA-23678 were predicted to be FAS, SMAD2, ABCA1, ABCCB, ABHDS, ABL2, AIPL1 and BB59. Predicted target genes of tsRNA-23761 were COL1A, A1CF, ABI3, BEST1, ABCA1, ACAP3, ACE and B3GNT7, and those of tsRNA-23727 were SMAD2, TGFBR1, ABCC8, ABHD5, ABL2, BEAN1, AIPH and BAALC. SMAD2 was a target gene for both tsRNA-23678 and tsRNA-23727. In the target gene prediction, the miRanda scores of tsRNA23727, tsRNA23761 and tsRNA23678 were 165, 150 and 157, respectively.

Table VI

Target genes for the three differentially expressed tsRNAs.

Yingbio_ID	Target gene
tsRNA-23678	FAS, SMAD2, ABCA1, ABCCB, ABHDS, ABL2, AIPL1, BB59
tsRNA-23727	SMAD2, TGFBR1, ABCC8, ABHD5, ABL2, BEAN1, AIPH, BAALC
tsRNA-23761	COL1A1, A1CF, ABI3, BEST1, ABCA1, ACAP3, ACE, B3GNT7

tsRNA, tRNA-derived small RNAs.

Validation with RT-qPCR

The differentially expressed tsRNA-23678, tsRNA-23727 and tsRNA-23761 were chosen for further analysis to verify the RNA-seq results. Amplification and dissociation curves for these tsRNAs were generated (Fig. 4), and the 2−ΔΔCq values of the tsRNAs were calculated according to the relative quantitative method. The 2−ΔΔCq analysis revealed upregulation of tsRNA-23678 (1.9452-fold), tsRNA-23727 (1.6684-fold) and tsRNA-23761 (1.6997-fold) in hypertrophic scar, which is consistent with the RNA-seq results (Fig. 5; Table VII).

Figure 4

Amplification plots and dissociation curves for selected tsRNAs. Amplification plots for (A) tsRNA-23678, (B) tsRNA-23727, (C) tsRNA-23761. In the amplification curves, the x-axis represents the cycle number, and the y-axis represents the real-time fluorescence signal intensity of the corresponding cycle number. In the dissociation curves, the x-axis represents the temperature of the RT-qPCR products, and the y-axis represents the real-time fluorescence-signal-intensity change rate with increasing temperature. Differently colored curves correspond to different RT-qPCRs. tsRNAs, tRNA-derived small RNA; RT-qPCR, reverse transcription-quantitative PCR. Amplification plots and dissociation curves for selected tsRNAs. Amplification plots for (D) U6. Dissociation curves for (E) tsRNA-23678, (F) tsRNA-23727. In the amplification curves, the x-axis represents the cycle number, and the y-axis represents the real-time fluorescence signal intensity of the corresponding cycle number. In the dissociation curves, the x-axis represents the temperature of the RT-qPCR products, and the y-axis represents the real-time fluorescence-signal-intensity change rate with increasing temperature. Differently colored curves correspond to different RT-qPCRs. tsRNAs, tRNA-derived small RNA; RT-qPCR, reverse transcription-quantitative PCR. Amplification plots and dissociation curves for selected tsRNAs. Amplification plots for (G) tsRNA-23761 and (H) U6. In the amplification curves, the x-axis represents the cycle number, and the y-axis represents the real-time fluorescence signal intensity of the corresponding cycle number. In the dissociation curves, the x-axis represents the temperature of the RT-qPCR products, and the y-axis represents the real-time fluorescence-signal-intensity change rate with increasing temperature. Differently colored curves correspond to different RT-qPCRs. tsRNAs, tRNA-derived small RNA; RT-qPCR, reverse transcription-quantitative PCR.

Figure 5

Verification of the expression of significantly dysregulated tsRNAs by reverse transcription-quantitative PCR. Relative quantitative levels of three tsRNAs were detected in HA and HB groups. Data are presented as the mean ± SD (n=3); *P<0.05 as indicated. tsRNAs, tRNA-derived small RNA; HA, hypertrophic scar tissues; HB, normal tissues.

Table VII

Comparison of RT-qPCR data with RNA-seq data for tsRNAs.

Method	HA/HB ratio
	tsRNA-23678	tsRNA-23727	tsRNA-23761
RT-qPCR	1.9452	1.6684	1.6997
RNA-seq	5.2361	8.8052	5.9145

RT-qPCR, reverse transcription-quantitative PCR; RNA-seq, RNA sequencing; tsRNA, tRNA-derived small RNA; HA, hypertrophic scar; HB, normal skin.

Expression of the several tsRNAs and miRNAs between hypertrophic scar tissue and normal skin tissue

In order to verify the aforementioned results, RT-qPCR assays were performed on selected tsRNAs (tsRNA-23678, tsRNA-23727 and tsRNA-23761) and miRNAs (miR-3135b and miR-29b-1-5p). Compared with normal skin tissue, the expression levels of tsRNA-23678, tsRNA-23727, tsRNA-23761 and miR-3135b were significantly upregulated in hypertrophic scar tissue, whereas the expression of miR-29b-1-5p was significantly downregulated in hypertrophic scar tissue (Fig. 6). These results demonstrate that the changes observed in tissue were consistent with the bioinformatics results.

Figure 6

Reverse transcription-quantitative PCR analysis for several tsRNAs and miRNAs in hypertrophic scar tissue and normal skin tissue, to verify the in vitro experimental results. Data are presented as the mean ± SD (n=3); *P<0.05 as indicated. tsRNAs, tRNA-derived small RNA; miRNA, microRNA; HA, hypertrophic scar tissues; HB, normal tissues.

Construction of coexpression networks

The functions of most tsRNAs are not currently annotated. The functional prediction of tsRNAs is based on the annotation of coexpressed miRNAs, and three differentially expressed tsRNAs in fibroblasts were chosen in the present study according to the degree of correlation. The coexpression network (Fig. 7) showed that one tsRNA might be associated with one or more miRNAs. A total of 40 miRNAs were associated with the three tsRNAs. Furthermore, the coexpression networks indicated these tsRNAs to be involved in a number of biological processes, including cell adhesion, proliferation, differentiation and metastasis. The network analysis also demonstrated that tsRNA-23678 is associated with miR-29b-1-5p, miR-222-3p and miR-423-5p, which had the same trends in expression. This finding aids in the identification of regulatory relationships between tsRNAs and miRNAs in hypertrophic scars.

Figure 7

Coexpression networks of three significantly dysregulated tsRNAs with their associated miRNAs in hypertrophic scars. Red indicates upregulated genes, and green indicates downregulated genes. tsRNAs, tRNA-derived small RNA; miRNA, microRNA.

Construction of ceRNA networks

The ceRNA network hypothesis provides a new mechanism for tsRNA-miRNA-mRNA interactions. miRNAs are known to cause gene silencing by binding to mRNA, and tsRNAs may regulate gene expression by competitively binding to miRNAs. Thus, tsRNAs can be considered as ceRNAs. According to the ceRNA hypothesis, numerous non-coding RNAs may function as ceRNAs, which compete for the same microRNA response elements (MREs) and regulate each other (21). The analysis of ceRNA interactions aids in the functional characterization of such noncoding transcripts. A tsRNA-miRNA-mRNA network associated with hypertrophic scars was established in the present study using high-throughput sequencing data (Fig. 8). In this network, tsRNA-23761 was positively associated with miR-3135b. Furthermore, the network indicates that tsRNA-23761 is a ceRNA of miR-3135b that targets angiotensin-converting enzyme (ACE) and PYGO2, and tsRNA-23678 is a ceRNA of miR-133a-3p that targets FXR2 and PNAGS.

Figure 8

tsRNA-miRNA-mRNA network in hypertrophic scars, based on tsRNA-miRNA and mRNA-miRNA interactions. In this network, tsRNA aRNA; miRNA, microRNA.

Discussion

tsRNA, which is derived from tRNA, is a newly discovered class of small molecular RNAs that are produced by the cleavage of the tRNA ring by Dicer or angiogenin enzymes (11,30). tsRNA can be classified into five different types: 5′- and 3′-tRNA fragments (tRFs); 5′-and 3′-halves; and 3′U tRFs (31,32). There is growing evidence that that tsRNAs are associated with the development of tumors, cell proliferation and viral replication (16), regulation of cell viability (33), inhibition of protein translation (34,35), regulation of cancer progression (36), offspring metabolism (37,38) and numerous other processes. It has also been reported that tsRNA may have a regulatory function similar to that of miRNA, which can act like a sponge to regulate mRNA stability and participate in gene transcription and translation (33). However, the role of tsRNA in hypertrophic scars has not yet been reported. Hypertrophic scar is a fibrotic disorder, mainly due to the response of the body to injury, and caused by the excessive proliferation of fibroblasts and excessive production of ECM (39). Understanding the relationship between hypertrophic scars and tsRNA may help to elucidate the pathogenesis and pathophysiology of hypertrophic scars and identify potential therapeutic targets for their clinical treatment.

In the present study, high-throughput sequencing was first applied to evaluate the differentially expressed profiles of tsRNAs and miRNAs between hypertrophic scar fibroblasts and normal skin fibroblasts. Differentially expressed tsRNAs and miRNAs were identified, including 27 upregulated tsRNAs, such as tsRNA-23678, tsRNA-23727, tsRNA-23761, tsRNA-04998, tsRNA-06150 and tsRNA-14697, and 40 downregulated tsRNAs, such as tsRNA-05068, tsRNA-05070, tsRNA-06213, tsRNA-06214, tsRNA-10278 and tsRNA-10279. In addition, 120 upregulated miRNAs, including miR-3135b, miR-10a-5p, miR-10b-5p, miR-129-5p and miR-133a-3p, and 29 downregulated miRNAs, including miR-29b-1-5p, miR-1275, miR-148b-5p, miR-151b and miR-18a-5p were identified. Furthermore, GO functional enrichment and KEGG pathway enrichment analyses were performed. These differentially expressed tsRNAs were found to be enriched in a number of important biological processes, including 'nervous system development', 'cell adhesion', 'focal adhesion', 'protein binding', 'angiogenesis' and 'actin binding'. Some biological signaling pathways, such as the 'Ras signaling pathway', the 'Rap1 signaling pathway' and the 'cGMP-PKG signaling pathway' were also significantly enriched. These pathways are associated with the formation of new tissue, remodeling and inflammation, which are important processes in wound healing (40-43). Subsequently, tissue experiments were performed to detect the expression of several tsRNAs (tsRNA-23678, tsRNA-23727 and tsRNA-23761) and miRNAs (miR-3135b and miR-29b-1-5p). The present study revealed that tsRNAs-23678, tsRNA-23727, tsRNA-23761 and miR-3135b were upregulated while miR-29b-1-5p was downregulated in hypertrophic scar compared with normal skin tissue, which is consistent with the in vitro cell experiments. Our data showed that tsRNAs-23678, tsRNA-23727 and tsRNA-23761 belong to the 5′-tsRNA class of tsRNAs. Previous studies have found that 5′-tsRNAs can regulate the differentiation of stem cells and are involved in translational repression (44). Subsequently, the possible roles of tsRNAs-23678 and tsRNA-23761 in hypertrophic scars were revealed by target prediction, coexpression network and ceRNA analyses.

The results showed a significant upregulation of tsRNA-23761, which targets collagen type I alpha 1 chain (COL1A1). COL1A1 encodes the pro-a1 chains of type I collagen, the triple helix of which comprises two a1 chains and one a2 chain. COL1 is a fibril-forming collagen that is found in most connective tissues and is abundant in tendons and dermis (45). Previous studies have confirmed that COL1 is the main structural element of the ECM and that its dysregulated accumulation leads to scarring (45-47). COL1 serves a critical role in the development and progression of hypertrophic scars, and its expression is increased in hypertrophic scars (46,47). The present study revealed that tsRNA-23761 and COL1A1 exhibit common expression trends in hypertrophic scarring. Overall, the dysregulation of tsRNA-23761 may contribute to the formation of hypertrophic scars.

In addition, TargetScan was used to predict potential targets of tsRNA-23727, which included SMAD2. A previous study has reported that SMAD2 gene expression in hyper-trophic scar fibroblasts is significantly upregulated (48), and the present study revealed that the expression of SMAD2 was positively associated with that of tsRNA-23727. The SMAD protein family is homologous to the gene products of the Drosophila mothers against decapentaplegic (Mad) and the Caenorhabditis elegans Sma genes (49). SMAD proteins are signal transducers and transcriptional modulators that mediate multiple signaling pathways (50). In a previous study, enhancement of SMAD expression promoted transforming growth factor (TGF)-β1 secretion, and TGF-β1 induced α-smooth muscle actin (α-SMA) expression, collagen synthesis and the formation of hypertrophic scars (51). SMAD2 regulates numerous cellular processes, including proliferation, survival, apoptosis and dormancy (52-55). Altogether, the results of the present study suggest that increased expression of tsRNA-23727 promotes the development of hypertrophic scars by targeting SMAD2. This finding provides suggests that the function of tsRNAs in hypertrophic scar fibroblasts may involve the regulation of key mRNAs.

Significantly differentially expressed tsRNAs, such as tsRNA-23678, tsRNA-23761 and tsRNA-23727, in hypertrophic scars were selected for construction of an miRNA-tsRNA coexpression network in the present study according to the degree of correlation. Functional prediction of tsRNAs can be based on annotation of the function of coexpressed miRNAs. In the present study, analyses using TargetScan and miRanda databases were performed to identify miRNAs with complementarity to the selected tsRNAs. tsRNA-23678 was complementary to multiple miRNAs, including miR-1268b, miR-1303, miR-1275, miR-17-5p, miR-15b-5p, miR-29b-1-5p, miR-222-3p, miR-423-5p, miR-491-5p, miR-423-5p, miR-3135b, miR-143-5p and miR-133a-3p. These miRNAs are implicated in a number of biological processes, such as cell proliferation, the cell cycle, apoptosis and angiogenesis. Some important related miRNAs are known to be associated with hypertrophic scars, such as miR-29b. Indeed, a previous study reported that miR-29b has anti-fibrosis activity, and demonstrated that the overexpression of miR-29b significantly reduced the expression levels of COL1A1 and α-SMA, inhibited cell proliferation and induced apoptosis in hypertrophic scar fibroblasts compared with normal fibroblasts (56). The miR-29b mimic remlarsen may be an effective treatment for preventing fibrotic scars (hypertrophic scars or keloids) from forming, or for preventing skin fibrosis such as scleroderma (57). miR-29 can directly inhibit YY1 and TGF-β3, improve skeletal muscle atrophy and reduce renal fibrosis by the lowering levels of YY1 and TGF-β pathway proteins (58). Consistent with this, Bi et al (45) reported that the miR-29 family plays a role in hypertrophic scars by regulating the translation of ECM mRNA. In addition, the data from the present study demonstrated that tsRNA-23678 was upregulated in hypertrophic scar compared with normal skin, while the expression of miR-29b-1-5p, a member of the miR-29b family, was markedly lower, which indicated opposing effects. The coexpression network of miRNA-tsRNA indicated that tsRNA-23678 was complementary to miR-29b-1-5p. Functional prediction of tsRNA-23678 can be based on annotation of the function of the coexpressed miR-29b-1-5p. These findings suggest that tsRNA-23678 may be an important regulator of hypertrophic scars. A tsRNA-miRNA-mRNA expression network based on RNA-seq data was constructed to further investigate the mechanisms underlying tsRNAs. This network indicated the predicted tsRNA-miRNA-mRNA associations in hypertrophic scars and included three tsRNAs, six miRNAs and 40 mRNAs. Recent studies have shown that tsRNAs can be used directly as miRNAs to identify target mRNA, both acting on mRNA and inhibiting target genes through translation (59,60). The present study indicated that tsRNA acts as a ceRNA. It is known that miRNA can lead to gene silencing by binding to mRNA, whereas ceRNA regulates gene expression by competing with miRNA (61). The ceRNA hypothesis indicates that RNA molecules harboring MREs can communicate with each other by competing with miRNA. Understanding this novel RNA crosstalk may provide insight into the function of tsRNAs (62). The ceRNA network in the present study showed that tsRNA-23761 may act as a competing endogenous RNA that binds miR-3135b to regulate the expression of miR targets, including ACE. The present high-throughput sequencing analysis demonstrated tsRNA-23761 to be upregulated in hypertrophic scars; associated miR-3135b and ACE were also upregulated, indicating a similar trend. The ACE gene, located in the 17q23.3 position of chromosome 17, encodes angiotensin I-converting enzyme. The full-length ACE gene comprises 44,769 bp, and the full-length ACE mRNA comprises 4,195 nt; the ACE gene harbors a total of 48 exons and encodes 1,306 amino acid residues. Angiotensin II activation by ACE is a significant mediator in wound healing and collagen production. Previous studies have shown that the expression of ACE in scars is significantly higher than that in normal skin, inducing the conversion of angiotensin I to angiotensin II, leading to excessive deposition of ECM and fibrosis, and regulating the activity and proliferation of fibroblasts in proliferative scars (63). ACE is an important regulatory factor for the renin angiotensin aldosterone system and the bradykinin system, which affect many physiological functions of the human body (64). tsRNA-23761 competes with the miRNA pool to regulate ACE expression. It is speculated that tsRNA-23761 may also affect the activity and proliferation of fibroblasts in hypertrophic scars via the regulation of certain mechanisms of ACE (63-65).

In summary, with the development of bioinformatics technology, it has been discovered that tsRNAs serve important roles in the occurrence and development of disease. However, the functions of the majority of tsRNAs have not been annotated to date. In the present study, tsRNAs that were differentially expressed in fibroblasts and tissue derived from hypertrophic scar tissue compared with normal skin tissue were identified. The data indicate that tsRNA-23678, tsRNA-23727 and tsRNA-23761 were significantly upregulated in fibroblasts and tissue derived from hypertrophic scar compared with normal skin tissue. tsRNA-23678 may serve critical functions through interaction with miR-29b-1-5p to modulate hypertrophic scars. In addition, tsRNA-23761 may interact with mRNA, such as ACE, to affect the activity and proliferation of fibroblasts in hypertrophic scars. As hypertrophic scar development is a dynamic process that is affected by many factors, the results of the present study are preliminary and limited. Further in vivo experiments are required to confirmed the results. The significantly differentially expressed tsRNAs identified in the present study may contribute to our understanding of the potential mechanisms of hypertrophic scars, as well as therapeutic strategies.