Altered Expression of Transfer-RNA-Derived Small RNAs in Human With Rheumatic Heart Disease

Rheumatic heart disease (RHD) remains a severe public health problem in developing countries. Atrial fibrillation (AF) is a medical complication of RHD. Although the understanding of disease pathogenesis has advanced in recent years, the key questions need to be addressed. Transfer RNA–derived small RNAs (tsRNAs) are a novel type of short non-coding RNAs with potential regulatory functions in various physiological and pathological processes. The present study used tsRNAs sequencing to investigate the relationship between RHD and atrial fibrillation (AF). Three paired cardiac papillary muscles were taken from six rheumatic RHD patients with AF (3 cases) or without AF (3 cases) from January 2016 to January 2017 in Xiangya Hospital, Central South University. A total of 219 precisely matched tsRNAs were identified, and 77 tsRNAs (fold change > 2.0 and P < 0.05) were differently changed. Three tsRNAs (AS-tDR-001269, AS-tDR-001363, AS-tDR-006049) were randomly selected and confirmed by qRT-PCR. The results of qRT-PCR were consistent with tsRNAs sequencing, suggesting the tsRNAs sequencing was reliable. Subsequently, we predicted the target mRNAs of the three tsRNAs. Moreover, we verified the functions of tsRNAs targeting mRNAs in vitro. Finally, bioinformatics analysis indicated that the target genes were abundant in regulation of transcription, DNA binding, intracellular. Most of the genes were predicted to interplay with cytokine-cytokine receptor by KEGG analysis. Our findings uncover the pathological process of AF in RHD through tsRNAs sequencing. This research provides a new perspective for future research on elucidating the mechanism of AF in RHD and offers potential new candidates for the treatment and diagnosis.


INTRODUCTION
Rheumatic heart disease (RHD) is a chronic autoimmune valvulitis, resulting from an autoimmune response to a group A streptococcal infection (1,2). RHD remains a neglected disease and is a major cause of morbidity and mortality in many developing countries (2,3). It is currently estimated that 40.5 million individuals around the world live with RHD (4). There were 306,000 deaths due to RHD during 2019 worldwide (4). Atrial fibrillation (AF) is a medical complication of RHD. AF occurs in approximately one-fifth of patients with RHD (5,6). However, the pathogenesis of AF in RHD and the underlying signaling pathways are still poorly understood.
In recent years, a growing number of evidence indicated that the small non-coding RNAs play important roles in the pathophysiological mechanism of AF (7,8). With the advent of multiple high-throughput sequencing technologies, numerous novel classes of small RNA have emerged. New classes of "non-micro-short" RNAs named transfer-RNA-derived small RNAs (tsRNAs, <50 nucleotides) attracted our attention (9)(10)(11)(12)(13)(14)(15)(16)(17). tsRNAs are the second abundant class of small noncoding RNAs (18) and can regulate biological processes, such as proliferation, apoptosis, and epigenetic inheritance (19,20). tsRNAs are abundant small ncRNAs that account for 4-10% of all cellular RNA (21). Generally, tsRNAs are divided into two main types, tRNA-derived fragment (tRF)s and tRNA-derived stress-induced RNA (tiRNA)s, based on their length and cleavage sites. tsRNAs regulate gene expression by directly inhibiting protein synthesis (22) or acting as the guide RNA in a miRNAfashion (14). Multiple innovative investigations have divulged that dysregulated tsRNAs are closely related to human diseases, such as neurological disorders (23), metabolic disorders (24), and cancer (25). Emerging evidence has proved that tsRNAs are detected in the heart. The cardiac pathophysiological conditions could be induced to tsRNAs biogenesis (26). A new study revealed that tsRNAs were potential therapeutic targets to cure myocardial ischemic injury (27). However, the expression of tsRNAs of RHD with AF is never discussed. Therefore, we intended to unveil the potential pathological mechanism of AF in RHD via tsRNAs. In addition, tsRNAs will provide novel approaches in grasping new therapeutic targets and understanding the underlying mechanisms of AF in RHD.
In the present study, we discovered the distinct difference in the expression of tsRNAs between RHD with AF and RHD without AF. Bioinformatics analysis identified potential targets genes and evaluated the putative biological functions. Our findings elucidate the molecular mechanism underlying RHD with AF and advance the knowledge of AF, which is of great clinical significance.

Patients Population
From January 2016 to January 2017, patients with RHD undergoing mitral valve replacement surgery (MVR) were enrolled from Xiangya Hospital, Central South University, Changsha, China. Written informed consents were acquired from all patients, and the study was passed by the ethics committee of the hospital and following the relevant guidelines and regulations. The number of Ethical Review is 201512546. Cardiac papillary muscles were obtained from patients who exhibited clinical characteristics of RHD with AF (n = 3) and without AF (n = 3). AF group: patients had permanent AF (documented arrhythmia more than 6 months) with mitral valve stenosis. The exclusion criteria contained a history of using anti-arrhythmic medications in the past 6 months, myocardial infarction, ischemic cardiomyopathy, heart failure, other types of arrhythmias, chronic hepatic or renal failure, and diabetes. The cardiac papillary muscles were immediately frozen in liquid nitrogen after surgical excision and stored at −80 • C before sequencing.

RNA Extraction
The total RNA was isolated from cardiac papillary muscles. Briefly, the tissues were homogenized with an electric homogenizer after adding TRIzol (Invitrogen life technologies). Then, Chloroform was added and centrifuged at 12,000 × g for 15 min to dissolve the RNA in the aqueous phase. Adding isopropanol to make RNA precipitated, and the resultant RNA pellet was then washed with 75% ethanol and dissolved in RNase-free water. Using NanoDrop ND-1000 identified the quality and concentration of RNA. The total optical densities at a 260/280 nm absorbance ratio of all total RNA samples ranged from 1.8 to 2.0. All RNA solutions were then stored at −80 • C.

tsRNAs Sequencing
The total RNA of the tissues for sequencing was pretreated to remove some RNA modifications. The following experiments were carried out to remove some RNA modifications that may disturb small RNA-sequencing library construction (28,29): 3 ′aminoacyl (charged) deacylation to 3 ′ -OH (hydroxyl group) for 3 ′ adaptor ligation, 3 ′ -cP (2 ′ , 3 ′ -cyclic phosphate) removal to 3 ′ -OH for 3 ′ adaptor ligation, 5 ′ -OH phosphorylation to 5 ′ -P for 5 ′ -adaptor ligation, and N1-methyladenosine and N3methylcytidine demethylation for efficient reverse transcription. All methods were conducted based on the rtStar tRF&tiRNA Pretreatment Kit (Arraystar, USA) protocols. The Shanghai BioChip Company constructed the small RNA library and carried out the Solexa high-throughput sequencing following their standard protocols. Briefly, the total RNA of each sample was sequentially ligated to 3 ′ and 5 ′ small RNA adapters. cDNA was then synthesized and amplified on Illumina's proprietary RT primers and amplification primers. Subsequently, ∼135-160 bp PCR amplified fragments were extracted and purified from the PAGE gel. The purified libraries were qualified with the NanoDrop TM ND-1000 Fluorometer (Thermofisher, ND-1000, German) and validated using the Agilent 2100 bioanalyzer (Agilent, G2938C, USA) to verify the insert size and figure out the molar concentration. Only the library that passed quality control was sequenced on an Illumina NextSeq 500/550 V2 kit (#FC-404-2005, Illumina, San Diego, CA, USA). Sequencing was carried out by running fifty cyclings.

Data Processing and Analysis
Sequencing quality was examined by FastQC software, and trimmed reads (pass Illumina quality filter, trimmed 3 ′ -adaptor bases by cut adapt) were aligned to mature-tRNA and pre-tRNA sequences from GtRNAdb (http://gtrnadb.ucsc.edu/) using NovoAlign software (v2.07.11). Only exactly matched reads were selected as tsRNAs. Moreover, tsRNAs expression levels were calculated and normalized as tag counts per million of total aligned tRNA reads (TPM). The expression profiling and differential expression analysis of tsRNAs were measured by the average TPM. The expression profiling and differential expression of tRNAs were calculated based on fold-change > 2.0 and P < 0.05 normalized TPM (30). Hierarchical clustering and volcano plots were conducted in the differentially expressed tsRNAs in the R environment for statistical computing and graphics.
Small RNA Real-Time Quantitative PCR tsRNAs were reverse transcription with specific primers (   total RNA was isolated from the transfected cells. The tsRNAtargeted genes were then measured by qRT-PCR. The specific primers were listed in Table 2, and the protocols were described as above.

Target Prediction and Bioinformatics
tsRNAs could target mRNA leading to mRNA degradation in a microRNA (miRNA) manner. Here we used two common algorithms to predict tsRNA targets, namely, TargetScan v6.0 (http://www.targetscan.org) and miRanda (http://www. microrna.org) (31,32). The overlapping target genes were applied to further bioinformatics. The biological process of the target genes was conducted by Gene Ontology (GO) annotations and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, through DAVID Bioinformatics Resources 6.8 (https:// david.ncifcrf.gov) (33). Cytoscape software (version 3.7.2, the Cytoscape Consortium, San Diego, CA, USA) was used to construct the network.

Statistical Analysis
Data are presented as mean ± standard error. Two-group difference analysis was used Student's t-test. The limma package in R software (version 4.0.4) was applied to determine the differential expression of tsRNAs and the pheatmap package in R was used to construct heat map. The ggplot2 built the figure of heat map. Fold change > 2.0 and P < 0.05 were considered to indicate a statistically significant difference in sequence analysis.

Expression Profiles of tsRNAs
The clinical characteristics of patients were summarized in Table 3, including age, gender, and color doppler echocardiography. To explore the involvement of small RNAs in RHD patients, the cardiac papillary muscles of the patients were processed for tsRNAs sequencing (tsRNA-seq). The RNA-seq data have been deposited into GEO (GSE185581).
The correlation analysis was based on the TPM counts of each sample. The correlation coefficient is a vital evaluation criterion of the reliability and reasonability of the sample selection (34,35). As shown in Figure 1A, the correlation coefficient of the compared samples in the same group was more than 0.9 and in the different groups was <0.75. That is to say, a distinguishable tsRNAs expression profiling was found among the two groups. A total of 219 tsRNAs were identified (211 in RHD without AF, 118 in RHD with AF) (Figure 1B). The tsRNAs of the tsRNAs' distributed at the length of 16-21 and 31-37 nt. RHD with AF group compared to RHD without AF group, the content of tsRNAs with different lengths was changed ( Figure 1C).

Changes of tsRNAs Expression
To investigate whether tsRNAs types were altered, we estimated the subtype numbers of tsRNA transcripts in both RHD without AF and RHD with AF groups. Over 99% of tsRNAs were originated from mature tsRNAs (tRF-1, tRF-3, tRF-5, i-tRF, tiRNA-3, tiRNA-5). Further analysis showed that most tsRNAs in both groups were tiRNA-5, and the proportion of RHD without AF and RHD with AF was 27.01 and 49.15%, respectively (Figures 2A,B). The pie chart demonstrated that the RHD with AF group mainly increased the expression of tiRNA-5 and decreased the expression of other tsRNAs (Figures 2A,B). In addition, the numbers of tsRNAs derived from the variable anticodon tRNAs are demonstrated in the stacked plots (Figures 2C,D). All results suggested that the types of tsRNAs were different in RHD without AF and RHD with AF groups.

Identification of Related tsRNAs and qRT-PCR Confirmation
We looked at changes in expression for the individual tsRNAs using a standard of fold change > 2.0 and P < 0.05 for significant changes in expression. A total of 77 tsRNAs were differentially expressed in the RHD with AF group compared to RHD without AF group. The volcano plot showed six tsRNAs up-regulated and 71 tsRNAs down-regulated ( Figure 3A). Fourteen tsRNAs (4 up-regulated: AS-tDR-001270, AS-tDR-001297, AS-tDR-001269, and AS-tDR-001289; 10 down-regulated: AS-tDR-000205, AS-tDR-000123, AS-tDR-007294, AS-tDR-007326, AS-tDR-007245, AS-tDR-000102, AS-tDR-006049, AS-tDR-001363, AS-tDR-000886, AS-tDR-000894) were significantly altered in RHD with the AF group compared with RHD without AF group ( Figure 3B, Table 1). We constructed a hierarchical clustering map to examine these differentially expressed tsRNAs. The RHD with the AF group clustered together in one group were primarily distinct from the RHD without the AF group ( Figure 3B). We used qRT-PCR to confirm the expression changes for the three tsRNAs (AS-tDR-001269, AS-tDR-001363, and AS-tDR-006049). The expression level of AS-tDR-001269 (P = 0.0023) was up-regulated with the statistical difference between RHD with the AF group and RHD without the AF group. The expression levels of AS-tDR-001363 (P = 0.0292) and AS-tDR-006049 (P = 0.0076) in RHD with the AF group AS-tDR-006049 (P = 0.0076) were statistically different between RHD with the AF group and RHD without the AF group. Data were present as mean ± SEM (n = 3 for each group). *P < 0.05 represent RHD with AF compared to RHD without AF; **P < 0.01 indicated RHD with AF compared to RHD without AF. qRT-PCR, quantitative real-time PCR; RNA-Seq, RNA sequencing.
were down-regulated (P < 0.05), compared to RHD without AF group ( Figure 3C). The qRT-PCR results of AS-tDR-001269, AS-tDR-001363, and AS-tDR-006049 were consistent with the RNA-seq, indicating that the results had higher reliability.

Prediction of Target Genes of tsRNAs and Validation of Target Genes
A growing number of evidence has revealed that tsRNAs contain some seed sequences that might match the seed regions of mRNA by antisense pairing, regulating the expression level of the target mRNA (16,23,36). Although different algorithms can get possible seed sequences and targets for tsRNAs, each methodology for tsRNA target prediction is referenced to miRNA target predictors (23). Therefore, the sequences of altered three tsRNAs were loaded to TargetScan and miRanda to acquire the targets genes. According to the above theory, miRanda and Targetscan, two algorithms were used for predicting the target genes. In total, 3,123 mRNA targets were predicted simultaneously for the three validated tsRNAs (Supplementary Figure 1). The target genes of AS-tDR-001269, AS-tDR-001363, and AS-tDR-006049 were 1,861, 1,179, and 336, respectively. In the current study, AS-tDR-001363 was reduced in the RHD with the AF group. The qRT-PCR was used to verify the relationship between tsRNAs and their relative mRNAs. The two mRNA genes (TNFRSF1B and CCL5)-the target genes of AS-tDR-001363-were selected. The binding site and seed sequence of these tsRNAs and their target mRNAs are displayed in Figure 4A. To confirm the relationship between the target genes and tsRNAs, we overexpressed AS-tDR-001363 in AC16 cells to identify the corresponding alterations in tsRNA target genes. After transfection with AS-tDR-001363 mimics, the expression of TNFRSF1B and CCL5 were predominantly downregulated ( Figure 4B). The qRT-PCR results of in vitro experiments explain the relationship between tsRNAs and target mRNAs. Therefore, the forecasted targets could be applied to further analysis.

Biological Function Analysis
We performed a bioinformatics analysis of the target genes with a context <−0.4. GO biological processes and KEGG pathway enrichment analysis was executed to explore the functions of 278 target genes by using the DAVID online analysis tool (Figure 5B). GO Table 4). According to KEGG enrichment analysis, cytokine-cytokine receptor interaction (hsa04060; 8 genes) and proteoglycans in cancer (has05205; 7 genes) were significantly detected ( Figure 5A, Table 4).
Identifying the non-coding RNAs changes profiles of serum (37), atrial appendages (38), atrium samples (39), and aortic valve (40) in RHD with AF has been well-documented in recent years. In this article, we chose the cardiac papillary muscle as the experimental sample. Because resection of this tissue can cause no damage to the surrounding structure during MVR surgery. Meanwhile, RHD can cause mitral pathologic change due to thickening of the papillary muscles (41,42). Therefore, the altered non-coding RNAs profiles of papillary muscle might elucidate the pathological process of RHD.
Currently, high-throughput sequencing and bioinformatics analysis, which may put a deep insight into disease occurrence at the molecular level, are frequently used by scientists. Many studies have discovered that non-coding RNAs are dysregulated in RHD with AF (37-40, 43, 44). Previous research focused on the aberrant non-coding RNAs in disease due to their disease-specific expression profiles (38,45). tsRNAs are abundant small non-coding RNA, constituting 4-10% of all cellular RNA (21). They are the fundamental components of the translation machinery. They deliver amino acids to the ribosome to translate the genetic information in an mRNA template into a corresponding polypeptide chain. Although the regulation of tsRNAs is similar to miRNAs regarding the related physiological and pathological processes, the higher stability and expression levels of tsRNAs place them as ideal biomarkers for diagnosing and prognosis in diseases (46). Recently, correlations between dysregulated tsRNAs expression and disease development have been reported (46). One of the principal drivers of the current tsRNA research is the discovery of abundant tsRNAs and tsRNA in mice and humans (47,48). Additionally, their capabilities as a potential biomarker for disease diagnosis and prognosis have been revealed in clinical studies (46). Thus, it is worth identifying the dysregulation of tsRNAs in RHD with AF. In this study, we explored the tsRNAs expression profiles between RHD with AF and RHD without AF. The results indicated the length of tsRNAs was from 16 to 21 and 31 to 37 nucleotides ( Figure 1C). The figures illustrated that the number, expression level, and type of tsRNAs have changed in the AF group (Figures 1, 2). All the results suggested the tsRNAs may be potential candidates for the pathophysiological process of AF.
Based on bioinformatics analysis, prior research demonstrated dilated cardiomyopathy, hypertrophic cardiomyopathy (38), and metabolic pathway (37) are the most important biological process of RHD with AF. In our research, one crucial pathway was enriched: cytokine-cytokine receptor interaction from KEGG pathway analysis. Cytokines and their receptor networks are an essential component of the body's signal transduction system (49). Cytokines can be widely involved in almost all physiological and pathological states of the body, affecting gene expression, cell membrane permeability, biological enzyme activity, and cytoskeletal protein function, leading to various physical effects on cells. In our study, eight target genes were enriched in cytokine-cytokine receptor interaction involving IL18R1, TNFRSF1B, CCR6, IL2RA, FASLG, TNFRSF14, EDAR, CCL5. CCL5 has been shown to orchestrate the recruitment to inflammatory sites of several inflammatory cell subsets, such as monocytes, neutrophils, dendritic cells, and lymphocytes through the binding to CCR1, CCR3, or CCR5 (50). CCL5 is increased in atherosclerosis (51), Myocardial infarction (50), and RHD (52). Treatment with anti-CCL5 mAb exerted cardioprotective effects (50). Nevertheless, no documents have addressed the relationship between CCL5 and AF. In the present research, CCL5 is a target gene of AS-tDR-001363, which is down-regulated in RHD with AF group. We overexpressed AS-tDR-001363 in AC16 cells to identify the corresponding alterations in tsRNA target genes. Although the tsRNAs mimics could not represent the actual tsRNAs, at present, mimics are usually used to explore the impact of tsRNAs on target genes (23,27,53). After transfection with AS-tDR-001363 mimics, the expression of CCL5 was predominantly downregulated ( Figure 4B). Therefore, our study may provide a novel perspective to treat RHD with AF.
Overall, the research firstly shows the altered expression patterns of tsRNAs in RHD with AF. Given the pathogenesis and prognosis of diseases, we chose myocardial papilla to reveal more regulator function of tsRNAs in RHD with AF. What's more, the validation of tsRNAs function is still primarily needed by future researchers and the future studies with larger sample sizes will be needed to confirm our present results.

DATA AVAILABILITY STATEMENT
The raw data of tsRNAs-sequencing in this article will be acquired from GEO (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE185581).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Ethic Committee of the Xiangya Hospital Central South University. The patients/participants provided their written informed consent to participate in this study.