Characterization of Transcription Termination-Associated RNAs: New Insights into their Biogenesis, Tailing, and Expression in Primary Tumors

Next-generation sequencing has uncovered novel classes of small RNAs (sRNAs) in eukaryotes, in addition to the well-known miRNAs, siRNAs, and piRNAs. In particular, sRNA species arise from transcription start sites (TSSs) and the transcription termination sites (TTSs) of genes. However, a detailed characterization of these new classes of sRNAs is still lacking. Here, we present a comprehensive study of sRNAs derived from TTSs of expressed genes (TTSa-RNAs) in human cell lines and primary tissues. Taking advantage of sRNA-sequencing, we show that TTSa-RNAs are present in the nuclei of human cells, are loaded onto both AGO1 and AGO2, and their biogenesis does not require DICER and AGO2 endonucleolytic activity. TTSa-RNAs display a strong bias against a G residue in the first position at 5′ end, a known feature of AGO-bound sRNAs, and a peculiar oligoA tail at 3′ end. AGO-bound TTSa-RNAs derive from genes involved in cell cycle progression regulation and DNA integrity checkpoints. Finally, we provide evidence that TTSa-RNAs can be detected by sRNA-Seq in primary human tissue, and their expression increases in tumor samples as compared to nontumor tissues, suggesting that in the future, TTSa-RNAs might be explored as biomarker for diagnosis or prognosis of human malignancies.


Introduction
In the last ten years, the technologies for RNA profiling strongly improved. These advances disclosed pervasive transcription of more than 70% of the human genome, with protein-coding genes accounting for less than 2% of the total RNA [1]. Thus, a dominant portion of the transcribed regions on the human genome originates from non-protein-coding genes (noncoding RNAs (ncRNAs)). Increasing evidences are being obtained, indicating that noncoding RNAs possess essential regulatory functions and could be one of the major contributors to the complex traits of the organisms [2].
sRNA profiling by next-generation sequencing (sRNA-Seq) is the method of choice for the identification of lowly abundant sRNA classes, in addition to the miRNA, siRNA, and piRNA families. Indeed, sRNAs with varying lengths of between 18 and 200 nucleotides have been reported to be derived from specific genomic regions in higher eukaryotes. At least three classes of sRNAs derived from regions mapping around the 5 ′ termini of genes have been described: transcription initiation RNAs (tiRNAs) [5], transcription start site-associated RNAs (TSSa-RNAs) [6][7][8], and promoterassociated small RNAs (PASRs) [9,10]. The origin and function of these RNAs are uncertain, but preliminary evidence suggests that they are involved in epigenetic control of gene expression. Recently it has been described a new family of sRNAs, termed DNA damage-response RNAs (DDRNAs) [11] or double-strand break-(DSB-) induced RNAs (diRNAs) [12], that are generated at sites of DNA damage and control the DNA damage response.
Over the years, sRNAs have become the focus of biomarker research, an approach that has been favorably used in the prediction and early detection of disease and in the investigation of response to treatment for several medical conditions. miRNAs are the most frequently assessed for their potential role as biomarkers, such as in cancer [16], liver [17], and cardiovascular disease [18] among many others. Moreover, dysregulation of other sRNA species such as piR-NAs, small nucleolar RNAs (snoRNAs), and small nuclear RNAs (snRNAs) is being found to have relevance to tumorigenesis, neurological, cardiovascular, developmental, and other diseases [19]. Thus, they might be potentially assessed as biomarkers of disease.
Here, we present a comprehensive characterization of TTSa-RNAs [15] in human cells. We show that TTSa-RNAs are present in both AGO1 and AGO2 complexes immunopurified from human cell nuclei. Moreover, TTSa-RNAs can be detected in both soluble and chromatin-associated nuclear extract. Biogenesis of this class of sRNAs is independent of both DICER and AGO2 endonucleolytic activity. TTSa-RNAs show a very strong bias against a G residue in the first position at 5′ end, which is a specific feature of AGO-bound sRNAs and a peculiar oligoA tail at 3′ end. Finally, AGObound TTSa-RNAs derive from genes involved in cell cycle progression regulation and DNA integrity checkpoints. Interestingly, this class of sRNAs can be detected in primary human tissues, and their expression increases in tumor samples as compared to nontumor tissues.
HeLaS3 AGO2KO monoclonal cell line was obtained by using two specific Zinc Finger Nucleases (ZNF1 and ZNF2, CompoZr® Knockout Zinc Finger Nucleases, Sigma). HeLaS3 cells were plated in a 6-well plate and transfected using Lipofectamine® 2000 (Thermo Fisher Scientific) with 1.5 μg of each ZNF expression plasmid. After two sequential transfections with the AGO2-specific ZNF expression plasmids, individual clones were isolated by serial dilution and assayed by Western Blot and RT-qPCR for loss of AGO2 expression.
2.3. RNA Isolation and RT-qPCR. RNA from HeLa S3, HeLa S3 AGO2KO, HCT 116 WT, and DICER Ex5 cells, as well as fractions from nuclear fractionation protocol, was isolated using Direct-zol™ RNA MiniPrep Kit (Zymo Research) according to the manufacturer's protocol.
TTSa-RNA (Table 1) expression was checked by using Custom TaqMan® Small RNA Assay (Thermo Fisher Scientific). RNA was retrotranscribed using SuperScript II Reverse Transcriptase (Thermo Fisher Scientific), and qPCR was performed using iTaq Universal Probes Supermix (Bio-Rad). Quantification was normalized to small nucleolar RNA U44 (for total RNA) or to ath-miR-159a spike-in (for SNE and CPE fractions), amplified by TaqMan Small RNA Assay (Thermo Fisher Scientific).
2.5. Small RNA Sequencing. sRNA-Seq was performed from total RNA isolated from HeLaS3 and AGO2KO HeLaS3 cells by IGA TECH (Udine, Italy). Library was prepared using TruSeq Small RNA Library Preparation Kit (Illumina). Sequencing was performed on HiSeq 2500 in a 50 bp singleread mode 10 million reads/sample. It is worth pointing out that the protocol we applied to generate sRNA libraries only captures 5 ′ phosphorylated RNA molecules and does not allow cloning of sRNAs harboring a 5 ′ Cap.
2.6. TTSa-RNA Data Analysis from Cell Lines. RNA Immunoprecipitation and next-generation sequencing were performed as previously published [6,22]. NGS data quality was checked using FastQC. tRNA and rRNA genomic coordinates on hg38 genome were retrieved from UCSC Repeat-Masker track [23] using Table Browser tool [24]. By using bedtools [25], tRNA and rRNA sequences in FASTA format were retrieved from hg38 genome sequence and used to build a bowtie1 [26] reference index for tRNA and rRNA. miRNA hairpin FASTA sequence was downloaded from miRbase (v21) [27], and human hairpins were extracted. Using bowtie-build, a bowtie1 index was obtained for miRNA hairpins.
Each FASTQ file was processed as follows: (1) Adapter sequences were removed using cutadapt; only reads longer than 17 nt were retained.
(3) Reads without any valid alignment in step 2 were aligned to the abovementioned tRNA and rRNA bowtie1 index with -n 2 -l 18 options.
(4) Reads without any valid alignment in step 3 were aligned to hg38 human reference with bowtie1 with -n 0 -l 18 -m 1 options. Only reads with a single valid alignment were investigated further.
TTSa-RNAs were defined as those reads aligned in step 4, with at least 1 nt overlapping a 5 nt-long window centered on genomic coordinates of polyA sites (GENCODE v24, [26]). Only polyA sites with a 1 nt length were taken into account in all the analyses (99.8% of the 48535 polyA sites reported by GENCODE v24 annotation).
Within each sample, RPM was defined as the number of reads divided by the total number of reads that gave at least one valid alignment in either step 2, 3, or 4 of the alignment procedure * 10 6 .

PAR-CLIP Data
Analysis. PAR-CLIP libraries were analyzed using the same pipeline outlined above.

Read Coverage Plots.
Cumulative coverage of sRNAs giving rise to a single alignment in step 4 was computed using bedtools in a window of 1000 nt centered on GENCODE v24 polyA sites. Each sRNA contributed to the coverage of the closest polyA site only. Three polyA sites (# 605212, 627827, and 582464) lying within 500 nt of loci giving rise to very large numbers of small RNA molecules (not-overlapping any polyA site) were omitted to obtain a clearer picture.
For PAR-CLIP datasets, only reads harboring a single mismatch to the reference genome were used to compute coverage.
2.9. GO Term Analysis. Ensembl GENE IDs were retrieved from GENCODE v24 annotation. Genes giving rise to TTSa-RNAs were defined as any Ensembl GENE overlapping a single TTSa-RNA. GO term analysis was performed on DAVID [29] by submitting the list of Ensembl GENE IDs overlapped by TTSa-RNA in each sample. As background, the entire set of Ensembl GENE IDs was used. The reported GO terms were obtained by selecting GO_term_BP_4 in the DAVID output and downloading the full table.
Significantly enriched GO terms (FDR < 0.01) were sorted for decreasing "Fold Enrichment" and the Fold Enrichment was plotted for the top 10 GO Terms.

Correlation of mRNA Expression with TTSa-RNA
Counts. mRNA expression in HeLaS3 was retrieved from ENA: EXR352930, EXR352926. Reads were pseudoaligned to the human transcriptome, and transcripts were quantified in transcript per million (TPM) using Kallisto v 0.43.1 [30]. Ensembl GENE IDs corresponding to the Ensembl TRANSCRIPT IDs were retrieved using BioMart, and the TPM for all transcripts of each gene were summed, giving a summarized TPM value for each gene. The genes with TPM below the first tertile were filtered out. mRNA expression in 44 adrenocortical carcinoma primary samples was retrieved from GEO repository (GSE49278) and belongs to the same samples whose sRNAseq were used for TTSa-RNA quantification as described below. The Affymetrix Human Gene 2.0 ST Array was analyzed in R (R package version 1.58.1) using Bioconductor packages oligo [31] and genefilter (using RMA as a background correction/normalization/summarization method). Genes with RMA expression values below the first tertile in more than 10 samples were filtered out; the PROBEIDs were mapped to Ensembl GENE IDs and the Ensembl GENE IDs mapped more than once (<10%) were excluded. mRNA expression in 10 kidney tumor samples and 10 matched controls was retrieved from SRA project SRP003901 and belongs to the same samples whose sRNAseq were used for TTSa-RNAs quantification. The digital gene expression (DGE) sequencing reads were trimmed from Illumina adapters using cutadapt [32], resulting in 17 nt reads, which were mapped using bowtie1 with -p 8 -v 0 -m 2 options, after adding the four nucleotides of NlaIII restriction site to the 5 ′ of each read. The DGE tags belonging to each transcript were quantified using htseq-count [33] (with intersection-strict option) and the htseq-counts were normalized using DESeq package [34]. Genes with no counts in more than the 10% of samples and genes with normalized counts below the first tertile were filtered out.
The correlations were computed using custom R scripts and the R 2 was computed from a Spearman correlation coefficient. In all cases, the genes with TTSa-RNA RPM in the first three quartiles were filtered out, because this corresponded to less than 4 TTSa-RNA reads per sample.
2.11. TTSa-RNAs Data Analysis from Primary Samples. A search in SRA was performed using miRNA and cancer/ tumor as keywords. We retrieved 17 datasets containing sRNA-seq data for tumor tissue and matched controls. However, only 5 datasets were suitable for our analysis. In fact, we discarded a dataset that was obtained in colorspace (it could not be analyzed using our computational pipeline), a dataset with a sudden drop in Phred quality (median Phred score = 2 for all positions > 16 in all samples), a dataset which contained only one single healthy tissue control sample, three datasets because all reads were shorter than 18 nt (multiple alignments on the genome did not allow specific identification of TTSa-RNAs), and six datasets in which abundant mRNA fragments with random size distribution were retrieved among sRNA-seq reads, suggesting that mRNA degradation does not allow proper identification of bona fide TTSa-RNAs.
The 5 datasets analyzed correspond to SRA projects: SRP028291, SRP014142, SRP045645, SRP003902, and SRP048750. Analysis of these datasets was performed as outlined above for AGO protein-bound sRNAs, except for the application of an 18 to 26 nt size filter as discussed in the results and discussion section. All p values were computed using Wilcoxon sum-rank test (a paired version of the test was used for datasets containing paired normal tissue and tumor tissue from each patient) using R function wilcox.test (https://www.r-project.org/).

Results and Discussion
3.1. Nuclear AGO1 and AGO2 Bind a Class of sRNAs Arising from Human Transcription Termination Sites. To obtain an overview of nuclear sRNA classes bound to AGO proteins in human cells, we reanalyzed a set of sRNA-seq libraries generated by immunopurifying AGO1 and AGO2 from nuclear extracts of HeLaS3 cell line [6]. We focused on less characterized classes of AGO bound sRNAs, that is reads not aligning to known human miRNAs (miRbase 21) or human tRNAs and rRNAs. All the reads surviving to this preliminary filtering were aligned to human hg38 reference genome requiring unambiguous mapping to a unique genomic position. Reads with multiple mappings on the hg38 reference genome were not further investigated (see Materials and Methods for details). We found that in nuclei of human cells, both AGO1 and AGO2 bind a class of sRNAs (21-24 nt long, Figure 1(a)) arising from TTSs of human mRNA transcripts and lying on the sense strand of the gene they arise from (Figure 1(b)). Notably, 3′ end of these sRNAs in most cases (60% to 66% depending on the dataset analysed) maps within 2 nt of GENCODE-annotated polyA sites. We assumed that these sRNAs belong to the class of TTSa-RNAs. Indeed, TTSa-RNAs have been described as AGO1/2-associated, miRNA-(22-24 nt) sized sRNAs whose 3′ends aligned with the annotated polyadenylation sites and mapping on the sense strand [15]. Another class of sRNA mapping around 3 ′ termini of genes has been described in human cells and in plants, namely TASRs [10]. However, human TASRs are more heterogeneous in size as compared to TTSa-RNAs. Moreover, TASRs, whose 3 ′ termini aligned with annotated polyadenylation sites, map antisense to mRNA transcription. In plants, TASRs are 23-24 nt in length and are bound by an AGO protein [13], similar to TTSa-RNAs. However, like human TASRs, plant TASRs seem to be a heterogeneous class mapping on the sense, on the antisense strand or on both strands of protein coding genes. Although we cannot exclude that TTSa-RNAs are a subclass of TASR, because of . Synthetic ath-miR159a was added to SNE and CPE fractions before RNA isolation and subsequently used for qRT-PCR normalization.
differences in size, orientation, and mapping with respect to polyadenylation sites, we believe that TTSa-RNAs are a distinct class of sRNAs. In our experimental design [6], we attained a sequencing depth of 100 million reads per sample, which is substantially higher than the depth employed by Valen and colleagues [15]. This prompted us to investigate peculiar features of this class of sRNAs in order to return a detailed characterization.
The estimated abundance of TTSa-RNAs is about 100 reads per million (RPM) in libraries generated from the nuclear fraction of sRNAs in HeLaS3, reaching about 300 RPM and 150 RPM in AGO1 and AGO2immunoprecipitated samples, respectively (Figure 1(b)). As shown in Figure 1(c), TTSa-RNAs arise from 2822 genes, and 37% of these genes give rise to TTSa-RNAs that are loaded on both AGO1 and AGO2. More interestingly, GO term analysis suggests that both genes giving rise to AGO1-and AGO2-bound TTSa-RNAs are significantly enriched for genes involved in cell cycle progression regulation and DNA integrity checkpoints (Supplementary Figure 1). A GO term analysis of genes giving rise to IgG-immunoprecipitated sRNAs lying on TTS did not highlight any significant enrichment for the same GO terms (data not shown).
A parallel analysis on a sRNA-seq library obtained by immunopurification of a FLAG-tagged AGO2 complexes from a HEK293 cell line expressing a FLAG-tagged AGO2 [22] yielded similar results, highlighting that TTSa-RNAs are found in different human cell lines (Figures 1(a)  and 1(b)).
Not only TTSa-RNAs show a miRNA-like size that is a specific feature for AGO loading but we also evidenced that AGO2 and AGO1 physically interact with this class of sRNA (Supplementary Figure 1A and B) by taking advantage of PAR-CLIP-Seq datasets [35]. Indeed, PAR-CLIP shows that TTSa-RNAs are associated with AGO2 and AGO1. Interestingly, the abundance of TTSa-RNAs in PAR-CLIP libraries is higher compared to the one observed in noncrosslinked AGO IPs, suggesting that TTSa-RNAs are directly associated with AGO proteins. Notably, only reads containing one single mismatch with reference (mainly T->C transitions, Supplementary Figure 1C and D) were used to generate Supplementary Figure 1A and B, thus ensuring that only bona fide AGO crosslinked molecules are contributing to the signal observed.
Furthermore, we asked whether TTSa-RNA expression levels are correlated with the abundance of the corresponding mRNA. We looked at the correlation between gene expression levels (as assessed by RNA-seq) of each gene and the abundance of the TTSa-RNAs arising from TTSs of the same gene in HeLaS3 cells for both AGO1 (R 2 <0.1) and AGO2 (R 2 <0.1, n = 2). The R 2 values observed do not support any correlation between gene expression and TTSa-RNA abundance.
Differently from Valen and colleagues [15], we showed that TTSa-RNAs can be detected in libraries obtained from nuclear extracts of HeLaS3 cells. We also investigated whether TTSa-RNAs are recruited on chromatin. We extracted HeLaS3 nuclei to separate soluble and loosely bound material (soluble nuclear fraction (SNF)) from the chromatin pellet (chromatin pellet extract (CPE)), which retains tightly bound factors [21] (Figures 1(d) and 1(e)). As a control of a chromatin-associated RNAs, we measured MALAT [36]. Moreover, in order to provide evidence that this procedure also preserves chromatin localization of sRNAs, we verified that the small nucleolar RNA RNU44 was enriched in the CPE as compared to SNF. On the other hand, we included in the analysis a soluble nuclearenriched mRNA, beta-actin (ACTB), and miR-21, which has been previously detected not only in the cytosol but also in the nucleus of human cells [37] (Figure 1(e)). We selected 5 abundant TTSa-RNAs arising from genes belonging to the enriched GO categories (YWHAG, WEE1, PSMB1, CDKNA1, and BUB3, see Table 1). We measured their abundance in different nuclear fractions by qRT-PCR using specific primers designed to selectively amplify sRNA molecules and not their precursors [38]. TTSa-RNAs display a profile of nucleoplasm/chromatin abundance very similar to the one of ACTB mRNA and miR-21 rather than MALAT and U44 (Figure 1(e)), arguing against their specific recruitment on chromatin.

DICER Is Not Involved in TTSa-RNA Biogenesis.
We next focused on TTSa-RNAs biogenesis. Previously, it has been speculated that the miRNA biogenesis pathway was not involved in TTSa-RNA processing [15]. We aimed to experimentally verify this hypothesis and we first investigated whether TTSa-RNAs are processed by DICER. To answer this question, we reanalyzed sRNA-seq libraries derived from sRNAs bound to nuclear AGO2 in HCT116 and HCT116 Dicer EX5 (a subclone of HCT116 cells in which DICER has been impaired) cell lines [20] (Figures 2(a) and 2(b)) that we previously generated [6]. If DICER endonucleolytic activity is required for TTSa-RNA biogenesis, one would expect a dramatic global reduction of TTSa-RNA abundance as it happens for miRNA [20] and for TSSa-RNAs [6]. On the contrary, we found that AGO2-bound TTSa-RNA abundance is considerably increased in HCT116 Dicer EX5 cells compared to parental HCT116 cells (Figure 2(c)). To further validate these findings by qRT-PCR, we analyzed the expression of the same five TTSa-RNAs mentioned above by qRT-PCR in whole cell extracts in HCT116 Dicer EX5 cells as compared to parental cells. None of the five TTSa-RNAs tested by qRT-PCR displayed any decrease in HCT116 Dicer EX5 cells (Figure 2(d)). In agreement with NGS data, qPCR validation showed that for some of the tested TTSa-RNAs, abundance increases in HCT116 Dicer EX5 compared to parental cells. Therefore, in order to exclude that this increase was due to a higher expression of the corresponding gene in HCT116 Dicer EX5 , we took advantage of the GEO dataset GSE6427 which reports mRNA expression levels in HCT116 and HCT116 Dicer EX5 cells. As shown in Figure 2(e), the expression of the YWHAG, WEE1, PSMB1, CDKNA1, and BUB3 mRNAs is not affected by DICER ablation. We hypothesize that higher levels of TTSa-RNAs in DICER hypomorphic cells as compared to control, observed both in NGS data and in qPCR validations, might be explained taking into account that lack of mature miRNA may promote loading on AGO1 and AGO2 of TTSa-RNAs, perhaps increasing their stability. Overall, our data highlights that TTSa-RNA biogenesis is not dependent on DICER.
Since it has been reported that miR-451 is processed in a DICER-independent manner by AGO2 endonucleolytic activity [39], we wondered whether TTSa-RNAs might be processed by AGO2. In case AGO2 was required for TTSa-RNA biogenesis, AGO2 depletion should have resulted in a major decrease of TTSa-RNA abundance in whole cell extracts. We therefore profiled sRNAs in parental and AGO2KO HeLaS3 cells (obtained by genome editing mutation of AGO2 gene, Figure 2(e)) by sRNA-Seq and we calculated TTSa-RNA abundance (Figure 2(g)).
Moreover, we evaluated the abundance of the five selected TTSa-RNAs by qRT-PCR (Figure 2(h)). Both approaches showed a slight decrease in TTSa-RNA abundance when AGO2 is depleted. We concluded that even though AGO2 is not required for TTSa-RNA biogenesis, it might participate to TTSa-RNA stabilization, recapitulating what was previously seen for miRNAs [40].
The fact that both DICER and AGO2 are dispensable for TTSa-RNAs biogenesis is in agreement with our observation that the regions flanking TTSa-RNAs are not biased toward the formation of hairpin structures compared to randomly picked genomic regions (data not shown).

Sequence Characteristics of TTSa-RNAs.
We set out to investigate whether any consensus nucleotide sequence could be found within TTSa-RNAs or in the surrounding genomic sequence. We could not identify a consensus sequence that might provide any hints on TTSa-RNA biogenesis. However, similarly to what was observed for other AGO-bound sRNAs [41,42], TTSa-RNAs have a bias against a G residue in   the first position at 5 ′ end if compared to other AGO2bound sRNAs (Fisher Exact test p value < 2.2e-16) in all libraries generated from AGO-immunoprecipitated samples (Figure 3(a)). We noted that for all libraries, this bias is stronger in TTSa-RNAs (3.25%) as compared to other AGObound sRNAs (8.10%). This can be explained by the low GC content in the regions upstream of human TTSs. We also checked for any tailing at 3′ end of TTSa-RNAs. We found that about 20% to 27% of TTSa-RNAs bound to endogenous AGO1 and AGO2 displayed a short 3′ tail (1 to 3 nt in most cases). Importantly, 3′ tailing of both AGO1and AGO2-bound TTSa-RNAs displays a similar pattern. Furthermore, 41% of the TTSa-RNAs bound to FLAGtagged AGO2 in HEK293T display a 3 ′ tail (Figure 3(b)). Since it has been reported that miRNAs are 3 ′ tailed [43,44], we compared TTSa-RNAs and miRNA 3 ′ tailing in the same samples. In both sRNA classes, we found abundant monoadenylation and monouridylation. On the other hand, our data highlight a significant (p value < 2.2E-16) increase of oligoA tails (4 or more A) in TTSa-RNAs compared to miRNAs. This finding suggests that mature polyadenylated mRNAs might be the precursors of TTSa-RNAs. Although this observation might suggest that TTSa-RNAs are mere degradation by-products of mRNAs, the fact that they arise specifically from the region of mRNAs immediately upstream of polyA sites, the lack of any correlation with corresponding mRNA abundance and the specific loading on AGO proteins argue against this hypothesis. In particular, it is difficult to conceive a general or widespread mRNA degradation pathway leading to the specific accumulation only of degradation fragments flanking polyA sites. Therefore, we hypothesized that these sRNAs are not simply mRNA degradation products but are specific, biologically generated species.

TTSa-RNAs Are Overexpressed in Human Tumors.
Up to now, the existence of sRNAs mapping around 3′ termini of mammalian genes was only reported for cultured cell lines [10,14,15]. Therefore, we extended our analysis of TTSa-RNAs to primary samples, and we set out to look for TTSa-RNAs in sRNA NGS libraries obtained from primary human samples. We searched the SRA repository for datasets in which microRNA profile was assessed by NGS. We retrieved five datasets suitable for TTSa-RNA analysis which contained publicly available sRNA expression profiles in tumor samples and matched control tissue. In all the five datasets analyzed, we could detect bona fide TTSa-RNAs displaying the expect size profiles (Figure 4(a)). A small number of reads mapping to TTS with a size of 35 nt were retrieved; therefore, we only counted molecules shorter than 28 nt as bona fide TTSa-sRNAs. Since AGO-bound TTSa-RNAs in cell lines do not correlate with the expression of the corresponding mRNA, we asked whether bona fide TTSa-RNAs that we identified in primary samples display any correlation with mRNA expression. To answer this question, we took advantage of the mRNA expression profiling available for two of the five datasets. As depicted in Figure 4(b), in both cases, we do not observe any significant correlation, highlighting that not only in cultured cell lines but also in human tissues, the amount of TTSa-RNAs is not correlated with the expression of the corresponding gene. Since deregulation of other sRNA classes has been reported in human tumors [16,19], we looked at TTSa-RNAs expression in tumor tissues and matched controls. In four out of five datasets analyzed, we found an increase of TTSa-RNAs in tumor samples compared to matched healthy tissues (Figure 4(c)). The lack of statistical significance of the result obtained in endometrial cancer dataset is likely due to the small number of samples in this dataset (n = 6). Interestingly, we found that the expression of most TTSa-RNAs is increased in cancer tissues, suggesting a general deregulation of this class of RNAs rather than overexpression of a few TTSa-RNAs (Supplementary Figure 3).
We concluded that TTSa-RNAs are not only detectable in human primary tissues, but their expression also increases in tumors. Accumulation of TTSa-RNAs in tumor samples might reflect the previously reported deregulation of the expression of RNAi components, such as DICER, AGO2, and AGO1 in cancer cells [45]. In particular, downregulation of DICER might impair miRNA processing, thus promoting TTSa-RNA loading on AGO proteins. On the other hand, AGO protein overexpression might stabilize TTSa-RNAs.

Conclusions
Here, we present a comprehensive study of sRNAs derived from TTS of expressed genes in human cell lines and primary tissues. Taking advantage of a sequencing depth of 100 million reads per sample, this characterization gets insights into previously unknown details of biogenesis, localization, sequence features, and expression of this new and poorly studied class. Our data demonstrate that TTSa-RNAs are a class of DICER-independent and AGO-bound sRNAs, which display an oligoA tailing at 3 ′ end and which are expressed not only in cultured cell lines but also in human primary tissues. Even though the function of TTSa-RNAs still remains elusive, the fact that genes that give rise to TTSa-RNAs are involved in regulation of proliferation and DNA damage response and that TTSa-RNAs are overexpressed in human tumors suggest that TTSa-RNA expression is linked to tumorigenesis and they might be explored as biomarker for diagnosis or prognosis of human malignancies in the future.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.