Distinct Profiles for Mitochondrial t-RNAs and Small Nucleolar RNAs in Locally Invasive and Metastatic Colorectal Cancer

Colorectal carcinoma is the third most common cancer and also the third leading cause of cancer-related deaths in the United States (1). The prognosis in advanced cases is poor with more than one third of the patients dying from progressive disease within 5 years (2). However, simultaneous evaluation of small RNA (sRNA), miRNA, and mRNA in matched pairs of colorectal carcinoma by RNA-Seq, a technology for high-throughput sequencing of the whole transcriptome with some advantages over microarrays, has not been reported.

miRNAs are noncoding RNA (ncRNA) molecules of 18 to 25 nucleotides (nts) in length and are known to repress mRNAs by inhibiting translation and stimulating mRNA degradation (3). Approximately 270 miRNAs have been shown to be dysregulated in colorectal carcinoma (3–5) by real-time PCR and microarray. sRNAs are approximately 60 to 300 nts in length, which include transfer RNAs (tRNA) and small nucleolar RNAs (snoRNA). Mitochondrial transfer RNAs (MT-tRNA), approximately 70 to 80 nts long, are structurally distinct from nuclear-encoded tRNAs by having one less loop than the nuclear canonical cloverleaf structure due to genome economization (6). Twenty-two MT-tRNAs and two ribosomal RNAs (rRNA) are required for the synthesis of 13 of the 84 essential protein subunits of NADH dehydrogenase, cytochrome reductase, cytochrome oxidase, and ATP synthase, which are essential for ATP production through oxidative phosphorylation in the mitochondria (7). Upregulated expression of such MT-tRNAs has been found in breast cancer (8). snoRNAs play a crucial role in ribosome biogenesis through methylation and pseudouridylation of rRNAs (9) and a linkage between snoRNAs and carcinogenesis has been established for non–small cell lung cancer (10). mRNAs are RNA templates for protein synthesis, and thus far, more than a dozen mRNAs have been identified as biomarkers for staging and prognosis of colorectal carcinoma (11).

Our study sought to investigate, using stringently collected and preserved tissue samples (12) and RNA-Seq analysis, whether expression of unique sRNA species correlated with and potentially contributed to the aggressive biologic behavior of advanced-stage colorectal cancers and cancers with known metastases, as compared with lower-grade tumors and healthy tissues.

Thirty-five paired tissues of pretreatment colorectal carcinoma were collected by Indivumed GmbH. Histologically, tumor content is 50% to 70% in tumors and 0% in normal tissues. Normal tissues were 5 cm away from tumors. Ischemia time to freeze was 6 to 11 minutes. The normal mucosa was collected from a site at a minimum of 5 cm from the tumor margin. Clinical and histopathologic characteristics of the patients are summarized in Supplementary Table S1. Written informed consent was obtained from all the participants involved in the study via the procuring laboratory at Indivumed.

All 35 pairs were sequenced at small RNA level, and the first 15 pairs were also sequenced at the miRNA level. The RNA quality was assessed using the Agilent 2100 Bioanalyzer. Samples with an RNA Integrity Number of 7 or higher were processed to generate libraries for small RNA sequencing following the Illumina TruSeq Small RNA Sample Preparation protocols. In brief, 3′ and 5′ RNA adapters, specifically modified to target the 3′ hydroxy group and 5′ phosphate group of most small RNA molecules, were ligated to the ends of sRNAs present in 1 μg of high quality total RNA. Reverse transcription was performed to generate the first-strand cDNAs followed by synthesis of the second-strand cDNAs. The small RNA libraries were loaded onto a 4% agarose gel after 11-cycle PCR amplification. DNA fragments of 145 to 160 bp containing miRNA inserts were excised from the gel and purified for micro RNA sequencing (miRNA-Seq), whereas those of 160 to 400 bp containing longer inserts were similarly purified for small RNA (including MT-tRNAs and snoRNAs) sequencing (sRNA-Seq). The yield of the small RNA library was quantified on an Agilent 2100 Bioanalyzer. Twelve to 24 small RNA libraries were pooled, denatured, and loaded onto one lane of a flow cell for cluster generation using the Illumina cBot. The flow cell was loaded onto an Illumina HiSeq 2500 sequencer and subjected to single-end, 50- and 100-cycle sequencing for miRNA and sRNA sequencing, respectively.

Fourteen pairs were further processed to generate mRNA-Seq libraries using an Illumina TruSeq “stranded mRNA sample prep kit.” In this method, poly (A) tailed RNA was purified from 0.5 μg total RNA, fragmented, and reverse-transcribed into cDNAs. Double-strand cDNAs were adenylated at the 3′ ends and ligated to indexed sequencing adaptors, followed by limited-cycle (13) amplification. Paired-end sequencing was carried out on HiSeq 2500 sequencer (Illumina) for 100 × 2 cycles.

The miRNA-Seq data analysis was performed on the NIH Biowulf supercomputer. miRDeep2 was used to trim the adapter sequences, map the reads to the miRNA database, miRBase (version19), and quantify miRNA expression levels. Novel miRNAs were predicted using miRDeeps2 for aligning the reads against the reference human genome (hg19). For mRNA sequencing, Tophat V.2.0.11 and Cufflinks V.2.2.1 were used to align the reads in fastq files to the RefSeq UCSC human hg19 transcript reference genome annotation database, and the quantification of relative abundance of each transcript was reported as Reads Per Kilobase of transcript per Million mapped reads (RPKM). sRNA sequencing data were analyzed using the software CLC Genomics Workbench 5 for the mapping of trimmed reads to Ensemble GRCh37.75 noncoding RNA database and quantitation of the sRNA expression levels. Statistical analyses and gene clustering were carried out using R 2.15.3 (The R Project for Statistical Computing Program). The sRNA expression levels were normalized to the expression levels of all identified, known sRNAs. Paired Student t test was conducted to identify sRNAs with differential expression between cancers and matched normal adjacent tissues (P < 0.05). Hierarchical clustering was conducted using Euclidean distance to cluster RNAs and samples. By standard detection criteria (14), 5,171 sRNAs, 1,009 miRNAs, and 13,580 mRNAs with 1 RPKM in at least one sample were selected for analysis. By paired t-test, 37 sRNAs, 204 miRNAs, and 665 mRNAs were identified to have significant differential expression (P < 0.05). The expression patterns of sRNAs, miRNAs, and mRNAs are presented on clustering analysis, which were performed by the R package “gplots.” Principal component analysis (PCA) was conducted to visualize differences between samples and cluster samples on the gene expression profiles, with the majority of the data variance explained by PC1, PC2, and PC3. Pearson correlation analysis was done with EXCEL 2003, and the correlation coefficient cutoff was 0.6 (r > 0.6). Plots were made with Prism 6 software and EXCEL 2003. Further information is provided in Supplementary Information S1 regarding the verification of RNA-Seq fidelity and identification of colonic tissue housekeeping genes.

To evaluate the sRNA profile of tumors that had documented lymph node (LN) and lymphatic vessel (LV) metastases, the average ratios (tumor/normal: T/N) of sRNAs in tumors with metastases (MET⁺; 9 tumors: T8_LN, T11_LV, T14_LN, T15_LN, T20_LV, T25_LV, T26_LN, T30_LV, and T35_LN; Supplementary Table S1) and tumors without documented metastases (MET⁻; 26 tumors) were compared. We found that 6 MtRNAs [MT-TI (isoleucine), MT-TL (leucine)1, MT-TE (glutamic acid), MT-TP (proline), MT-TF, and MT-TL2] and 4 sRNAs (snord19.2, snord86, snord77, and snora71D) were among the top 10 upregulated genes with 2.5- to 6.5-fold increases in tumors with LN/LV metastasis compared with 1.5- to 3-fold increases in tumors without LN/LV metastases (Fig. 1A). Next, we evaluated whether such genes were also upregulated in advanced-stage versus early-stage tumors by comparing the average ratios (T/N) of sRNAs among tumors at stage 1 (5 tumors), stage 2 (6 tumors), and stage 3 (21 tumors) compared with stage 4 (3 tumors). We found that there was a 6- to 28-fold increased expression of 7 MT-tRNAs [MT-TC (cysteine), MT-TI, MT-TL1, MT-TP, MT-TE, MT-TF, and MT-TS (serine)1] and 3 snoRNAs (snord19.2, snord86, and snord77) in stage 4 tumors compared with 1- to 3-fold increases in all lower stage tumors (Fig. 1B). Overall, MT-tRNAs were found to be highly expressed in high-grade tumors and in tumors with LN/LV metastases. Moreover, the upregulated expression of snord19.2, snord86, snord77, and snora71D in colorectal carcinoma has not been reported, and their potential roles in malignant transformation/progression not known.

Thus to determine whether deregulated MT-tRNAs and snoRNAs can be used as a signature to distinguish some metastatic and high-stage tumors, unsupervised clustering and PCA analysis of the selected 6 top upregulated MT-tRNAs [MT-TF, MT-TL1, MT-TY (tyrosine), MT-TE, MT-TH (histidine), and MT-TS1] and 2 top upregulated snoRNAs (snord43 and snord86) were performed. The hierarchical clustering distinguished MET⁺ T8_LN, T14_LN, T15_LN, T25_LV, T26_LN, T30_LV, and T35_LN (Fig. 1E), whereas the PCA plot distinguished a subset of such MET⁺ tumors, including T8_LN, T14_LN, T15_LN, and T26_LN (Fig. 1G). Among 9 MET⁺ tumors, 3 (T20_LV, T30_LV, and T35_LN) had only upregulation of snord43 and snord86 but not MT-tRNAs. Furthermore, PCA of 22 MT-tRNAs also distinguished metastatic T11_LV, T14_LN, T15_LN, and T26_LN (Supplementary Fig. S2A). Five MET⁻ tumors (T9, T12, T13, T18, T22, and T23) that clustered with MET⁺ tumors (Fig. 1E and G; Supplementary Fig. S2A) were moderately differentiated stage 2 or 3 tumors with high expression of MT-tRNAs (Supplementary Fig. S3A). Thus, patients with tumors bearing such MT-tRNA profiles may be at higher risk for progression and metastasis.

We next evaluated sRNAs that were most prominently downregulated in MET⁺ as well as stage 4 tumors and found that the snord113-114-116 cluster, snord64, snord71, snord93, snord107, scarna18, CTD2651B20.6 (CTD: comparative toxicogenomics database), and AP000318.1 (AP: annotation process) were the 10 most downregulated sRNAs with 1.2- to 10-fold decreases (Fig. 1C and D) Interestingly, the snord113-114-116 cluster, snord116-7, and snord64 are deleted in Prader–Willi syndrome (PWS), a severe metabolic disorder causing excess adiposity, which suggests a role for such dysregulation in energy storage in tumors (15). However, these snoRNAs have not previously been specifically linked to colorectal carcinoma malignant transformation or progression.

Because MET⁺ tumors clustered in multiple analyses of noncoding gene expression, we performed an unsupervised clustering analysis of 14 tumor pairs using 11 known metastatic mRNA biomarkers (13) including the following: liver-intestine cadherin 17 (CDH17), protocadherin 1 (PCDH1), transcription factor 5 (E2F5), matrix metalloproteases 3 and 7 (MMP3 and MMP7), transcription factor AP-4 (TFAP4), transcription factor 1 (E2F1), transforming growth factor, beta 1 (TGFB1), zinc finger transcriptional repressor 1 (SNAI1), and extracellular matrix receptor III (CD44). The resulting cluster analysis distinguished all MET⁺ tumors (T8_LN, T14_LN, T15_LN, and T11_LV) in 14 tumors (Fig. 1F).

Taken together, the coding RNA data support ncRNAs signatures that distinguish MET⁺ colorectal carcinoma. Thus, such MT-tRNA/snoRNA signatures have the potential to more precisely identify these more aggressive and deadly tumors and suggest novel therapeutic targets.

To evaluate the spectrum of sRNA biomarkers associated with tumors regardless of grade, stage, and LN/LV metastases, as distinguished from normal tissue, the average ratios (T/N) of sRNAs in all 35 tumors and all 35 normal samples were compared. We again found a predominance of 9 MT-tRNAs (MT-TI, MT-TL1, MT-TL2, MT-Y, MT-TC, MT-TE, MT-TF, MT-TH, and MT-TD) among the top 20 upregulated genes with 1.4- to 5-fold increases (Supplementary Fig. S2B and Supplementary Table S3A). Nineteen snoRNAs and AP000318.1 were among the 20 most downregulated genes in tumors with approximately 1.5- to 2.5-fold decreases relative to normal tissue (Supplementary Fig. S2C). Overall, 41 sRNAs were found upregulated or downregulated (Supplementary Table S2A–S2C; Supplementary Fig. S3A–S3C). PCA analysis of these 32 differentially expressed snoRNAs distinguished 28 of 35 tumors from normal controls (Supplementary Fig. S2D). Among 41 differentially expressed sRNAs (9 MT-tRNAs and 32 snoRNAs), we found that the upregulation of snord12B expression in tumors bore the greatest statistically significant difference from normal samples (P = 3.830e–12; Supplementary Fig. S4A). Consistent with upregulation of snord12B, the host gene of this snoRNA, zinc finger nuclear transcription factor x-box binding 1 anti-sense 1 (ref. 16), a long non-coding RNA, was upregulated in all examined colorectal carcinomas (Supplementary Fig. S4B). There was an association strength (r = 0.66) between the expression of snord12B and ZNFZ1 AS1 (Supplementary Fig. S4C). Paradoxically, the expression of ZNFZ1-AS1 is downregulated in mouse mammary tumors, and thus it may have distinct functions in different tissues or species. Although, snord12B was reported as upregulated in rectal cancer (17), the profound differences in snoRNA expression between colorectal carcinoma and normal mucosa have not hitherto been reported and suggest potentially important roles for these ncRNAs in malignant transformation, progression, or metastasis in colorectal carcinoma.

Key biologic premises that underlie relationships of genes regarding “drivers” and “passengers” consist of the following: highly coexpressed genes are more likely to be coregulated, and those genes that display prominent connectivity patterns tend to play biologically influential or regulatory roles in disease-related processes (18–20). Measures of node centrality in biologic networks may detect genes with critical functional roles. In gene coexpression networks, highly connected genes (i.e., candidate hubs) have been associated with key drivers of disease pathways, and gene connectivity has been shown to be a measure of functional relevance (18–20). Thus to explore whether there were potential cancer “driver” sRNAs, we performed correlational analyses among 41 differentially expressed sRNAs (9 MT-tRNAs and 32 snoRNAs), 204 miRNAs, and 665 mRNAs in 14 tumor pairs. Interestingly, we found that expression of some sRNAs, such as MT-TF, snord12B, and the snord114 cluster, correlated with expression of more than 100 to 450 miRNAs/mRNAs, whereas expression of some sRNAs, such as MT-TY, snord19B, and snord83.9, correlated with expression of less than 10 miRNAs/mRNAs (Table 1A and B and Supplementary Information S2). Moreover, the upregulated expression of MT-TF and snord12B sRNAs positively correlated with upregulated hypoxia, metastasis, and pentose phosphate pathways (PPP) genes, such as mir21, mir181, ECE2, PHF19, and transketolase (TKT), and negatively correlated with expression of tumor suppressor genes, such as let7c, mir139, CGRRF1, and SRSF5 (Fig. 2A–D). The downregulation of snord114-1 correlated with downregulation of tumor suppressors, such as let7c and AHNAK (Supplementary Information S2D). In addition, MT-TF, snord12B, and snord114-1 were commonly deregulated in all 35 tumors (Supplementary Table S2A–S2C and Supplementary Fig. S3A–S3C). Thus, our data suggest that deregulation of these 3 sRNAs may be causally linked to tumorigenesis, tumor proliferation, or aggressiveness. The regulatory mechanisms controlling expression of these 3 sRNAs are not known.

Because the prominent upregulation of MT-TF correlated with hypoxia biomarkers (mir21 and mir31; Fig. 2A), we examined whether enzymes to mitigate hypoxia were also affected. The glutathione antioxidative pathway consisting of 15 enzymes pertaining to glutathione S-transferase (GSTA), glutathione peroxidase (GPX), gamma-glutamyltransferase (GGT), glutathione reductase (GSR), microsomal/glutathione S-transferase (MGST), and catalase (CAT) in 14 colorectal carcinomas was examined, and downregulation of 9 out of 15 antioxidative enzymes (Fig. 3A) was found in all tumors as well as in 2 normal controls (N8 and N14), which were associated with MET⁺ tumors. Furthermore, we also examined hypoxia-induced miRNAs (21–23) and found 13 upregulated miRNAs pertaining to hypoxia in 14 of 15 tumors (Fig. 3B), the sole exception being a poorly differentiated stage 3 tumor. Because IL8 is known to be upregulated and glycerol-3-phosphate dehydrogenase 1-like (GPD1L) is known to be downregulated by hypoxia (24), we examined the mRNA expression profile of these 2 genes in the 14 tumor pairs to establish the validity of the hypoxia clustering analysis. Relative to healthy control tissue, IL8 was upregulated 2- to 35-fold, and GPD1L was downregulated 1.2- to 5.2-fold in all tumors examined (Fig. 3C). Thus, these data are indicative of hypoxia in colorectal carcinomas. Because hypoxia switches ATP production from oxidative phosphorylation to glycolysis in tumor cells (25–27), we analyzed the expression of 3 SIRTs, which are known to suppress glycolysis (28) and found that all tumors had downregulation of all 3 SIRTs (Fig. 3D). Further, a clustering analysis of 6 genes directly pertaining to glycolytic enzymes, including glyceraldehyde-3-phosphate dehydrogenase (GAPDH), glucose phosphate isomerase (GPI), pyruvate kinase2 (PKM2), enolase1 (ENO1), phosphofructo-2-kinase/fructose-2,6-biphosphatase 3 (PFKFB3), and phosphoglycerate kinase 1 (PGK1), clustered all 14 tumors together with upregulation of all 6 genes, except in one stage 1 and two stage 3 tumors (Fig. 3E), whereas a clustering analysis of 7 genes directly pertaining to TCA cycle enzymes, including malate dehydrogenase 1 (MDH1), oxoglutarate (alpha-ketoglutarate) dehydrogenase (OGDH), aconitase 2 (ACO2), succinate dehydrogenase (SDHD), citrate synthase (CS), fumarate hydratase (FH), and isocitrate dehydrogenase 1 (IDH1), grouped 12 of 14 tumors together (Fig. 3F). Among these 12 tumors, 6 high-stage (3 or 4) tumors did not have downregulation of all 7 TCA cycle enzymes. One poorly differentiated stage 3 tumor (T4) and one poorly differentiated stage 4 tumor (T8) with LN metastases grouped with normal samples in this analysis. These two tumors showed the most upregulated mRNA expression of CS, a rate-limiting enzyme in the TCA cycle, suggesting that although the glycolytic pathway is generally upregulated and the TCA pathway is generally downregulated in colorectal carcinoma, the most aggressive tumors may use both glycolysis and the TCA cycle for energy production. Alternatively, one pathway or the other may predominate in different cells within the heterogeneous tumor population or in different environmental niches. Furthermore, all 14 examined tumors also had evidence of upregulation of the pentose phosphate pathway as shown by upregulation of enzymes, including TKT, glucose 6-phosphate dehydrogenase (G6PD), and 6-phosphogluconate dehydrogenase (PGD); Fig. 3G). The pentose phosphate pathway utilizes glucose to generate ribonucleotides as well as NADPH and plays a pivotal role in anabolic processes and in combating oxidative stress in glycolytic cancer cells. Overall, these data suggest that there is increased utilization of glucose in colorectal carcinomas via anaerobic glycolysis. Intriguingly, one normal sample (N8), paired to MET⁺ T8, was adjacent to tumors in various clustering analyses (Figs. 1F and 3A, D, and F). Compared with other normal samples, this N8 had higher expression of oncogenes, such as SNAI1 and CD44, as well as lower expression of antioxidative enzymes and tumor suppressors, such as GST, GSR, CAT, SIRT3, and SIRT6 (Figs. 1F and 3A and D) at the mRNA level. These data suggest the potential presence of metastasis in this supposedly normal sample despite its designation as “normal” by histologic diagnosis (Supplementary Table S1).

Finally, to compare and contrast expression profiles in tumors across the histologic spectrum, we established individual tumor sRNA fingerprints in paired matched samples and evaluated commonality and dissimilarity in sRNA expression. Because both the magnitude of the ratio (T/N) and of the subtractive difference (T–N) in gene expression may have biologic significance, we established both relative ratio (T/N) and absolute difference (T–N) fingerprints for each tumor. We then selected the three most highly up- and downregulated sRNAs based on values of (T/N) or (T–N) to build an sRNA fingerprint for each tumor. Relative ratio fingerprints identified the MT-tRNA family, AP000318.1, snord19, CTD2651B20.6, snord64, and snord19 cluster as the most frequently upregulated genes, whereas the snord114 cluster, AP000318.1, and snord64 were identified as the most frequently downregulated genes (Table 2, and Supplementary Table S3A), largely reflecting earlier analyses (Supplementary Fig. S2B and S2C). As in the initial analysis of average fold increased or decreased expression of each sRNA in all MET⁺ tumors (Fig. 1A and C), the 3 most upregulated (8–20 fold) genes in 9 individual metastatic tumors were also either MT-tRNAs or the snord19 family, whereas the 3 most downregulated (2–25 fold) genes again were in the snord113/114/116 cluster (Supplementary Table S3A). For absolute difference fingerprints, we found that in metastatic tumors, snord43, MT-tRNAs, and the snord12 cluster were the most frequently upregulated genes, whereas snord71, snord59A, and the snord114 cluster were the most frequently downregulated genes among the top 3 up/downregulated genes (Table 3 and Supplementary Table S3B). Met⁺ tumors (T8_LN, T15_LN, T25_LV, and T35_LN) displayed the greatest level of absolute RPKM changes of snord43, snord12 cluster, and snord113/114/116 cluster expression. The deregulation of these snoRNAs in colorectal carcinoma has thus far not been reported, and their roles, if any, in tumorigenesis, proliferation, or metastasis are not known. Overall, each tumor had its own genomic signature, with some unique features and some common features. Based on the presence or absence of MT-tRNAs in the top 3 upregulated genes, 35 colorectal carcinomas can be placed into two groups, regardless of their histologic types and TNM classification (Table 4).

Colorectal carcinoma is a disease with variegated genetic and epigenetic profiles (29). Advanced molecular profiling of colorectal carcinoma may improve staging and grading of tumors, thus better predicting clinical course, as well as defining specific tumor survival and proliferation pathways, thereby elucidating novel treatment strategies. Currently, there are a number of drugs that target such specific oncogenic pathways (30). In this study, we delineated sRNA for 35 colorectal carcinomas. We found that all 35 colorectal carcinomas had distinct combinations of common as well as uniquely deregulated genes. This new understanding of oncogenic mechanisms may influence risk assessment, diagnostic categories, and therapeutic strategies, with increasing use of drugs and antibodies designed to counter the influence of specific molecular drivers.

Our study mainly focused on distinctions among tumors of high stage and tumors with metastases. In this study, we identified MT-tRNAs as the top upregulated sRNAs in LN/LV-positive and high-stage colorectal carcinoma, as well as a signature of sRNAs associated with such tumors and, in some cases, with their adjacent “normal” tissues which likely contain undetected metastasis. Furthermore, the distinct profiles of tumors with metastasis, as defined by both ncRNA and coding RNA signatures in colorectal carcinoma, revealed some novel elements. Especially notable was the stunning downregulation of the snord113/114/116 clusters that are deleted in PWS, a genetically imprinted disease marked by profound obesity, which may likewise enhance energy storage in cancer cells, but not described thus far in cancer. Whether these snoRNA clusters having tumor suppressor function needs to be studied.

A multitude of data support the role of hypoxia in tumor physiology as hypoxia has been shown to upregulate glycolytic metabolism and MT-tRNAs (31) and downregulate enzymes involved in the TCA cycle. More generally, hypoxia has long been recognized to play a role in promoting the invasive and metastatic behavior of cancer cells, and their ability to disseminate from established colonies and establish new colonies in distinct environments may depend on their ability to diversify metabolic pathways (21, 23). In this regard, the stringency of tumor harvest and preservation in our study is critical in assuring that these findings are not artifacts of prolonged ex-vivo time duration prior to storage (12). The superposition of upregulated MT-tRNAs, which may support generation of ATP through mitochondrial oxidative phosphorylation, on a profile of a shift to glycolytic metabolism, suggests that tumors with enhanced capacity for local and distant invasion are poised to use either oxidative phosphorylation or glycolytic metabolism, depending on the oxygenation circumstances in which they find themselves. Though downregulation of snoRNAs, such as snord50A, h5sn2, snord43, and snord44, has been found in breast, prostate, brain, and lung cancers (32, 33), the panel of downregulated snoRNAs, identified in this study, has not been reported previously and may be unique to colonic cancers. Interestingly, the host gene of the downregulated snord114 cluster (Supplementary Fig. S5A) is the maternally expressed 3 RNA gene, MEG3 (34, 35). MEG3 inhibits cancer cell proliferation by both p53-dependent and p53-independent pathways (34, 35). Consistent with this finding, our sRNA/mRNA sequencing data revealed a strong expression correlation (downregulation of expression) between snord114-23/snord114-26 and MEG3 (rs > 0.8; Supplementary Fig. S5B). Whether these imprinted snoRNAs have specific biologic functions distinct from their host gene in malignant transformation needs to be further evaluated.

In conclusion, the 35 paired samples analyzed in this study serve as a training set that helps to identify sRNA signatures for metastatic and advanced colorectal carcinoma. To strengthen our findings, these putative sRNA markers will be evaluated and potentially validated by an independent study with larger numbers of samples. The hope is that such ncRNAs and coding RNAs will be validated as biomarkers and therapeutic targets which will ultimately benefit patients.

No potential conflicts of interest were disclosed.

Conception and design: L. Xu, H. Juhl, A. Rosenberg

Development of methodology: L. Xu, W.W. Wu, H. Juhl, R. Wang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): W.W. Wu, R.-F. Shen, H. Juhl, R. Wang

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): L. Xu, J. Ziegelbauer, R. Wang, W.W. Wu, H. Juhl

Writing, review, and/or revision of the manuscript: L. Xu, J. Ziegelbauer, R. Wang, W.W. Wu, R.-F. Shen, H. Juhl, A. Rosenberg

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L. Xu, H. Juhl, Y. Zhang, R. Wang

Study supervision: A. Rosenberg

The authors thank Julia Berkson, Baolin Zhang, and Ashutosh Rao for critically reviewing and commenting on this article.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.