Profiling of tRNA Halves and YRNA Fragments in Serum and Tissue From Oral Squamous Cell Carcinoma Patients Identify Key Role of 5′ tRNA-Val-CAC-2-1 Half

Oral squamous cell carcinoma (OSCC) is the most common type of head and neck cancer and, as indicated by The Oral Cancer Foundation, kills at an alarming rate of roughly one person per hour. With this study, we aimed at better understanding disease mechanisms and identifying minimally invasive disease biomarkers by profiling novel small non-coding RNAs (specifically, tRNA halves and YRNA fragments) in both serum and tumor tissue from humans. Small RNA-Sequencing identified multiple 5′ tRNA halves and 5′ YRNA fragments that displayed significant differential expression levels in circulation and/or tumor tissue, as compared to control counterparts. In addition, by implementing a modification of weighted gene coexpression network analysis, we identified an upregulated genetic module comprised of 5′ tRNA halves and miRNAs (miRNAs were described in previous study using the same samples) with significant association with the cancer trait. By consequently implementing miRNA-overtargeting network analysis, the biological function of the module (and by “guilt by association,” the function of the 5′ tRNA-Val-CAC-2-1 half) was found to involve the transcriptional targeting of specific genes involved in the negative regulation of the G1/S transition of the mitotic cell cycle. These findings suggest that 5′ tRNA-Val-CAC-2-1 half (reduced in serum of OSCC patients and elevated in the tumor tissue) could potentially serve as an OSCC circulating biomarker and/or target for novel anticancer therapies. To our knowledge, this is the first time that the specific molecular function of a 5′-tRNA half is specifically pinpointed in OSCC.

Introduction

Head and neck squamous cell carcinomas (HNSCC) are a type of malignant tumors located in the epithelium covering the upper aerodigestive tract (1). HNSCC most commonly occurs at the oral cavity (OSCC), is most prevalent in older patients (2), and is strongly associated with tobacco and alcohol consumption (3). HPV infections, on the other hand, account for up to 25% of all HNSCC cases and up to 60% of the oropharyngeal carcinoma subset (4–6). Accumulating evidence has linked the presence of specific non-coding RNAs to cancer diagnosis and prognosis (7–11) and underscored their role as potential regulators of tumorigenesis (12–14). MicroRNAs (miRNAs) are potent regulators of gene expression and are currently the most studied sncRNAs (15). Circulating miRNAs found in a variety of bodily fluids have been demonstrated to have excellent potential as biomarkers of diseases including cancer (11, 16–18) and of physiological conditions such as aging, calorie restriction and dwarfism (19–22). Tumor-derived miRNAs can regulate the expression of oncogenes and tumor suppressor genes in the local microenvironment (7, 23) and by traveling in the circulation as complexes with other proteins and lipids or inside exosomes, can also contribute to priming distant tissues for metastasis (24–27).

Small RNA sequencing used in our previous and recent work revealed significant association between pathologic and physiologic states and changes in the expression of sncRNAs derived from tRNAs and YRNAs (17, 19, 21, 28–33). Additional evidence supports an active role for these sncRNAs in cell function, disease development and progression (34–38). Importantly, we have previously reported that serum levels of 5′ tRNA halves and 5′ YRNA fragments are altered in patients with breast cancer and head and neck squamous cell carcinoma (17, 31). More recently, we found that levels of 5′ tRNA halves change in seminal fluid and tumor tissues from patients with prostate cancer (33). We also previously presented evidence that aging, calorie restriction and dwarfism altered the expression of circulating 5′ tRNA halves in addition to miRNAs (20, 29, 30). In the present study, we report that serum and tumor tissue levels of 5′ tRNA halves and 5′ YRNA fragments change in association with an OSCC diagnosis. Additionally, we suggest novel cancer-related functions for the 5′ tRNA-Val-CAC-2-1 half.

Results

tRNA Halves and YRNA Fragments Detected in the Human Circulation and Tissue Are Mostly Derived From the 5′-End of the Respective Full-Length Molecule

To investigate sncRNAs as potential circulating biomarkers of OSCC, we used small RNA-Seq to measure differential expression of tRNA- and YRNA-derived small RNAs in serum samples collected from subjects without known cancers and patients with OSCC. To reveal whether tRNA- and YRNA-derived small RNAs are also differentially expressed in solid OSCC tissue besides serum, we analyzed a small RNA-Seq generated from tumor and matched normal tissues. We have previously found in our studies, that circulating tRNA- and YRNA-derived small RNAs are 30–34 nucleotides long and almost all of the sequencing reads map to the 5′ end of tRNA or YRNA genes (17, 28, 29, 31, 33). Once again, we find in the present study that 99.97% of the reads with length 30–34 nt derived for tRNA align with the 5′ end of tRNA genes. Similarly, 99.93% of the reads with length 30–34 nt derived for YRNAs align with the 5′ end of YRNA genes. Hence, only reads that are 30–34 nucleotides long and mapping to the 5′ end of tRNA or YRNA genes were used for differential expression analysis in this study.

Levels of 5′ tRNA Halves Decreased in Serum of HNSCC Cancer Patients

We compared the levels of 5′ tRNA halves in serum from subjects with and without OSCC using the edgeR exact test which is applicable to experiments with a single factor (39). We first assessed the biological coefficient of variation (BCV) between samples and found that the cancer samples are markedly separated from the normal samples (Figure 1A, BCV distance 1). Moreover, the normal samples appear more homogeneous than the cancer samples (Figure 1A, BCV distance 2). The differential expression analysis revealed that OSCC decreased the serum levels of 22 types of 5′ tRNA halves with average read counts-per-million (CPM) > 1,770, fold decrease > 5, and FDR < 0.01% (Table 1). The decreased 5′ tRNA halves are derived from five isodecoders of tRNA-Arg, -Glu, -His, -Lys, and -Val. We then examined the tissue expression levels of these 22 types of 5′ tRNA halves and found that they are substantially abundant in tumor and normal tissues (Table 1; see CPM under “Tumor tissues”). However, only one of the 22 types of 5′ tRNA halves (i.e., the 5′ tRNA half derived from the tRNA-Val-CAC-2-1 gene) significantly (FDR < 5%) changed expression in both serum and tissue (Table 1, Figures 1E,F).

FIGURE 1

Figure 1. Small non-coding RNA-Seq analysis in OSCC. (A–D) Multi-dimensional scaling (MDS) analysis of the expression levels of 5′ tRNA-halves and 5′ YRNA fragments. The MDS analysis generates distances that represent the biological coefficient of variation (BCV) between samples. MDS plots show serum expression differences of 5′ tRNA-halves between patients with and without OSCC (A) and between tumor and adjacent normal tissues (B). Also, MDS plots show serum expression differences of 5′ YRNA fragments between patients with and without OSCC (C) and between tumor and adjacent normal tissues (D). Dimension 1 of the MDS plot depicts the cancer effect on the expression levels of 5′ tRNA-halves or 5′ YRNA fragments, while dimension 2 represents the homogeneity between biological replicates. Samples N1–N5 represent Normal while samples C1–C5 represent Cancer. (E,F) Differential expression of a 5′ tRNA-half derived from the tRNA-Val-CAC-2-1 gene in serum from patients with and without OSCC (E), and in solid OSCC tumor relatively to adjacent normal tissue (F). The UCSC genome browser screenshots illustrate the alignment of reads to the tRNA-Val-CAC-2-1 gene. (E) The alignment (number of reads, y-axis) shows that the numbers of reads mapping to the 5′ end of tRNA-Val-CAC-2-1 gene are significantly lower in serum from patients with (red) than without (blue) OSCC. (F) In contrast, the alignment shows that the numbers of reads mapping to the 5′ end of tRNA-Val-CAC-2-1 gene are significantly higher in solid OSCC (red) relatively to adjacent normal tissue (blue). Shown at the bottom is the tRNA-Val-CAC-2-1 gene annotation from the tRNA genes track “Transfer RNA Genes Identified with tRNAscan-SE” associated with the human GRCh38/hg38 genome.

TABLE 1

Table 1. Twenty-two 5′ tRNA halves significantly (FDR < 5%) decreased in serum of patients with OSCC.

Levels of 5′ tRNA Halves Increased in OSCC Tissues Relatively to Adjacent Healthy Tissues

Since most of the 5′ tRNA halves (21 out of 22) that decreased expression in the serum of cancer patients did not significantly change in solid tumor tissue, we sought to identify other types of 5′ tRNA halves that may change in solid tumor tissue independently of the changes in the serum levels. We used a statistical test appropriate for paired designs in EdgeR (39) to analyze the expression of 5′ tRNA halves in tumor and matched normal tissues. Here, even though the biological coefficient of variation (BCV) showed separation between normal and cancer samples, the homogeneity between replicates was less marked than in the serum samples (Figures 1A,B). The differential expression analysis detected six 5′ tRNA halves that are highly expressed (CPM > 12,630) and significantly increased (Fold change ≥ 1.8 and FDR < 4.0%) in tumor compared to matched normal tissues (Table 2). These differentially expressed 5′ tRNA halves are derived from two isodecoders of tRNA-Gly and tRNA-Val. Upon examination of the serum levels of these six 5′ tRNA halves that increased expression in OSCC tissues, we found that all five 5′ tRNA halves derived from tRNA-Gly remain unchanged in the serum of cancer patients (Table 1). As reported above, only the 5′ tRNA half derived from tRNA-Val-CAC-2-1 gene significantly changed expression in both serum and solid tumor (Tables 1, 2). However, the expression changes of 5′ tRNA-Val-CAC-2-1 half are in opposite directions: it is upregulated in solid tumor and downregulated in serum (Figures 1E,F).

TABLE 2

Table 2. Six 5′ tRNA halves significantly (FDR < 5%) increased in tumor relatively to healthy adjacent tissue in patients with OSCC.

Levels of 5′ YRNA Fragments Decreased in Serum of OSCC Patients But Did Not Significantly Change in Solid Tumor Relatively to Adjacent Healthy Tissue

Similar to the analysis for 5′ tRNA halves, we compared the levels of 5′ YRNA fragments in serum from subjects with and without cancer and found marked separation between normal and cancer samples (Figure 1C, BCV distance 1) and adequate homogeneity of the replicates (Figure 1C, BCV distance 2). The differential expression analysis revealed that OSCC decreased the serum levels of 4 types of 5′ YRNA fragments with average read CPM > 150,000, a fold decrease > 2.5, and FDR < 0.01% (Table 3). In contrast, we found no significant separation between normal and cancer samples (Figure 1D) and no significant differential expression for these specific four types of 5′ YRNA fragments between tumor and matched normal tissues.

TABLE 3

Table 3. 5′ YRNA fragments significantly changed in serum of patients with OSCC.

5′ tRNA-Val-CAC-2-1 Half Is Part of a Key miRNA-Enriched Regulatory Coexpression Module That Highly Significantly Correlate With the Cancer Trait

To gain additional insight into potential functions of the 5′ tRNA halves and 5′ YRNA fragments, we implemented a novel approach based on weighted gene coexpression network analysis (WGCNA) of the pooled small RNA-Seq data including 5′ tRNA halves (this study), 5′ YRNA fragments (this study) and miRNAs from the same samples as previously published (18). WGCNA was originally developed to build gene networks based on a gene coexpression measure that describes the correlation among genes with similar patterns of expression (40, 41). Those sets of genes with similar (correlated) expression patterns are denominated modules and represent non-linear pathways with specific biological function. Since its inception in 2005, this methodology has been applied mainly in the gene expression context, however, the applicability of the method has also been extended to the study of other genetic features such as protein, miRNA, and lncRNA networks (42–46). In this study, we further extend the applicability of WGCNA by applying it to sncRNA sequence data that includes 5′ tRNA halves, 5′ YRNA fragments, and miRNAs. We reasoned that miRNAs and 5′ tRNA halves and/or 5′ YRNAs fragments that are coregulated and involved in common biological pathways will belong to common coexpression modules detected by the WGCNA method. This would allow us, using a concept reminiscent of “guilt by association,” to assign/predict specific biological functions for novel cancer-relevant 5′ tRNA halves and 5′ YRNA fragments based on the larger body of knowledge accumulated for miRNAs.

Our analysis identified a total of 14 coexpressed modules (Figures 2A–D). Two modules, namely Green (with 43 features upregulated in OSCC tissue) and Red (with 27 features downregulated in OSCC tissue), demonstrated a representative coexpressed pattern that significantly correlated with the cancer trait (positively for the Green module with r = 0.83, P = 0.003 and negatively for the Red, with r = −0.78, P = 0.01, Figures 2B, 3A,B). No association was detected for either the Green or the Red module with the age at diagnosis (Figure 2B). In addition, module membership (MM) of the individual genetic features in the Green module demonstrated a highly significant correlation with their respective gene significances (GS) for the cancer trait (r = 0.75, P = 7.1 × 10⁻⁹, Figure 3C). MM quantifies how coexpressed, in absolute terms, a given genetic feature (i.e., miRNA, 5′ tRNA halves, or 5′ YRNA fragments) is with respect to the other features in a particular module (47). The Green module mainly consisted of miRNAs but importantly included six 5′ tRNA halves, one of which (5′ tRNA-Val-CAC-2-1 half, shown enclosed in an oval in Figure 3C) displayed both high membership in the module (MM = 0.83) and high GS (GS = 0.67) to the cancer trait. Although to a lower extent, membership of the genetic features in the Red module also significantly correlated with the GS to the cancer trait (r = 0.45, P = 0.019, Figure 3D). The Red module was mainly composed of miRNAs (unlabeled dots in the plots), with a single 5′ tRNA half (labeled red dot) belonging to the module but with low module membership and low GS to the cancer trait (Figure 3D). No 5′ YRNA fragments were detected in any of the two modules with statistically significant correlations to the cancer trait. However, multiple other 5′ tRNA halves and 5′ YRNA fragments belonged to modules (i.e., Black, Brown, GreenYellow, and Turquoise) that displayed marginally significant trends (0.05 < P < 0.1) with the cancer trait (Figure 2B, Supplementary Figures 1, 2).

FIGURE 2

Figure 2. Identifying cancer-relevant modules by weighted gene coexpression network analysis (WGCNA). The network was created from the weighted correlation matrix generated by the WGCNA package in the R environment. First, an adjacency matrix is calculated, then the topological overlap (TO) to hierarchically cluster genes into coexpression modules (see section Materials and Methods). Final module assignments were made based on module membership. (A) Cluster dendrogram groups genetic features into distinct modules. The y-axis represents a dissimilarity distance (1–TO). Dynamic tree cutting was used to determine modules, by dividing the dendrogram at significant branch points (identifying cancer-relevant modules). (B) Correlations between the module eigengenes (representative module expression pattern) and the cancer tratit. Green and Red modules display strong, highly significant correlations to the cancer trait (P-values between parenthesis, scale bar represents the range of correlation coefficients). (C) Barplot of the average gene significance (GS) for each detected module, equivalent to the average correlation between the module genetic features and the cancer trait. (D) Barplot of the average –log P-value of GS. Two modules: Green and Red have a mean GS P < 0.05 (–log₁₀(P) > 1.3). Color-coding of the plots' bars represent module names, unless specified as for the scale bar in (B).

FIGURE 3

Figure 3. Characteristic coexpression patterns (module eigengenes), gene significance (GS) and module membership for genetic features in OSCC-relevant modules. (A,B) Plots of the module eigengenes (representative coexpression pattern) for the cancer-relevant Green and Red modules (C,D). Scatter plots of correlations between gene significance (GS) and module membership for cancer-relevant Green and Red modules. Color-coding is equivalent to module names. Non-miRNA genetic features are labeled with black fonts in the scatter plots. Oval: most relevant 5′ tRNA half.

5′ tRNA-Val-CAC-2-1 Half Is Part of a Key miRNA-Enriched Regulatory Coexpression Module That Targets Genes Involved in the Negative Regulation of the G1/S Transition of the Mitotic Cell Cycle

As commented above, we reasoned that miRNAs and 5′ tRNA halves and/or 5′ YRNA fragments that are coregulated and involved in common biological pathways will belong to common coexpression modules detected by the WGCNA method. Significant coexpression among miRNAs and 5′ tRNA halves was evident in several modules detected by our analysis (modules containing multiple types of sncRNA, Figure 3, Supplementary Figures 1, 2). Because of its particular relevance to the cancer trait, we focused on the Green module to further explore the biological functions of the module and the “guilt by association” functions of the most relevant 5′ tRNA half belonging to this module (namely, the 5′ tRNA-Val-CAC-2-1 half). For this, we conducted “miRNA-overtargeting” analysis for the top 10 miRNAs (Supplementary Table 1) belonging to the Green module, which demonstrated high module membership and high gene significance to the cancer trait. By definition, significantly overtargeted genes in a particular network are genes that experimentally demonstrate significantly higher number of interactions with select miRNAs than expected by chance (48). By identifying enriched biological processes and functions for the network of significantly overtargeted genes, we could then assign specific functions to the related module(s) and to novel genetic features (e.g., 5′ tRNA halves and/or 5′ YRNA fragments) with high module membership. This analytical approach is supported by our previous experiences and those from others indicating that global miRNA-driven regulatory events commonly occur in a coordinated/cooperative fashion (17, 20, 42, 43, 48–52).

In this study, our analysis detected 135 validated miRNA-targeted genes that are significantly overtargeted by the set of top 10 miRNAs from the Green module (Figure 4A, Supplementary Table 2). Simulation analysis with 100,000 iterations demonstrated that finding such a number (or larger) of transcripts being significantly overtargeted by a set of 10 random miRNAs has a significantly low probability equal to 0.01137. This underscores the relevance of the specific gene network targeted by the cancer-associated Green module miRNAs and tRNA halves. KEGG pathway enrichment analysis on the list of targeted genes identified multiple cancer-related pathways, the cell cycle, and viral infection-related pathways among the most significantly enriched (P < 0.01, FDR < 0.05, Figure 4B, Supplementary Table 3). Interestingly, the enrichment in six viral-related pathways suggests a role for Green module miRNAs and tRNA-halves in the transduction of viral signals that may contribute to the development of OSCC. Supporting our reasoning, accumulating evidence confirm the virus-miRNA-cancer axis (53–59) and recent reports highlight the modulation of tRNA-halves in response to viral infections, including in the context of cancer (60–62). Similar to the KEGG pathway analysis, further analysis of gene ontology (GO) annotations for biological processes enriched among the significantly overtargeted genes highlighted the involvement of the Green coexpression module in the regulation of the cell cycle, specifically in the negative regulation of the G1/S transition of the mitotic cell cycle, and in negative regulation of cell differentiation (Figure 5A, Supplementary Table 4). Key cell cycle and differentiation-related transcripts regulated by members of the Green module are: WEE1, E2F1, RB1, and FBXO31 (Figure 4A, Supplementary Table 4). Remarkably, we were able to experimentally validate the results from our bioinformatic analysis. Using qPCR on tissue RNA, we demonstrated the significant downregulation of FXBO31 and WEE1 (Figures 5B,C) in the tumor tissue compared to healthy adjacent mucosal tissue. Although non-significantly, RB1 and E2F1 also displayed downregulation trends (Figures 5D,E). An additional interesting observation from our network analysis is that the AGO1 transcript is significantly overtargeted by the larger number of cancer-correlated miRNAs (AGO1 is targeted by four out of the top 10 Green module miRNAs, Figure 4A, Supplementary Table 2).

FIGURE 4

Figure 4. Discovering molecular functions for the cancer-associated Green coexpression module. (A) Network of genes that are significantly overtargeted by the top 10 miRNAs in the Green module. Blue circles represent overtargeted genes, which are potentially downregulated in tumor tissue (validation of gene expression subset in next figure). Purple rectangles represent the top 10 miRNAs with highest membership in the Green module and high gene significance to the cancer trait (these miRNAs are upregulated in the OSSC tissue). (B) KEGG pathway enrichment among the Green module-overtargeted genes shown in (A). Companion table for the enrichment analysis is presented in Supplementary Table 3.

FIGURE 5

Figure 5. Confirming cell cycle-related functions for the cancer-associated Green coexpression module. (A) Enrichment for gene ontology (GO, biological processes) annotations among the Green module-overtargeted genes shown in Figure 4A. Companion table for the enrichment analysis is presented in Supplementary Table 4. (B–E) Validation of gene expression changes in tumor tissue (compared to adjacent normal/healthy tissue) by quantitative real time PCR (qPCR). Relative quantification was implemented using the 2^−ΔΔCt method for calculation of fold changes in gene expression. Expression of housekeeping gene β2-microglobulin was used as endogenous reference for normalization. Statistical analysis using the Student's t-test was implemented on the log transformed fold change data.

Discussion

The association between changes in the expression of 5′ tRNA halves and 5′ YRNA fragments and cancer are currently largely descriptive. It is still unclear how changes in the expression of these unusual sncRNAs contribute to the process of carcinogenesis. Nonetheless, there are instances where the abundance and function of particular tRNA halves have been linked to cancer (17, 31, 33, 34, 37, 38). In this study, we are able to observe more changes in 5′ tRNA halves in serum than in tissue samples, indicating a wide array of circulating 5′ tRNA halves available for biomarker investigations in cancer patients.

Interestingly, we found that 5′ tRNA-Val-CAC-2-1 half significantly changed expression in both serum and tissue samples, although it was down-regulated in serum and up-regulated in tissue samples. Similar to 5′ tRNA halves, the 5′ YRNA fragments were more affected in serum than solid tissue samples, further confirming their suitability for blood-based diagnosis. We also found that the 5′ tRNA-Val-CAC-2-1 half belonged to a key miRNA-enriched regulatory coexpression module (denominated Green module) that was upregulated in the cancer tissue and targets genes involved in the negative regulation of the G1/S transition of the mitotic cell cycle and cell differentiation, namely FBXO31, WEE1, RB1, and E2F1. We also noticed in the network of Green module-overtargeted genes that AGO1 was the transcript overtargeted by the larger number of cancer-correlated miRNAs, which suggest a central, hub-like role for this gene in mediating the network effects of the Green module. Interestingly, AGO1 has been reported to directly bind tRNA precursors and to be potentially involved in the biogenesis of a subset of 5′ tRNA-derived fragments. Such 5′ tRNA fragments were predicted to be involved in gene silencing (63). More recently, tRNA fragments that act as tumor suppressors were reported to downregulate RNA-binding protein YBX1, which consequently led to destabilization of many pro-oncogenic transcripts (64). Our results now suggest pro-oncogenic 5′ tRNA half-mediated regulation in the opposite direction. Altogether, the evidence suggests that the miRNAs and 5′ tRNA halves (in particular 5′ tRNA-Val-CAC-2-1 half) belonging to the oncogenic Green module may function, at least in part, by downregulating the expression of AGO1 and consequently downregulating the expression of a subset of tumor suppressors 5′ tRNA halves. Therefore, our results and the existent literature support the notion that cancer regulatory mechanisms involve the mutually exclusive action of oncogenic and tumor suppressor subsets of sncRNAs comprised of coregulated miRNAs and 5′ tRNA halves, among other genetic players. Notably, the transcriptional regulatory activity of the oncogenic Green module, including 5′ tRNA-Val-CAC-2-1 half, appears to result in the activation of the cell cycle by downregulating the expression of genes involved in the negative regulation of the G1/S transition (downregulation of these genes in our study samples was experimentally confirmed by qPCR). Supporting our findings, Green module-overtargeted FBXO3 has been indeed reported to be downregulated in multiple cancers and to suppress cell proliferation partly through the degradation of cyclin D1 (65). On the other hand, the nuclear protein Wee1 (the product of the Green module-overtargeted WEE1 gene) was recently recognized as a tumor suppressor that regulates the progression of pancreatic cancer cells via delaying G2, which allows for longer time to repair post-replication errors in genomically unstable cancers (66). Similarly, the RB1 gene is mutated in a variety of human cancers and represents a classical prototype of tumor suppressor, inhibiting cell proliferation primarily through targeting of the E2F transcription factor and, additionally, through transcription-independent mechanisms (67). Although the E2F is hyperactive in many human cancers, its effects are reported to be context dependent, dichotomous and contradictory in almost all cases (68). Supporting our results, Zhong and colleagues recently reported that loss of E2F1 improves survival and accelerates oral tumor growth in mice (69). The tumor-suppressive effect of E2F1 has been ascribed to occur at the cancer promotion stage and may involve the induction of autophagy and apoptosis (68, 70).

Additional evidence support the role of tRNA-derived fragments, including 5′ tRNA halves, in the regulation of mRNA stability and translation, thereby affecting cancer cell proliferation and metastasis (64, 71–73). For example, 5′ tRNA halves induced by hypoxic stress were shown to suppress breast cancer metastasis, which was further shown to only proceed when specific 5′ tRNA halves were antagonized (73). Likewise, a specific 5′-half of tRNA-Asp-GUC was shown to be required for the proliferation of prostate carcinoma cells and knockdown of the 5′ tRNA-Asp-GUC half, but not its 3′ counterpart, interfered with cell proliferation (74). Moreover, dysregulation of other 5′ tRNA halves has been correlated with cancer staging and grading levels (73, 75). Expression of 5′ tRNA-Leu-CAG half was altered in tumor tissue and serum from patients with non-small cell lung cancer with a positive correlation between tumor development stage and the 5′ tRNA half expression levels (73). Expression of 5′ tRNA-4-Val-AAC was altered in renal cell carcinoma tissue samples with inversed correlation between expression levels and tumor staging and grading (75). It was also reported that specific 5′ tRNA halves were produced more abundantly in androgen receptor (AR)-positive prostate cancer and in estrogen receptor (ER)-positive breast cancer than in the respective negative forms (74). That investigation demonstrated that biogenesis of specific 5′ tRNA halves may be sex hormone-dependent (74). Additionally, given the well-established role of sex hormones in the control of proliferation and progression of prostate and breast cancers (76, 77), the observed link between the biogenesis of 5′ tRNA halves and alterations in the abundance of sex hormones may imply that 5′ tRNA halves take part in tumorigenesis of hormone-sensitive cancers. Similarly, we previously reported a significant association between alterations in the circulating levels of specific 5′ tRNA halves and the ER status and other clinicopathologic characteristics of breast cancer including PR, HER2, inflammation, and subsequent relapse (31).

In addition to our findings in this study that serum levels of 5′ YRNA fragments change in patients with OSCC, we have previously reported significant changes in serum and tissues levels of specific YRNA-derived fragments in patients with breast cancer (31) and in serum of patients with HNSCC (17). This association between changes in the abundance of YRNA-derived fragments and cancer suggests a possible function related to carcinogenesis. Supporting evidence includes not only the presence of YRNAs or their derivatives in several types of cancer (78, 79), but also their involvement in processes closely related to cancer such as cell proliferation (80, 81), apoptosis (82–84), senescence (85–87), and stress responses (88). Moreover, it was suggested that YRNA fragments may play a role in the mechanisms by which cancer cells establish a favorable microenvironment for proliferation and invasion (83). Furthermore, the expression of full length YRNAs was found to be dysregulated in prostate (89) and bladder (90) cancers, and associated with poor prognosis in both cancer types.

Our findings with differential expression of 5′ tRNA halves and 5′ YRNA fragments in OSCC patients is in alignment with previous studies with miRNAs (18). Volinia and colleagues (7) had identified 21 miRNAs commonly regulated in solid lung, breast, stomach, prostate, colon, and pancreatic tumors. In previous work of our group, we observed that four of these 21 miRNAs were also differentially regulated between healthy and tumor tissue in OSCC patients (18), suggesting that these miRNAs are important and consistently associated with different types of cancer. In the same way, a previous case-control study from our group identified several serum miRNAs deregulated in patients with HNSCC (17). Several of those were also deregulated in the tissue of OSCC (18), further suggesting these miRNAs as biomarkers for non-invasive diagnosis. Based on these, we suggest that besides miRNAs, the 5′ tRNA halves and 5′ YRNA fragments found deregulated in the current study should be further studied and can be targets not only in HNSCC, but other types of cancers as well.

In conclusion, we found that serum levels of 5′ tRNA halves and 5′ YRNA fragments changed in association with OSCC and that levels of 5′ tRNA halves changed in cancer tissues relatively to adjacent healthy tissues. These and comparable previous findings imply that 5′ tRNA halves and 5′ YRNA fragments may play a role in tumorigenesis. By additionally conducting weighted gene coexpression network analysis of 5′ tRNA halves, 5′ YRNA fragments, and miRNAs measured in the same samples and assessing the functional enrichment in a network of relevant overtargeted genes, we uncovered a potential novel mechanism by which coregulated modules of miRNAs and tRNA halves regulate OSCC. Our study has important limitations including the small sample size and the fact that all tumor samples covered a single TNM classification (i.e., pT2 N1 M0), which limits the generalizability of our findings and warrants validation in independent cohorts (ongoing efforts from our group and others). However, we were able to identify significant differences that underscore the important roles of novel sncRNAs such as 5′ tRNA-Val-CAC-2-1 half in OSSC and were able to experimentally validate our main claims from the functional network analysis. Of note, this is the first time, to our knowledge, that the specific molecular function of a 5′-tRNA half is specifically pinpointed in OSCC. This particular sncRNA is suggested to have potential for use as both, biomarker and therapeutic target in OSCC.

Materials and Methods

Patients

Tumor samples from five patients (2 males and 3 females, 64.4 ± 5.6 years old) diagnosed with HNSCC were used. All patients were qualified for the primary surgical resection by the multidisciplinary team from the institution. Recurrences and patients initially treated with other therapeutic approaches were excluded, as well as HPV positive tumors. All tumors were confirmed by a pathologist and have been collected intraoperatively as it was previously published (91, 92). All samples were classified based on the guidelines from the International Union Against Cancer using standard TNM classification. Three of the five samples were located at the floor of the mouth and the remaining two in the base of the tongue. All samples were classified as pT2 N1 M0. This study was carried out in accordance with the recommendations and approval by Institutional Review Board of the University of Medical Sciences in Poznan. All subjects gave written informed consent in accordance with the Declaration of Helsinki.

Sample and Tissue Collection

Blood samples were collected prior surgical intervention and centrifuged for serum separation. Additionally, two separate tissue samples were collected during surgical resection, one from the tumor and one from a healthy mucosa adjacent to the tumor in the same patient as described previously (18). Samples were immediately frozen in liquid nitrogen and stored ay −80°C.

RNA Isolation and Small RNA Library Construction

As previously described in reference (18), tissues samples were removed from the −80°C freezer and homogenized with Qiazol (Qiagen, Valencia, CA, USA) using 0.5 mm zirconium oxide beads in the Bullet Blender 24 (Next Advance, Averill Park, NY, USA). Total RNA was extracted using a commercial column purification system (miRNeasy Mini Kit, Qiagen) and on-column DNase treatment (RNase-free DNase Set, Qiagen) following manufacturer's instructions. Serum samples were extracted using the miRNEasy Serum/Plasma kit (Qiagen) following manufactures instructions but adjusting for an initial volume of 300 μL of serum. MicroRNA libraries were prepared using the TruSeq Small RNA Sample Preparation Kit (Illumina Inc., San Diego, CA, USA) following the manufacturer's instructions and adjusted by Matkovich et al. (93). Briefly, small RNAs from 1 μg of total RNA or from 300 μL of serum (not quantified) were ligated with 3′ and 5′ adapters, followed by reverse transcription to produce single stranded cDNAs. Samples were then amplified by PCR in 14 cycles (94°C for 30 s, 14 cycles of 94°C for 15 s, 62°C for 30 s, and 70°C for 15 s, and a final extension of 70°C for 5 min) using indexes to allow all individual libraries to be processed in a single flowcell lane during sequencing. The amplified libraries were size-selected and purified in a 6% agarose gel. The quantity and quality of miRNA libraries was determined using BioAnalyzer and RNA Nano Lab Chip Kit (Agilent Technologies, Santa Clara, CA, USA), and the samples were combined in a single microtube and submitted to sequencing on a HiSeq. 2500 instrument (Illumina Inc.) using the Illumina HiSeq v4 kit in a single read 50 bp (1 × 50) run.

Quantitative Real Time PCR (qPCR)

Confirmation of tissue gene expression changes in genes targeted by cancer-associated miRNAs and tRNA halves was conducted by qPCR as described by Allen and collaborators (94) with some modifications. In summary, RNA extracted from tumor and normal adjacent tissue as described in previous section was converted to cDNA using the iScript cDNA synthesis kit (Bio-Rad) following the manufacturer instructions. qPCR reactions were set up in a MicroAmp^® Fast Optical 96-well reaction plate (Applied Biosystems) with 2 μL of diluted cDNA, 0.2 μL each of forward and reverse primer, 12.6 μL of nuclease free water, and 5 μL of Fast SYBR Green Master Mix (Applied Biosystems) per well. Amplification was performed in a 7900 HT Fast system (Applied Biosystems) at 95°C for 20 s, followed by 45 cycles of 1 s denaturation at 95°C and 20 s annealing/extension at 62°C. Primer sequences used were: FBXO31 (forward: 5′-GCCGTGAGGAGTATGGTGTTT-3′ and reverse: 5′-GTACATCCACCCGATGATGAAC-3′), WEE1 (forward: 5′-CTTGGGGACTTCTGCATGA-3′ and reverse: 5′-GCTTGGGGACTATCACCACT-3′), RB1 (forward: 5′-CTCTCGTCAGGCTTGAGTTTG-3′ and reverse: 5′-GACATCTCATCTAGGTCAACTGC-3′), and E2F1 (forward: 5′-CTACGTGACGTTCAGGACC-3′ and reverse: 5′-CCGGAGTTCCCGATCTAC-3′). β2-microglobulin (forward primer: 5′-GAGTATGCCTGCCGTGTGAA-3′ and reverse primer: 5′-CGGCATCTTCAAACCTCCAT-3′) was used as the housekeeping gene to normalize the qPCR data. Relative quantification was implemented using the 2^−ΔΔCt method for calculation of fold changes in gene expression. Statistical analysis using the Student's t-test was implemented on the log transformed fold change data.

Bioinformatics and Statistical Analyses of the Circulating tRNA- and YRNA-Derived Small RNAs in Patients With and Without Cancer

The sequencing reads obtained from control and cancer serum samples were trimmed by removing adaptor sequences with FASTX-Toolkit (hannonlab.cshl.edu) and mapped to the GRCh38/hg38 human genome with Bowtie v1.1.2 [35]. The obtained aligned reads were analyzed to annotate and count the reads that map to tRNA genes from Genomic tRNA Database (95, 96), YRNA genes from the UCSC GENCODE v24 track (97), and non-coding RNA genes from Ensembl GRCh38 release 86.

To assess differences in the abundance of circulating tRNA- and YRNA-derived small RNAs between cancer and control groups, we analyzed the read count using the negative binomial model in the Bioconductor package edgeR (39). After normalization of the read counts with the trimmed mean of M-values (TMM) method, the differential expression was assessed by the exact test which is applicable to experiments with a single factor. P-values were adjusted for multiple testing using the Benjamini and Hochberg method to control the false discovery rate (FDR). Differences in expression were considered significant below an FDR of 5%.

Bioinformatics and Statistical Analyses of tRNA- and YRNA-Derived Small RNAs in Tumor and Matched Normal Tissues

Sequencing reads from tissue of tumor-normal pairs were aligned to the GRCh38/hg38 human genome using the same parameters used for the analysis of the small RNA-Seq from the serum samples described above. To detect statistically significant differential expression of tRNA- and YRNA-derived small RNAs, read counts from tissue of tumor-normal pairs were analyzed with the Bioconductor package edgeR using a statistical test appropriate for paired designs (39). The test measures differential gene expression between tumor and normal tissue. To adjust for differences between the patients, the test uses an additive linear model with “Patient” as the blocking factor. P-values were adjusted for multiple testing using the Benjamini and Hochberg method to control the FDR. Differences in expression were considered significant below an FDR of 5%. Multi-dimensional scaling (MDS) analysis of the expression levels of 5′ tRNA-halves in serum and tissues was performed using the edgeR function plotMDS (39).

Weighted Gene Coexpression Network Analysis (WGCNA)

Traditional genetic approaches focusing on the action of single genes have contributed to the identification of important candidate genes for many diseases, particularly monogenetic diseases. However, these approaches have limited effectiveness in identifying genes and other genetic features (e.g., sncRNAs) that contribute to complex diseases such as cancer. On the other hand, the analysis of entire gene and non-coding RNA networks using modern genomics and bioinformatic tools have proven invaluable for the identification of biological pathways underlying disease (98, 99). WGCNA is a network analysis tool based on a coexpression measure that describes the correlation among genetic features with similar patterns of expression. Those sets of features with similar (correlated) expression patterns are denominated modules and represent non-linear pathways of related biological functions. The general framework of WGCNA has been described in detail elsewhere (40, 41, 47).

In this study, WGCNA was implemented as previously reported by Nunez et al. (43) and Gorini et al. (42) with some modifications as described below. Normalized small RNA-Seq data for tissue (cancer vs. normal) 5′ tRNA halves, 5′ YRNA fragments (this work) and miRNA (18) were pooled into a single data table for WGCNA analysis (47) in the R biostatistical computing environment. WGCNA was performed on all tissue samples and included 540 genetic features (i.e., 5′ tRNA halves, 5′ YRNA fragments, and miRNAs) with median absolute deviation (MAD, a robust measure of variability) <0 (out of a total of 1,429 features). A soft power of β = 6 was chosen to produce a network with approximate scale-free topology (soft R² = 0.8). Parameter deepSplit was set to 3. Because the number of genetic features in our analysis is lower than the number of genes usually detected in a microarray study (the data type for which WGCNA was originally developed), we reduced the minimum module size parameter to 10. All other parameters were used as defined by default. Gene modules corresponding to the branches cut off of the gene tree were labeled in unique colors. Unassigned genes were assigned to the Gray module. miRNA-overtargeted gene (validated interaction) networks were constructed for the Green module as described below. The trait-based gene significance (GS) measure is defined in WGCNA as the absolute value of the correlation between the expression profile of a specific node (genetic feature such as miRNA, tRNA half, or YRNA fragment in our case) and the sample trait (cancer vs. normal tissue). The higher the value of GS for a specific genetic feature, the more biologically significant the specific genetic feature is in the context of the specific trait. On the other hand, module membership, also known as eigengene-based connectivity, correlates the expression profile of the specific genetic feature with the module eigengene (equivalent to a representative expression pattern) of a given module. The higher the value of MM for a specific genetic feature, the higher the contribution of the specific genetic feature to the overall pattern and possibly function of the specific module (pathway). By focusing on modules with high correlation to sample traits of interest (cancer in our case), and subsequently identifying within those modules the genetic features with higher GS and MM, one can identify the more important pathways and genetic features contributing to the disease trait of interest.

Functional Network Analysis: miRNA-Overtargeted Analysis, Network Construction, and Enrichment Analysis of Gene Ontology Annotations

Functional enrichment analysis was performed as previously described (48). Significantly overtargeted genes in a particular network are defined as those genes that experimentally demonstrate significantly higher number of interactions with select miRNAs than expected by chance (as determined by hypergeometric tests of network proportions compared to the respective proportions in the miRNA–target interaction universe/background) (48). By identifying enriched biological processes and functions for the significantly overtargeted genes, we could assign specific functions to the related module(s) and to novel genetic features with high module membership as calculated using WGCNA. This miRNA-overtargeting analysis approach is supported by our previous experiences and those from others indicating that global miRNA-driven regulatory events commonly occur in a coordinated/cooperative fashion (17, 20, 42, 43, 48–52). In short, significance of the overtargeting effect (significant higher number of miRNA-targeting events than expected by chance on a given gene) was assessed by comparing module-specific network proportions and relevant number of events against a collection of 100,000 simulated equivalent random networks. The list of validated targets supported by strong experimental evidence (i.e., reporter assay or Western blot) used for this analysis was downloaded from miRTarBase (file: miRTarBase_SE_WR.xls) using the SpidermiR package (100). Interaction networks were constructed using Cytoscape 3.5.1 (101). Enrichment of gene ontology annotations among sets of overtargeted genes was assessed using the GOCluster_Report function of the systemPipeR package (102) in the R environment.

Data Availability Statement

Datasets in this study are available at: https://www.ncbi.nlm.nih.gov/bioproject with the accession BioProject: PRJNA561150.

Ethics Statement

This study was carried out in accordance with the recommendations and approval by Institutional Review Board of the University of Medical Sciences in Poznan. All subjects gave written informed consent in accordance with the Declaration of Helsinki.

Author Contributions

Individual contributions to the present work was as follows: JD, YN, and MM: conceptualization and methodology. JD, KB, TM, YN, AS, and BV: software. JD, YN, AS, BV, TS, KB, TM, HA, WS, PG, WG, and MM: validation, data curation, manuscript writing, reviewing, and editing. JD, YN, AS, and MM: formal analysis. JD, AS, BV, TS, KB, TM, HA, WS, PG, and WG: experimental investigation. JD, YN, PG, WG, and MM: resources. JD and MM: supervision, project administration, and funding acquisition.

Funding

This work was supported by the National Institute on Aging (NIA) grants R15 AG059190 and R03 AG059846.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2019.00959/full#supplementary-material

References

1. Rezende TM, De Souza Freire M, Franco OL. Head and neck cancer: proteomic advances and biomarker achievements. Cancer. (2010) 116:4914–25. doi: 10.1002/cncr.25245

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Gugic J, Strojan P. Squamous cell carcinoma of the head and neck in the elderly. Rep Pract Oncol Radiother. (2012) 18:16–25. doi: 10.1016/j.rpor.2012.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Rossini AR, Hashimoto CL, Iriya K, Zerbini C, Baba ER, Moraes-Filho JP. Dietary habits, ethanol and tobacco consumption as predictive factors in the development of esophageal carcinoma in patients with head and neck neoplasms. Dis Esophagus. (2008) 21:316–21. doi: 10.1111/j.1442-2050.2007.00769.x

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Chaturvedi AK, Engels EA, Pfeiffer RM, Hernandez BY, Xiao W, Kim E, et al. Human papillomavirus and rising oropharyngeal cancer incidence in the United States. J Clin Oncol. (2011) 29:4294–301. doi: 10.1200/JCO.2011.36.4596

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Gillison ML, Chaturvedi AK, Anderson WF, Fakhry C. Epidemiology of human papillomavirus-positive head and neck Squamous Cell Carcinoma. J Clin Oncol. (2015) 33:3235–42. doi: 10.1200/JCO.2015.61.6995

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Golusinski P, Pazdrowski J, Szewczyk M, Misiolek M, Pietruszewska W, Klatka J, et al. Is immunohistochemical evaluation of p16 in oropharyngeal cancer enough to predict the HPV positivity? Rep Pract Oncol Radiother. (2017) 22:237–42. doi: 10.1016/j.rpor.2017.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Volinia S, Calin GA, Liu CG, Ambs S, Cimmino A, Petrocca F, et al. A microRNA expression signature of human solid tumors defines cancer gene targets. Proc Natl Acad Sci USA. (2006) 103:2257–61. doi: 10.1073/pnas.0510565103

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Mizoguchi M, Guan Y, Yoshimoto K, Hata N, Amano T, Nakamizo A, et al. Clinical implications of microRNAs in human glioblastoma. Front Oncol. (2013) 3:19. doi: 10.3389/fonc.2013.00019

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Prokopi M, Kousparou CA, Epenetos AA. The secret role of microRNAs in cancer stem cell development and potential therapy: a notch-pathway approach. Front Oncol. (2014) 4:389. doi: 10.3389/fonc.2014.00389

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Kolenda T, Guglas K, Rys M, Bogaczynska M, Teresiak A, Blizniak R, et al. Biological role of long non-coding RNA in head and neck cancers. Rep Pract Oncol Radiother. (2017) 22:378–88. doi: 10.1016/j.rpor.2017.07.001

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Nunez Lopez YO, Victoria B, Golusinski P, Golusinski W, Masternak MM. Characteristic miRNA expression signature and random forest survival analysis identify potential cancer-driving miRNAs in a broad range of head and neck squamous cell carcinoma subtypes. Rep Pract Oncol Radiother. (2018) 23:6–20. doi: 10.1016/j.rpor.2017.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Calin GA, Croce CM. MicroRNA signatures in human cancers. Nat Rev Cancer. (2006) 6:857–66. doi: 10.1038/nrc1997

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Kolenda T, Przybyla W, Kapalczynska M, Teresiak A, Zajaczkowska M, Blizniak R, et al. Tumor microenvironment - unknown niche with powerful therapeutic potential. Rep Pract Oncol Radiother. (2018) 23:143–53. doi: 10.1016/j.rpor.2018.01.004

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Pardini B, Calin GA. MicroRNAs and long non-coding RNAs and their hormone-like activities in cancer. Cancers. (2019) 11:378. doi: 10.3390/cancers11030378

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Bartel DP. Metazoan MicroRNAs. Cell. (2018) 173:20–51. doi: 10.1016/j.cell.2018.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Shiiba M, Uzawa K, Tanzawa H. MicroRNAs in head and neck squamous cell carcinoma (HNSCC) and oral squamous cell carcinoma (OSCC). Cancers (Basel). (2010) 2:653–69. doi: 10.3390/cancers2020653

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Victoria Martinez B, Dhahbi JM, Nunez Lopez YO, Lamperska K, Golusinski P, Luczewski L, et al. Circulating small non-coding RNA signature in head and neck squamous cell carcinoma. Oncotarget. (2015) 6:19246–63. doi: 10.18632/oncotarget.4266

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Schneider A, Victoria B, Lopez YN, Suchorska W, Barczak W, Sobecka A, et al. Tissue and serum microRNA profile of oral squamous cell carcinoma patients. Sci Rep. (2018) 8:675. doi: 10.1038/s41598-017-18945-z

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Dhahbi JM, Spindler SR, Atamna H, Yamakawa A, Guerrero N, Boffelli D, et al. Deep sequencing identifies circulating mouse miRNAs that are functionally implicated in manifestations of aging and responsive to calorie restriction. Aging. (2013) 5:130–41. doi: 10.18632/aging.100540

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Victoria B, Dhahbi JM, Nunez Lope YO, Spinel L, Atamna H, Spindler SR, et al. Circulating microRNA signature of genotype-by-age interactions in the long-lived Ames dwarf mouse. Aging Cell. (2015) 14:1055–66. doi: 10.1111/acel.12373

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Dhahbi JM, Atamna H, Li R, Yamakawa A, Guerrero N, Lam HT, et al. MicroRNAs circulate in the hemolymph of drosophila and accumulate relative to tissue microRNAs in an age-dependent manner. Genomics Insights. (2016) 9:29–39. doi: 10.4137/GEI.S38147

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Schneider A, Dhahbi JM, Atamna H, Clark JP, Colman RJ, Anderson RM. Caloric restriction impacts plasma microRNAs in rhesus monkeys. Aging Cell. (2017) 16:1200–3. doi: 10.1111/acel.12636

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Calin GA, Croce CM. MicroRNA-cancer connection: the beginning of a new tale. Cancer Res. (2006) 66:7390–4. doi: 10.1158/0008-5472.CAN-06-0800

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Becker A, Thakur BK, Weiss JM, Kim HS, Peinado H, Lyden D. Extracellular vesicles in cancer: cell-to-cell mediators of metastasis. Cancer Cell. (2016) 30:836–48. doi: 10.1016/j.ccell.2016.10.009

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Lopatina T, Gai C, Deregibus MC, Kholia S, Camussi G. Cross talk between cancer and mesenchymal stem cells through extracellular vesicles carrying nucleic acids. Front Oncol. (2016) 6:125. doi: 10.3389/fonc.2016.00125

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Blackwell RH, Foreman KE, Gupta GN. The role of cancer-derived exosomes in tumorigenicity & epithelial-to-mesenchymal transition. Cancers. (2017) 9:105. doi: 10.3390/cancers9080105

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Vlaeminck-Guillem V. Extracellular vesicles in prostate cancer carcinogenesis, diagnosis, and management. Front Oncol. (2018) 8:222. doi: 10.3389/fonc.2018.00222

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Dhahbi JM, Spindler SR, Atamna H, Boffelli D, Mote P, Martin DI. 5′-YRNA fragments derived by processing of transcripts from specific YRNA genes and pseudogenes are abundant in human serum and plasma. Physiol Genomics. (2013) 45:990–8. doi: 10.1152/physiolgenomics.00129.2013

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Dhahbi JM, Spindler SR, Atamna H, Yamakawa A, Boffelli D, Mote P, et al. 5′ tRNA halves are present as abundant complexes in serum, concentrated in blood cells, and modulated by aging and calorie restriction. BMC Genomics. (2013) 14:298. doi: 10.1186/1471-2164-14-298

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Dhahbi JM. Circulating small noncoding RNAs as biomarkers of aging. Ageing Res Rev. (2014) 17:86–98. doi: 10.1016/j.arr.2014.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Dhahbi JM, Spindler SR, Atamna H, Boffelli D, Martin DI. Deep sequencing of serum small RNAs identifies patterns of 5′ tRNA half and YRNA fragment expression associated with breast cancer. Biomark Cancer. (2014) 6:37–47. doi: 10.4137/BIC.S20764

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Dhahbi JM. 5′ tRNA halves: the next generation of immune signaling molecules. Front Immunol. (2015) 6:74. doi: 10.3389/fimmu.2015.00074

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Dhahbi JM, Atamna H, Selth LA. Data mining of small RNA-seq suggests an association between prostate cancer and altered abundance of 5′ transfer RNA halves in seminal fluid and prostatic tissues. Biomark Cancer. (2018) 10:1179299X18759545. doi: 10.1177/1179299X18759545

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Selitsky SR, Baran-Gale J, Honda M, Yamane D, Masaki T, Fannin EE, et al. Small tRNA-derived RNAs are increased and more abundant than microRNAs in chronic hepatitis B and C. Sci Rep. (2015) 5:7675. doi: 10.1038/srep07675

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Kaudewitz D, Skroblin P, Bender LH, Barwari T, Willeit P, Pechlaner R, et al. Association of MicroRNAs and YRNAs with platelet function. Circ Res. (2016) 118:420–32. doi: 10.1161/CIRCRESAHA.114.305663

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Chiou NT, Kageyama R, Ansel KM. Selective export into extracellular vesicles and function of tRNA fragments during T cell activation. Cell Rep. (2018) 25:3356–70 e3354. doi: 10.1016/j.celrep.2018.11.073

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Shen Y, Yu X, Zhu L, Li T, Yan Z, Guo J. Transfer RNA-derived fragments and tRNA halves: biogenesis, biological functions and their roles in diseases. J Mol Med. (2018) 96:1167–76. doi: 10.1007/s00109-018-1693-y

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Zhao C, Tolkach Y, Schmidt D, Muders M, Kristiansen G, Muller SC, et al. tRNA-halves are prognostic biomarkers for patients with prostate cancer. Urol Oncol. (2018) 36:503 e501–3 e507. doi: 10.1016/j.urolonc.2018.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Robinson MD, Mccarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Zhang B, Horvath S. A general framework for weighted gene co-expression network analysis. Stat Appl Genet Mol Biol. (2005) 4:Article17. doi: 10.2202/1544-6115.1128

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Horvath S, Dong J. Geometric interpretation of gene coexpression network analysis. PLoS Comput Biol. (2008) 4:e1000117. doi: 10.1371/journal.pcbi.1000117

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Gorini G, Nunez YO, Mayfield RD. Integration of miRNA and protein profiling reveals coordinated neuroadaptations in the alcohol-dependent mouse brain. PLoS ONE. (2013) 8:e82565. doi: 10.1371/annotation/6286be0f-d729-495a-8a72-78995e9ceda7

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Nunez YO, Truitt JM, Gorini G, Ponomareva ON, Blednov YA, Harris RA, et al. Positively correlated miRNA-mRNA regulatory networks in mouse frontal cortex during early stages of alcohol dependence. BMC Genomics. (2013) 14:725. doi: 10.1186/1471-2164-14-725

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Zhao F, Ge YZ, Zhou LH, Xu LW, Xu Z, Ping WW, et al. Identification of hub miRNA biomarkers for bladder cancer by weighted gene coexpression network analysis. Onco Targets Ther. (2017) 10:5551–9. doi: 10.2147/OTT.S146479

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Alaei S, Sadeghi B, Najafi A, Masoudi-Nejad A. LncRNA and mRNA integration network reconstruction reveals novel key regulators in esophageal squamous-cell carcinoma. Genomics. (2018) 111:76–89. doi: 10.1016/j.ygeno.2018.01.003

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Oliveira GB, Regitano LCA, Cesar ASM, Reecy JM, Degaki KY, Poleti MD, et al. Integrative analysis of microRNAs and mRNAs revealed regulation of composition and metabolism in Nelore cattle. BMC Genomics. (2018) 19:126. doi: 10.1186/s12864-018-4514-3

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. (2008) 9:559. doi: 10.1186/1471-2105-9-559

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Nunez Lopez YO, Retnakaran R, Zinman B, Pratley RE, Seyhan AA. Predicting and understanding the response to short-term intensive insulin therapy in people with early type 2 diabetes. Mol Metab. (2019) 20:63–78. doi: 10.1016/j.molmet.2018.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Laurent LC, Chen J, Ulitsky I, Mueller FJ, Lu C, Shamir R, et al. Comprehensive microRNA profiling reveals a unique human embryonic stem cell signature dominated by a single seed sequence. Stem Cells. (2008) 26:1506–16. doi: 10.1634/stemcells.2007-1081

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Lewohl JM, Nunez YO, Dodd PR, Tiwari GR, Harris RA, Mayfield RD. Up-regulation of microRNAs in brain of human alcoholics. Alcohol Clin Exp Res. (2011) 35:1928–37. doi: 10.1111/j.1530-0277.2011.01544.x

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Xiao Y, Xu C, Guan J, Ping Y, Fan H, Li Y, et al. Discovering dysfunction of multiple microRNAs cooperation in disease by a conserved microRNA co-expression network. PLoS ONE. (2012) 7:e32201. doi: 10.1371/journal.pone.0032201

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Shao T, Wang G, Chen H, Xie Y, Jin X, Bai J, et al. Survey of miRNA-miRNA cooperative regulation principles across cancer types. Brief Bioinform. (2018). doi: 10.1093/bib/bby038. [Epub ahead of print].

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Pfeffer S, Voinnet O. Viruses, microRNAs and cancer. Oncogene. (2006) 25:6211–9. doi: 10.1038/sj.onc.1209915

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Boix L, Bruix J. MicroRNA analysis for outcome prediction in hepatitis B virus-related liver cancer. Gastroenterology. (2010) 138:1624–6; discussion: 1626–7. doi: 10.1053/j.gastro.2010.02.021

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Dudda JC, Salaun B, Ji Y, Palmer DC, Monnot GC, Merck E, et al. MicroRNA-155 is required for effector CD8⁺ T cell responses to virus infection and cancer. Immunity. (2013) 38:742–53. doi: 10.1016/j.immuni.2012.12.006

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Jimenez-Wences H, Peralta-Zaragoza O, Fernandez-Tilapa G. Human papilloma virus, DNA methylation and microRNA expression in cervical cancer (review). Oncol Rep. (2014) 31:2467–76. doi: 10.3892/or.2014.3142

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Pandya D, Mariani M, Mchugh M, Andreoli M, Sieber S, He S, et al. Herpes virus microRNA expression and significance in serous ovarian cancer. PLoS ONE. (2014) 9:e114750. doi: 10.1371/journal.pone.0114750

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Wang X, Liu S, Zhou Z, Yan H, Xiao J. A herpes simplex virus type 2-encoded microRNA promotes tumor cell metastasis by targeting suppressor of cytokine signaling 2 in lung cancer. Tumour Biol. (2017) 39:1010428317701633. doi: 10.1177/1010428317701633

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Song F, Wei M, Wang J, Liu Y, Guo M, Li X, et al. Hepatitis B virus-regulated growth of liver cancer cells occurs through the microRNA-340-5p-activating transcription factor 7-heat shock protein A member 1B axis. Cancer Sci. (2019) 110:1633–43. doi: 10.1111/cas.14004

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Ruggero K, Guffanti A, Corradin A, Sharma VK, De Bellis G, Corti G, et al. Small noncoding RNAs in cells transformed by human T-cell leukemia virus type 1: a role for a tRNA fragment as a primer for reverse transcriptase. J Virol. (2014) 88:3612–22. doi: 10.1128/JVI.02823-13

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Deng J, Ptashkin RN, Chen Y, Cheng Z, Liu G, Phan T, et al. Respiratory syncytial virus utilizes a tRNA fragment to suppress antiviral responses through a novel targeting mechanism. Mol Ther. (2015) 23:1622–9. doi: 10.1038/mt.2015.124

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Taxis TM, Bauermann FV, Ridpath JF, Casas E. Analysis of tRNA halves (tsRNAs) in serum from cattle challenged with bovine viral diarrhea virus. Genet Mol Biol. (2019) 42:374–9. doi: 10.1590/1678-4685-gmb-2018-0019

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Burroughs AM, Ando Y, De Hoon MJ, Tomaru Y, Suzuki H, Hayashizaki Y, et al. Deep-sequencing of human Argonaute-associated small RNAs provides insight into miRNA sorting and reveals Argonaute association with RNA fragments of diverse origin. RNA Biol. (2011) 8:158–77. doi: 10.4161/rna.8.1.14300

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Goodarzi H, Liu X, Nguyen HC, Zhang S, Fish L, Tavazoie SF. Endogenous tRNA-Derived Fragments Suppress Breast Cancer Progression via YBX1 Displacement. Cell. (2015) 161:790–802. doi: 10.1016/j.cell.2015.02.053

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Wang Z, Liu P, Inuzuka H, Wei W. Roles of F-box proteins in cancer. Nat Rev Cancer. (2014) 14:233–47. doi: 10.1038/nrc3700

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Chang Q, Chandrashekhar M, Ketela T, Fedyshyn Y, Moffat J, Hedley D. Cytokinetic effects of Wee1 disruption in pancreatic cancer. Cell Cycle. (2016) 15:593–604. doi: 10.1080/15384101.2016.1138188

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Dyson NJ. RB1: a prototype tumor suppressor and an enigma. Genes Dev. (2016) 30:1492–502. doi: 10.1101/gad.282145.116

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Meng P, Ghosh R. Transcription addiction: can we garner the Yin and Yang functions of E2F1 for cancer therapy? Cell Death Dis. (2014) 5:e1360. doi: 10.1038/cddis.2014.326

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Zhong R, Bechill J, Spiotto MT. Loss of E2F1 extends survival and accelerates oral tumor growth in HPV-positive mice. Cancers (Basel). (2015) 7:2372–85. doi: 10.3390/cancers7040895

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Pierce AM, Schneider-Broussard R, Gimenez-Conti IB, Russell JL, Conti CJ, Johnson DG. E2F1 has both oncogenic and tumor-suppressive properties in a transgenic model. Mol Cell Biol. (1999) 19:6408–14. doi: 10.1128/MCB.19.9.6408

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Lee YS, Shibata Y, Malhotra A, Dutta A. A novel class of small RNAs: tRNA-derived RNA fragments (tRFs). Genes Dev. (2009) 23:2639–49. doi: 10.1101/gad.1837609

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Maute RL, Schneider C, Sumazin P, Holmes A, Califano A, Basso K, et al. tRNA-derived microRNA modulates proliferation and the DNA damage response and is down-regulated in B cell lymphoma. Proc Natl Acad Sci USA. (2013) 110:1404–9. doi: 10.1073/pnas.1206761110

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Shao Y, Sun Q, Liu X, Wang P, Wu R, Ma Z. tRF-Leu-CAG promotes cell proliferation and cell cycle in non-small cell lung cancer. Chem Biol Drug Des. (2017) 90:730–8. doi: 10.1111/cbdd.12994

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Honda S, Loher P, Shigematsu M, Palazzo JP, Suzuki R, Imoto I, et al. Sex hormone-dependent tRNA halves enhance cell proliferation in breast and prostate cancers. Proc Natl Acad Sci USA. (2015) 112:E3816–25. doi: 10.1073/pnas.1510077112

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Nientiedt M, Deng M, Schmidt D, Perner S, Muller SC, Ellinger J. Identification of aberrant tRNA-halves expression patterns in clear cell renal cell carcinoma. Sci Rep. (2016) 6:37158. doi: 10.1038/srep37158

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Lonergan PE, Tindall DJ. Androgen receptor signaling in prostate cancer development and progression. J Carcinog. (2011) 10:20. doi: 10.4103/1477-3163.83937

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Haldosen LA, Zhao C, Dahlman-Wright K. Estrogen receptor beta in breast cancer. Mol Cell Endocrinol. (2014) 382:665–72. doi: 10.1016/j.mce.2013.08.005

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Schotte D, Chau JC, Sylvester G, Liu G, Chen C, Van Der Velden VH, et al. Identification of new microRNA genes and aberrant microRNA profiles in childhood acute lymphoblastic leukemia. Leukemia. (2009) 23:313–22. doi: 10.1038/leu.2008.286

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Meiri E, Levy A, Benjamin H, Ben-David M, Cohen L, Dov A, et al. Discovery of microRNAs and other small RNAs in solid tumors. Nucleic Acids Res. (2010) 38:6234–46. doi: 10.1093/nar/gkq376

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Christov CP, Gardiner TJ, Szuts D, Krude T. Functional requirement of noncoding Y RNAs for human chromosomal DNA replication. Mol Cell Biol. (2006) 26:6993–7004. doi: 10.1128/MCB.01060-06

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Krude T, Christov CP, Hyrien O, Marheineke K. Y RNA functions at the initiation step of mammalian chromosomal DNA replication. J Cell Sci. (2009) 122:2836–45. doi: 10.1242/jcs.047563

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Rutjes SA, Van Der Heijden A, Utz PJ, Van Venrooij WJ, Pruijn GJ. Rapid nucleolytic degradation of the small cytoplasmic Y RNAs during apoptosis. J Biol Chem. (1999) 274:24799–807. doi: 10.1074/jbc.274.35.24799

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Chakrabortty SK, Prakash A, Nechooshtan G, Hearn S, Gingeras TR. Extracellular vesicle-mediated transfer of processed and functional RNY5 RNA. RNA. (2015) 21:1966–79. doi: 10.1261/rna.053629.115

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Hizir Z, Bottini S, Grandjean V, Trabucchi M, Repetto E. RNY (YRNA)-derived small RNAs regulate cell death and inflammation in monocytes/macrophages. Cell Death Dis. (2017) 8:e2530. doi: 10.1038/cddis.2016.429

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Clancy RM, Alvarez D, Komissarova E, Barrat FJ, Swartz J, Buyon JP. Ro60-associated single-stranded RNA links inflammation with fetal cardiac fibrosis via ligation of TLRs: a novel pathway to autoimmune-associated heart block. J Immunol. (2010) 184:2148–55. doi: 10.4049/jimmunol.0902248

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Meredith BK, Berry DP, Kearney F, Finlay EK, Fahey AG, Bradley DG, et al. A genome-wide association study for somatic cell score using the Illumina high-density bovine beadchip identifies several novel QTL potentially related to mastitis susceptibility. Front Genet. (2013) 4:229. doi: 10.3389/fgene.2013.00229

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Reed JH, Sim S, Wolin SL, Clancy RM, Buyon JP. Ro60 requires Y3 RNA for cell surface exposure and inflammation associated with cardiac manifestations of neonatal lupus. J Immunol. (2013) 191:110–6. doi: 10.4049/jimmunol.1202849

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Nicolas FE, Hall AE, Csorba T, Turnbull C, Dalmay T. Biogenesis of Y RNA-derived small RNAs is independent of the microRNA pathway. FEBS Lett. (2012) 586:1226–30. doi: 10.1016/j.febslet.2012.03.026

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Tolkach Y, Niehoff EM, Stahl AF, Zhao C, Kristiansen G, Muller SC, et al. YRNA expression in prostate cancer patients: diagnostic and prognostic implications. World J Urol. (2018) 36:1073–8. doi: 10.1007/s00345-018-2250-6

PubMed Abstract | CrossRef Full Text | Google Scholar

90. Tolkach Y, Stahl AF, Niehoff EM, Zhao C, Kristiansen G, Muller SC, et al. YRNA expression predicts survival in bladder cancer patients. BMC Cancer. (2017) 17:749. doi: 10.1186/s12885-017-3746-y

PubMed Abstract | CrossRef Full Text | Google Scholar

91. Zhi X, Lamperska K, Golusinski P, Schork NJ, Luczewski L, Golusinski W, et al. Expression levels of insulin-like growth factors 1 and 2 in head and neck squamous cell carcinoma. Growth Horm IGF Res. (2014) 24:137–41. doi: 10.1016/j.ghir.2014.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

92. Zhi X, Lamperska K, Golusinski P, Schork NJ, Luczewski L, Kolenda T, et al. Gene expression analysis of head and neck squamous cell carcinoma survival and recurrence. Oncotarget. (2015) 6:547–55. doi: 10.18632/oncotarget.2772

PubMed Abstract | CrossRef Full Text | Google Scholar

93. Matkovich SJ, Hu Y, Dorn GW II. Regulation of cardiac microRNAs by cardiac microRNAs. Circul Res. (2013) 113:62–71. doi: 10.1161/CIRCRESAHA.113.300975

PubMed Abstract | CrossRef Full Text | Google Scholar

94. Allen B, Schneider A, Victoria B, Nunez Lopez YO, Muller M, Szewczyk M, et al. Blood serum from head and neck squamous cell carcinoma patients induces altered MicroRNA and target gene expression profile in treated cells. Front Oncol. (2018) 8:217. doi: 10.3389/fonc.2018.00217

PubMed Abstract | CrossRef Full Text | Google Scholar

95. Chan PP, Lowe TM. GtRNAdb: a database of transfer RNA genes detected in genomic sequence. Nucleic Acids Res. (2009) 37:D93–7. doi: 10.1093/nar/gkn787

PubMed Abstract | CrossRef Full Text | Google Scholar

96. Chan PP, Lowe TM. GtRNAdb 2.0: an expanded database of transfer RNA genes identified in complete and draft genomes. Nucleic Acids Res. (2016) 44:D184–9. doi: 10.1093/nar/gkv1309

PubMed Abstract | CrossRef Full Text | Google Scholar

97. Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, Kokocinski F, et al. GENCODE: the reference human genome annotation for The ENCODE Project. Genome Res. (2012) 22:1760–74. doi: 10.1101/gr.135350.111

PubMed Abstract | CrossRef Full Text | Google Scholar

98. Barabási A-L, Gulbahce N, Loscalzo J. Network medicine: a network-based approach to human disease. Nat Rev Genet. (2011) 12:56–68. doi: 10.1038/nrg2918

PubMed Abstract | CrossRef Full Text | Google Scholar

99. Anastasiadou E, Jacob LS, Slack FJ. Non-coding RNA networks in cancer. Nat Rev Cancer. (2017) 18:5. doi: 10.1038/nrc.2017.99

PubMed Abstract | CrossRef Full Text | Google Scholar

100. Cava C, Colaprico A, Bertoli G, Graudenzi A, Silva TC, Olsen C, et al. SpidermiR: an R/bioconductor package for integrative analysis with miRNA data. Int J Mol Sci. (2017) 18:274. doi: 10.3390/ijms18020274

PubMed Abstract | CrossRef Full Text | Google Scholar

101. Demchak B, Hull T, Reich M, Liefeld T, Smoot M, Ideker T, et al. Cytoscape: the network visualization tool for GenomeSpace workflows. F1000Res. (2014) 3:151. doi: 10.12688/f1000research.4492.2

PubMed Abstract | CrossRef Full Text | Google Scholar

102. Backman T, Girke T. systemPipeR: NGS workflow and report generation environment. BMC Bioinformatics. (2016) 17:388. doi: 10.1186/s12859-016-1241-0

CrossRef Full Text | Google Scholar