key: cord-0277314-hltr4voz
authors: Yao, Jun; Wu, Douglas C.; Nottingham, Ryan M.; Lambowitz, Alan M.
title: Identification of protein-protected mRNA fragments and structured excised intron RNAs in human plasma by TGIRT-seq peak calling
date: 2020-06-26
journal: bioRxiv
DOI: 10.1101/2020.06.25.171439
sha: 0d590c07f427db14b5f7dfd89acbdb3e9a6c603c
doc_id: 277314
cord_uid: hltr4voz

Human plasma contains >40,000 different coding and non-coding RNAs that are potential biomarkers for human diseases. Here, we used thermostable group II intron reverse transcriptase sequencing (TGIRT-seq) combined with peak calling to simultaneously profile all RNA biotypes in apheresis-prepared human plasma pooled from healthy individuals. Extending previous TGIRT-seq analysis, we found that human plasma contains largely fragmented mRNAs from >19,000 protein-coding genes, abundant full-length, mature tRNAs and other structured small non-coding RNAs, and less abundant tRNA fragments and mature and pre-miRNAs. Many of the mRNA fragments identified by peak calling correspond to annotated protein-binding sites and/or have stable predicted secondary structures that could afford protection from plasma nucleases. Peak calling also identified novel repeat RNAs, miRNA-sized RNAs, and putatively structured intron RNAs of potential biological, evolutionary, and biomarker significance, including a family of full-length excised introns RNAs, subsets of which correspond to mirtron pre-miRNAs or agotrons.

Extracellular RNAs in human plasma have been avidly pursued as potential biomarkers for cancer and other human diseases (1-7). In healthy individuals, plasma RNAs arise largely by apoptosis or secretion from cells in blood, bone marrow, lymph nodes, and liver, while in cancer and other diseases, plasma RNAs may arise by necrosis or secretion from tumors or other damaged tissues, potentially providing diagnostic information (8) (9) (10) (11) (12) (13) (14) . As plasma contains active RNases (15, 16) , the extracellular RNAs that persist there are thought to be protected from degradation by bound proteins, RNA structure, or encapsulation in extracellular vesicles (EVs).

Although high-throughput RNA-sequencing (RNA-seq) has identified virtually all known RNA biotypes in human plasma, studies aimed at identifying disease biomarkers have focused mostly on plasma mRNAs or miRNAs. mRNAs in plasma and blood have been profiled by RNA-seq methods that enrich for or selectively reverse transcribed poly(A)-containing RNAs (10, 17) or by sequencing mRNA fragments with (13, 14, 18, 19) or without (4, 20) size selection, and it remains unclear which methods might be optimal for biomarker identification. miRNAs in plasma are analyzed by methods that enrich for small RNAs and neglect pre-miRNAs or longer transcripts (9, 19, (21) (22) (23) . Additionally, almost all RNA-seq studies of plasma RNAs have used retroviral reverse transcriptases (RTs) to convert RNAs into cDNAs for sequencing on highthroughput DNA sequencing platforms. Retroviral RTs have inherently low fidelity and processivity and even those that are highly engineered have difficulty reverse transcribing through stable RNA secondary structures or post-transcriptional modifications, resulting in under-representation of 5'-RNA sequences and aborted reads that can be mistaken for RNA fragments (24) . Thus, we still have an incomplete understanding of the biology of plasma RNAs and how to optimize their identification as biomarkers for human diseases.

As a substitute for retroviral RTs, we have been developing RNA-seq methods using thermostable group II intron reverse transcriptases (TGIRTs) (20, 24, 25) . In addition to high fidelity, processivity, and strand-displacement activity, group II intron RTs have a proficient template-switching activity that enables efficient, seamless attachment of RNA-seq adapter sequences to target RNAs without RNA tailing or ligation (25, 26) . Further, unlike retroviral RTs, which tend to dissociate from RNAs at post-transcriptional modifications that affect base pairing, TGIRT enzymes pause at such modifications but eventually read through by characteristic patterns of mis-incorporation that can be used to identify the modification (20, 24, (27) (28) (29) (30) (31) . This combination of activities enables TGIRT enzymes to give relatively uniform 5'-and 3'-sequence coverage of mRNAs when initiating conventionally from an annealed oligo(dT) primer (25) and to give full-length, end-to-end reads of tRNAs and other structured small ncRNAs, when initiating by template switching to the 3' end of the RNA (20, 24, 25, 27) . In size selected RNA preparations, TGIRT-seq also profiles human miRNAs with bias equal to or less than alternative methods (25, 32) .

The comprehensive TGIRT-seq method used in this work to analyze human plasma RNAs employs TGIRT-template-switching for 3'-RNA-seq adapter addition followed by a singlestranded DNA ligation for 5' RNA-seq adapter addition (Fig. 1A) . In a validation study using rRNA-depleted chemically fragmented human reference RNAs with External RNA Control Consortium spike-ins, TGIRT-seq gave better quantitation of mRNAs and spike ins, more uniform 5' to 3' coverage of mRNA sequences, detected more splice junctions, particularly near the 5' ends of mRNAs, and had higher strand specificity when compared to benchmark TruSeq-v3 datasets (24) . This study also showed that TGIRT-seq enables simultaneous sequencing of chemically fragmented mRNAs together with tRNAs and other structured sncRNAs, which were poorly represented in the TruSeq datasets, even after chemical fragmentation (24) . A subsequent study of chemically fragmented human cellular RNAs with customized spike ins confirmed these findings and showed that TGIRT-seq enables simultaneous quantitative profiling of mRNA fragments and sncRNAs of >60 nt, but under-represents smaller RNAs due in part to differential loss during library clean-up to remove adapter dimers (33) . Recent TGIRT-seq profiling of structured RNAs in unfragmented cellular RNA preparations from human cells revealed previously unannotated sncRNAs, including novel snoRNAs, tRNA-like RNAs, and tRNA fragments (tRFs) (34) .

The ability of TGIRT-seq to simultaneously profile mRNAs and non-coding RNAs from small amounts of starting material without size selection is advantageous for the analysis of extracellular RNAs, which are present in low concentrations in human plasma or EVs secreted by cultured human cells. In a previous study, TGIRT-seq showed that human plasma from a healthy male individual contained predominantly full-length tRNAs and other structured small ncRNAs, which could not be seen in RNA-seq studies using retroviral RTs, together with RNA fragments derived from large numbers of protein-coding genes and lncRNAs, with higher proportions of intron and antisense RNAs than in cellular RNA datasets (20) . Similarly, TGIRT-seq showed that highly purified EVs and exosomes secreted by cultured human cells also contained predominantly full-length tRNAs, Y RNAs, and Vault RNAs, together with low concentrations of mRNAs, including 5' terminal oligopyrimidine (5' TOP) mRNAs (30) . The same TGIRT-seq workflow using template switching for facile 3'-RNA-seq adapter enabled single-stranded (ss) DNA-seq profiling of human plasma DNA, including analyses of nucleosome positioning and DNA methylation sites that inform tissue-of-origin (35) , as shown previously for other ssDNA-seq methods (36, 37) .

Since these earlier studies, TGIRT-seq methods have undergone improvements including:

(i) the use of modified RNA-seq adapters that substantially decrease adapter-dimer formation;

(ii) computational and biochemical methods for remediating residual end biases, enabling profiling of miRNAs with accuracy equivalent to the least biased current methods; and (iii) modified reaction conditions that increase the efficiency of RNA-seq adapter addition by TGIRT template switching, enabling the more efficient capture of target RNAs (26, 32) .

Here, we used an updated TGIRT-seq method incorporating these improvements together with RNA-seq adapters that add a unique molecular identifier (UMI) to deconvolute duplicate reads to comprehensively profile RNAs present in commercial human plasma from healthy individuals prepared by apheresis, a method used to obtain large volumes of plasma for clinical purposes, including convalescent plasma from recovered COVID-19 patients (38) . Additionally, we compared the fragmented mRNAs detected in plasma by TGIRT-seq with the polyadenylated mRNAs detected in the same plasma by ultra-low input SMART-Seq v4 and introduced the use of a peak-calling algorithm for analyzing TGIRT-seq datasets. We thus identified numerous mRNA fragments corresponding to annotated binding sites for ~100 different RNA-binding proteins. We also identified a wide variety of discrete structured RNAs and RNA fragments, including abundant repeat and transposable element RNAs, unannotated miRNA-sized RNAs, a family of putatively structured full-length, excised intron RNAs, some corresponding to mirtron pre-miRNAs and/or agotrons, and putatively structured intron RNA fragments, including a family corresponding to conserved structured segments of retrotransposed mRNAs that inserted within long introns.

TGIRT-seq of human plasma nucleic acids. The TGIRT-seq datasets in this study were obtained from commercial human plasma pooled from multiple healthy individuals and prepared by apheresis with EDTA as the anticoagulant (IPLA-N-K2E; Innovative Research). For each dataset, nucleic acids were extracted from 4 mL of plasma by using a QIAamp cfDNA/RNA kit (Qiagen) and treated with NaOH to obtain RNA-free plasma DNA (n=4) or with DNase I (n=12), exonuclease I (Exo I; n=1) or Exo I + DNase I (n =2) to obtain plasma RNA with minimal residual DNA (referred to collectively as DNase-treated plasma RNA; n = 15). Each preparation gave ~10 ng of nucleic acid, which yielded ~8 ng DNA or ~2 ng RNA, as judged by Bioanalyzer analysis (Fig. S1 ). Although the QIAamp cfDNA/RNA kit gave a lower yield of plasma RNA than other kits, it provided a uniform method for obtaining both RNA and DNA from the same plasma preparations.

To prepare TGIRT-seq libraries, we used the workflow outlined in Fig. 1A , which is based on the previously described TGIRT total RNA-seq method (20, 24) with the following improvements. First, the initial TGIRT-template switching reaction was done at a lower salt concentration (200 instead of 450 mM) that increases the efficiency of 3'-RNA-seq adapter addition (26, 35) . Second, we used a modified R2R adapter that substantially decreases adapterdimer formation (32) , thereby improving the representation of very small RNAs, such as miRNAs and short tRNA fragments (tRFs). Finally, we used a modified R1R adapter with a 6-nt UMI to deconvolute duplicate reads (35) . The libraries were sequenced on an Illumina NextSeq 500 instrument to obtain 10 to 27 million 75-nt paired-end reads, which were mapped to the hg19 human genome reference sequence. Mapping rates for the DNase-treated plasma RNA datasets (n=15) ranged from 82-96% (Supplementary File), with pairwise Pearson's correlation coefficients r = 0.67-1.00 (Fig. S2) . Although, we were concerned that the lower salt concentration used for template switching might increase the frequency of multiple sequential template switches, the percentages of soft-clipped and fusion/discordant read pairs, which include reads from multiple template switches, were relatively low (1.5-2.5% and 0.5-4.1%, respectively; Supplementary File).

The read-span distributions for the untreated and NaOH-treated plasma DNA samples showed two broad peaks, one at 130-180 nt, corresponding to plasma DNA protected in chromatosomes (nucleosomes plus linker histone) or trimmed mononucleosomes, and the other at 50-100 nt corresponding to nicked DNA strands that are protected by transcription factors or other bound proteins (Fig. 1B ), in agreement with previous studies (35) (36) (37) . Bioanalyzer traces using a High Sensitivity DNA kit showed that all of the dsDNA in the untreated plasma nucleic acid sample was sensitive to DNase I (Fig. S1A) . The read-span distributions for plasma RNA samples obtained after treatment with DNase I (n=12) or Exo I ± DNase I (n=3) showed a broad peak of 40-80 nt (Fig. 1B) , which corresponded to a similarly sized peak that was detected and confirmed to be NaOH-sensitive in Bioanalyzer traces obtained by using an RNA 6000 Pico Kit (Fig. S1B) .

Overview of RNA biotypes detected in human plasma. Fig. 1C -F show stacked bar graphs comparing the proportion of read pairs mapping to different genomic features in combined datasets obtained from plasma nucleic acids before and after different treatments. For NaOHtreated plasma DNA samples (n=4), the overall distribution of read pairs was similar to that simulated for random sampling of human genomic DNA (WGS-sim; simulated by ART). Most of these plasma DNA read pairs mapped to genomic repeats (46%; RepBase), unannotated sequences (26%), or protein-coding genes (21%; Fig. 1C ), with the read pairs mapping to protein-coding genes equally distributed between the sense and antisense strands (Fig. 1E ) and most of the bases in these read pairs located within introns (Fig. 1F ).

The plasma RNA samples after treatment with DNase I (n=12), Exo I ± DNase I (n=3), or DNase I + T4 polynucleotide kinase to remove 3' phosphates (n=3) gave read pairs mapping to 26,000-43,000 different genes (Fig. 1C, Fig. S2 and Supplementary File). For the DNase Itreated plasma RNA, high proportions of the read pairs mapped to cellular or mitochondrial (Mt) tRNAs (44% and 29%, respectively) followed by cellular and Mt rRNAs (16 and 5%, respectively), other sncRNAs (4.6%), repeats (RepBase; 0.8%), protein-coding gene transcripts (0.6%), other ncRNAs (0.1%), and miRNAs (0.006%; Fig. 1C ). Most of the rRNA reads corresponded to 18S and 28S rRNA (69%), with the remainder to 5S and 5.8S rRNAs (31%; Fig.   1C ). Treatment of plasma nucleic acids with T4 polynucleotide kinase to remove 3' phosphates increased the proportions of total reads mapping to 18S and 28S rRNAs (to 33%), repeats (to 9%), and protein-coding genes (to 4.1%; Fig. 1C ; lane 6). In agreement with previous findings (20) , the proportions of protein-coding gene reads mapping to the antisense strand (17-36%; Fig.   1E ) and bases within introns (63-73%; Fig. 1F ) were higher than those in TGIRT-seq datasets of cellular RNAs (10% and 24%, respectively (24)), suggesting that intron and antisense RNAs may be selectively exported or stable in plasma.

The most abundant sncRNAs in the DNase I-treated plasma RNA samples were tRNAs (36, 766 counts per million (CPM) for Val-CAC, the most abundant isodecoder) followed by Y RNAs (15,247 CPM for RNY4, the most abundant species), and 7SL RNA (158-1,504 CPM for the three 7SL RNA genes; Fig. 1D and Supplementary File). As found previously (20) , TGIRTseq showed that most of the cellular and Mt tRNAs detected in plasma were full-length mature tRNAs, which cannot be fully reverse transcribed by retroviral RTs (Fig. S3 and S4 , respectively). The plasma preparations also contained lower concentrations of specific cellular 5'-tRNA halves (0.02-68 CPM), 3'-tRNA halves (0.02-9,036 CPM), and shorter tRNA fragments (5'-tRFs, 0.02-39 CPM; 3'-tRFs, 3-640 CPM; and internal 2-tRFs, 0.02-7 CPM), with the proportion of reads corresponding to some 5'-tRNA halves and 5'-tRFs increasing after 3'phosphate removal (e.g., 5' half Gly-CCC, 2.9-fold increase; 5'-tRF Cys-GCA, 2-fold increase; see Supplementary File). Notably, TGIRT-seq identified 5' ends of 3'-tRFs that extend beyond the m 1 A modification site in the TψC loop, which is a strong stop for retroviral RTs and can be mistaken for the 5' end of 3'-tRFs ( Fig. S3) (39, 40) . Other structured sncRNAs in plasma, including Vault RNAs, snoRNAs, snRNAs, 7SK RNA, 7SL RNA, 5S rRNA and 5.8S rRNA, were also present largely as full-length mature transcripts (Fig. S5 ), in agreement with previous findings (20) .

The improved TGIRT-seq method detected 263 miRNAs annotated as high confidence miRNAs in miRbase v20 in the DNase I-treated plasma RNAs datasets (n=12; Supplementary File). As found previously using TGIRT-seq (20) , these were a mixture of pre-and mature miRNA, passenger strands, and longer transcripts, whose proportions varied for different miRNA species (Fig. S6) . The most abundant miRNA detected by TGIRT-seq, miR-223, was present at only 2.89 CPM (Fig. S6 ), considerably lower than the most abundant similarly sized Only a very low proportion of reads (<0.001%) mapped to piRNA loci ( Fig. 1D ; lanes 4-6), and none corresponded to a mature piRNA, even though under the lower salt conditions used for library construction TGIRT-III can template switch to and reverse transcribe RNAs with 2'-Omethylated 3' ends characteristic of piRNAs (26) .

To identify bacterial and viral RNAs in the commercial plasma preparations, we used the taxonomic classification program Kraken2 to analyze non-human reads in combined DNasetreated plasma RNA datasets (n=15; 3.3 million of ~110 million deduplicated reads including reads mapping to E. coli, which were filtered from the datasets and reintroduced for this analysis; see Methods). Kraken2 assigns reads to bacterial, viral, and archaeal genomes based on matches to unique k-mer sequences in a database of non-redundant bacterial, viral, and archaeal sequences compiled from Refseq. We thus identified 2.83 million RNA reads assigned to bacteria, with the most abundant phyla being Proteobacteria, Actinobacteria, and Firmicutes and the most abundant genera being Pseudomonas, Salmonella, Acinetobacter, Halomonas, and Bacillus, consistent with other studies (41, 42) . A much smaller number of non-human RNA reads (~52,000) were assigned to viruses, the most abundant being Parvoviridae and Mason-Pfizer monkey virus, a retrovirus (Fig. S7A ). The combined NaOH-treated plasma DNA datasets (n=4; 46 million deduplicated reads) contained only 12,000 non-human reads, much lower than 1% reported in another study using plasma prepared by a different method (43) . Most (93.3%) of the non-human DNA reads were assigned by Kraken2 to bacteria, mostly Pseudomonas spp., with only 64 reads assigned to viruses, mostly human β herpesvirus 6B ( Fig. S7B and Supplementary File). The commercial plasma preparations had been screened by the supplier for hepatitis B virus, hepatitis C virus, HIV-1, HIV-2, and syphilis, and we found no reads mapping to any of these pathogens.

TGIRT-seq versus ultra low input SMART-Seq v4. Our previous TGIRT-seq analysis of plasma from a healthy male individual prepared immediately from freshly drawn blood revealed that the mRNAs in plasma are largely fragmented, consistent with other studies of plasma RNAs using different RNA-seq methods with or without size selection (4, 13, 14, 18, 19) . In the commercial plasma from multiple healthy individuals prepared by apheresis, the proportion of the TGIRT-seq reads mapping to protein-coding genes was relatively low (0.6-4.1%; Fig. 1C , lanes 4-6), but nevertheless sufficient to the identify mRNAs originating from 15,000-19,000 different protein-coding genes, with the highest number of different mRNAs detected after 3' phosphate removal with T4 polynucleotide kinase ( Fig. S8A and B) . The identities and relative abundances of different protein-coding gene transcripts in plasma pooled for multiple healthy individuals were similar to those in the previous TGIRT analysis using plasma RNA from a single healthy male individual (20) (r = 0.61-0.80; Fig. S8C ).

To investigate the relationship between the mRNAs detected in plasma by TGIRT-seq and the polyadenylated mRNAs detected by other methods (e.g., (10, 17) ), we compared the proteincoding gene transcripts identified in the TGIRT-seq datasets for DNase I-treated plasma RNAs (n=12) with those identified in identically prepared RNA from the same plasma preparations by ultra low input SMART-Seq v4 (Takara Bio; n=4; Fig. 2 ). The latter is a highly sensitive method for profiling polyadenylated mRNAs in which an engineered retroviral RT (SMARTScribe) initiates cDNA synthesis from an anchored oligo(dT) primer with an appended PCR primerbinding site and then template switches from the 5' end of that mRNA to an acceptor oligonucleotide containing a second primer-binding site, enabling PCR amplification of the resulting cDNAs. In our experiments, the resulting dsDNAs were fragmented by Covaris sonication and used to prepare DNA-seq libraries by using an NEBNext Ultra II DNA Library Prep (New England Biolabs). The SMART-seq libraries were then sequenced on an Illumina NextSeq 500 to obtain 2.7-4.7 million 2 x 75 nt reads that mapped to >16,000 protein-coding genes with a transcript per million value (TPM) >0.1. Unlike TGIRT-seq, ultra-low input SMART-Seq does not preserve RNA strand information nor does it use an UMI to deconvolute duplicate reads. Bioanalyzer traces showed that the amplified dsDNAs obtained by SMART-Seq prior to Covaris sonication were substantially shorter than those from the cellular RNA control provided with the kit (peaks at 100-200 and 400-700 nt for plasma RNA compared to 700-7,000 nt for the cellular control RNA; Fig. 2A ), indicating that the poly(A)-containing mRNAs detected in plasma by SMART-Seq are enriched in 3' fragments.

A scatter plot comparing the relative abundance of transcripts originating from different genes showed that most of the polyadenylated mRNAs detected in DNase I-treated plasma RNA by ultra low input SMART-Seq were also detected by TGIRT-seq at similar TPM values when normalized for protein-coding gene reads (r=0.61), but with some, mostly lower abundance mRNAs undetected either by TGIRT-seq or SMART-Seq, and with SMART-seq unable to detected non-polyadenylated histone mRNAs, which are relatively abundant in plasma ( Fig. 2B and D). Similar correspondences were found for TGIRT-seq of plasma RNAs after 3'-phosphate removal (n=3; r=0.58) or after additional chemical fragmentation followed by 3'-phosphate removal (n=4; r=0.68; Fig. S8D and E). Heat maps comparing protein-coding gene transcripts detected by TGIRT-seq and SMART-Seq to primary tissue and platelet transcriptome data (44) (45) (46) indicated that the detected mRNA originated largely from hematopoietic tissues, including bone marrow, lymph nodes, and spleen, with contributions from erythrocytes, white blood cells, and platelets (Fig. S8F ), in agreement with other studies (10, 14) . The highly represented mRNAs detected by both TGIRT-seq and SMART-Seq included hemoglobins and other blood cell mRNAs (e.g., pro-platelet basic protein (PPBP) and S100 calcium-binding proteins S100A8 and S100A9), which are involved in regulating immune responses (47) (Fig. 2B and heat maps aminoacyl-tRNA synthetase (aaRS) mRNAs, which have potential intercellular signaling associated with different protein isoforms generated by alternative splicing (49) (Fig. 2B ). Both 5' TOP and aaRS mRNAs have been found in EVs secreted by cultured human cells (30, 49) , a potential vehicle by which these mRNAs might enter plasma.

Plots of normalized coverage over all detected mRNAs showed that SMART-Seq underrepresented 5'-RNA sequences and over-represented 3'-RNA sequences, as expected for 3'-mRNA fragments, while TGIRT-seq coverage, with or without 3'-phosphate removal or with chemical fragmentation followed by 3' phosphate removal, was relatively uniform across most of the length of the detected mRNAs (Fig. 2C ). Integrated Genomics Viewer (IGV) alignments for representative protein-coding genes confirmed the relatively uniform coverage of both polyadenylated and non-polyadenylated mRNAs by TGIRT-seq and showed that TGIRT-seq simultaneously detected an intron-encoded snoRNA in a ribosomal protein gene (RPL10; upper right that was invisible to SMART-Seq (Fig. 2D) .

Notably, although 3'-phosphate removal increased both the proportion of reads mapping to protein-coding genes and the number of mRNAs detected by TGIRT-seq (see above), neither the 5'-to 3'-coverage nor the abundance of most protein-coding transcripts detected in plasma by TGIRT-seq were strongly affected (Fig. 2C, Fig. S8A and B; r values = 0.93 and 0.84, respectively). This reflects that many of the additional protein-coding gene transcripts detected after 3'-phosphate removal had relatively low read counts and that the abundance of only a small proportion of protein-coding gene transcripts increased significantly after this treatment (0.41% by differential expression analysis using DESeq2; adjusted p-value < 0.01). Because 3' phosphates inhibit RNA-seq adapter addition by TGIRT template switching (25) , these findings are consistent with the suggestion from the previous TGIRT-seq analysis (20) that many mRNA fragments in plasma have 3' OH termini, as expected for RNA fragments generated by cellular ribonucleases that function in mRNA processing or turnover (50, 51) .

Use of peak-calling for analysis of plasma RNAs. As DNA fragments that persist in plasma were found to be packaged in nucleosomes or associated with transcription factors or other DNA binding proteins that afford protection from plasma nucleases (35) (36) (37) , we wondered whether many of the mRNA fragments in plasma might be similarly protected by bound proteins and whether such protein-protected RNA fragments could be detected by peak calling, as done for ChIP-seq (52) . To test this idea, we combined the TGIRT-seq datasets for all DNase-treated plasma samples (n=15) and removed reads from human blacklist regions, which are known to produce artifactual peaks (53) , as well as reads corresponding to rRNAs, Mt RNAs, and annotated sncRNAs, which were analyzed above. We were left with ~1.6 million deduplicated read pairs (~43% of the total) that mapped to the hg19 human genome reference sequence to use as input for the ChIP-seq narrow peak caller, MACS2 (52) . As a base line control, we used the ~38 million deduplicated read pairs that mapped to human genome from the NaOH-treated plasma DNA samples (n=4).

After peak calling using a read coverage cutoff of ≥5, a false discovery rate cutoff of 0.05 (q-value assigned by MACS2), and a requirement that the peak be detected in at least 5 of the 15 samples to avoid batch effects, we were left with 1,036 peaks that were enriched in the DNasetreated plasma RNA over the base line control. After further filtering to remove reads with low mapping quality (MAPQ <30) or ≥ 5 mismatches from the mapped locus and examining peaks containing such reads individually to be sure that nothing significant was discarded, we were left with 950 high confidence peaks ranging in size from 59 to 1,207 nt (Supplementary File). The percentage and number of the peaks mapping to different annotated features are shown by the pie charts in Fig. 3 . Among the detected peaks, 922 had one or more annotations for a gene (Ensemble), RBP-binding site (ENCODE K-562 and HepG2 cells eCLIP datasets (54)), or repeat sequence (RepBase) on the sense strand; 25 had one or more annotations on the antisense strand;

and three had no overlap with any of these annotations (Fig. 3A) .

To understand the origins of these peaks, we first extracted those with sense-strand annotations and plotted the proportions corresponding to different annotated genomic features Because some peaks corresponded to annotated binding sites for more than one RBP, for purposes of this count, we assigned the peak to the best matched RBP (or RBPs in case of ties) scored by the number of overlapping bases. EFTUD2, an RNA splicing factor (Fig. 4D) . Farther down the list were PRPF8, a spliceosomal protein that provides a structural scaffold for the assembly of snRNAs at splice sites (Fig. 4E ); SND1, a transcriptional co-activator; and UCHL5, an RBP that modulates mRNA expression ( Fig. 4F ). Also prevalent in plasma were RNA fragments containing binding sites for DExH-box RNA helicase DHX30; GRWD1, a histone-binding protein that may also play a role in ribosome assembly; ZNF622m, a zinc-finger protein; BUD13, a spliceosome component; PPIG, a protein involved in protein-folding; and SERBP1, a mRNA-binding-protein. Although these findings do not prove that the identified RBP was associated with the RNA fragments identified in plasma, collectively they suggest that many of the mRNA fragments that persist in plasma are protected from plasma nucleases by bound RBPs.

Peak calling against a human genome reference sequence might miss RBP-binding sites that overlap splice junctions. To address this possibility, we mapped the reads in the combined DNase-treated datasets (n=15) to a human transcriptome reference sequence (Ensemble human cDNA references) and obtained ~0.26 million deduplicated read pairs that were used as input for MACS2 without a control dataset (52) . Using the same peak-calling parameters as above, we identified 865 peaks that were enriched over the base line. After further filtering to remove annotated sncRNAs and low confidence peaks with RNAs in these peaks are described below along with those of additional intron peaks that were not annotated as an RBP-binding site.

Peaks mapping to exons and pseudogenes. Two hundred and thirteen peaks that mapped to exons or pseudogenes in the human genome reference sequence did not correspond to ENCODE eCLIP-annotated RBP-binding sites or genomic repeats, and an additional 33 such peaks were identified by mapping to the human transcriptome (Figs. S10-S13 and Supplementary File).

These peaks ranged in size from 52-1,207 nt and could correspond to RNAs containing RBPbinding sites that are not annotated in the two eCLIP datasets, unannotated structured sncRNAs, or simply reflect uneven sequence coverage across mRNA sequences. Twenty five of these 246 peaks were classified by Infernal/Rfam analysis (57) into four types of known RNA structures: 21 histone 3'-UTR stem-loops required for 3'-end processing of histone mRNAs; an iron response element (IRE I) in a ferritin light chain (FTL) exon; a selenocysteine insertion sequence 1, which directs the translation of UGA as SelCys; and two potential pre-miRNA stem-loop structures, with the closest match for both being mouse miR-692, which has no annotated human homolog ( Fig. S10 and Supplementary File).

We inspected IGV plots for all of the remaining 221 exon and pseudogene peaks and found that most were comprised of RNA fragments extending across protein-coding gene exons ( 

Almost a quarter of the RNA peaks identified by peak calling (213 peaks; 23.1%) corresponded to RepBase-annotated sequence repeats, including 154 comprised of or containing short tandem repeats (also referred to as simple repeats or microsatellite sequences) and 59 to transposable element RNAs (Fig 3B and Supplementary File). TGIRT enzymes are advantageous for the analysis of tandem repeat RNAs because of their ability to reverse transcribe through stable higher-order RNA structures formed by interactions between the repeat units, enabling them to more completely reverse transcribe and better quantitate these RNAs than can retroviral RTs (58) . The peaks for the short tandem repeat RNAs in human plasma had a surprisingly narrow length distribution (sharp peak ~80 nt) with >50% of the short tandem repeat RNAs having complete or nearly complete read coverage across the annotated repeat region ( Fig. S14A and B). Of the top 15 RepBase-annotated repeat types with the highest peak count, 11 were low complexity sequences, such as polypurine or polypyrimidine tracts or simple repeat RNAs (e.g., (CCG)n and (TTC)n), and four were transposable element RNAs (AluSx, AluJb, FLAM_C and LINE-1 element HAL1-3A_ME) ( We noticed that some RepBase simple repeat sequences (e.g., (CATTC)n and (TTAGGG)n) were called by MACS2 at some loci with roughly equal numbers of reads on both strands, possibly reflecting mismapping of DNA or RNA reads. To address this issue, we used an Empirical Bayes method (Methods and Fig. S14C ) to identify the top 15 repeat sequences with the largest differences in the percentage of sense (+) orientation reads in DNase-treated plasma RNA compared to NaOH-treated plasma DNA and found that all corresponded to tandem repeats of short (3-6 nt) sequences, including centromeric and telomeric repeats (Fig. 3E, Fig. 5C and D, and Fig. S14C and D). We also detected potential HSATII RNAs, which had been identified previously as a repeat RNA present in plasma (12) , but with similar numbers of (+) and (-) orientations reads (2,327 and 2,414, respectively) and not called as peaks by MACS2.

Differential expression analysis for transposable element RNAs found that those with highest enrichments in plasma RNA compared to plasma DNA included endogenous retroviral, LINE-1, and Alu RNAs (Fig. S14E ). In many but not all cases, the transposable element peaks corresponded to RNA fragments with peak lengths in the same narrow range found for simple repeats ( Fig. 5G-J, Fig. S14A ). The detection of short tandem repeat and transposable element RNAs as relatively discrete peaks in plasma could reflect RNA structural features or bound proteins that protect from plasma nucleases.

Peak calling identifies a family of putatively structured, full-length excised intron RNAs. As indicated above, 29 peaks identified as containing RBP-binding sites were full-length excised intron RNAs. An additional 14 such peaks were found among those mapping to long RNAs but not containing an annotated RBP-binding site, and another seven (denoted by IDs beginning with M) were not called as peaks but were identified as containing annotated miRNAs sequences. As discussed further below, seven of these 50 intron peaks corresponded to the same intron found in the PKD1 gene and six of its pseudogenes. Counting these seven peaks as a single intron, we AG with no indication of impediments that might be indicative of a lariat RNA. In several cases, the reads corresponding to full-length excised intron RNAs had 1 or 2 non-templated U or A residues at their 3' end, reminiscent of non-templated 3' A-tails found for yeast linear introns that accumulate in stationary phase cells (59) .

Thirteen of the full-length excised intron RNAs that we detected in plasma, including the most abundant (ID#468 corresponding to DOCK6 intron 25, 389 deduplicated reads, 7 CPM), corresponded to annotated agotrons, intron RNAs that were identified as binding to argonaute-2 protein (AGO2) in CLIP-seq datasets and shown to repress mRNA translation in reporter assays (60, 61) (Fig. 3F, Fig. 6A and Supplementary File). Six other possible agotrons (peak ID#s 2, 404, 416, 522, 844, and 846) were identified by intersecting the full-length excised intron RNAs with AGO1-4 PAR CLIP datasets (62) . Agotrons had been hypothesized to be full-length linear intron RNAs based on Northern hybridization experiments and CLIPseq 5'-end sequences (61), but to our knowledge, this has not been confirmed previously by RNA-seq, possibly because the retroviral RTs used in other RNA-seq methods are unable to fully reverse transcribe these structured RNAs.

Ten of the excised intron RNAs that we detected in plasma corresponded to annotated mirtrons, miRNAs that are processed by DICER from debranched structured intron RNAs (examples shown in Fig. 6B and C) (63) (64) (65) . Seven of these ten mirtron pre-miRNAs were also annotated as agotrons (including four of the top 15; green with asterisk in Fig Fig. S15A) , part of the predicted stem-loop structure for the called peak corresponded to a separate 19-nt RNA (red), which comprised a major component of the peak. In one of these peaks, the 19-nt RNA (red) was accompanied by a complementary 22-nt antisense RNA (blue; Fig. 8B ), and in the other two , the 19-nt RNA (red) was part of a longer 48-nt tandem repeat unit (green) within the predicted stem-loop structure ( Fig. 8C and Fig.   S15A ). Two other peaks, one corresponding to an annotated binding for the dsRNA-binding protein ILF3 (peak ID#677), contained complementary segments of long (46 and 65 bp) inverted repeats ( Fig. 8D and E) . Another of the putatively structured intron RNA fragment peaks (ID#731) was identified by SnoGPS (66) as being able to form secondary structures resembling an H/ACA-box snoRNA, but these did not correspond to the most stable secondary structure predicted by RNAfold for this peak (Fig. S15C ).

The six intron peaks that could not be folded into stable secondary structures had other features that might contribute to nuclease resistance in plasma. Three of these peaks consisted of AG-rich sequences or tandem repeats ( Fig. 8F and G, and ID#737), including one with tandem AGAA repeats identified as an annotated binding site for TRA2A, a protein that controls alternative splicing (Fig. 8G ) (67) . Two others contained one arm of a long inverted repeat sequence, whose complementary arm lies outside of the called peak (ID#s 756 and 761), and the remaining peak was a highly AU-rich RNA (74% AU; ID# 964; Supplementary File).

Finally, we found that the nine additional peaks (71-295 nt) mapping to internal segments of intron RNAs were part of retrotransposed mRNAs sequences (432-1,745 nt) that had integrated into seven different long introns (8 to 138 kb; Fig. 8H and I, Fig. S16 ). These retrotransposed mRNA sequences were identified by a BLAT search for >95% identity to the called peak in the human genome reference sequence and originated from six different highly expressed mRNAs (ribosomal protein mRNAs RPS3, RPL18A, and RPL41; translation elongation factor mRNAs EEF1A1 and EEF1G; and β-actin mRNA ( Fig. 8H and I, Fig. S16 ). All of these retrotranposed mRNA sequences had a short poly(A) tail (9-22 nt) and were flanked by short direct repeats (target site duplications (TSDs; Fig. S16 ), hallmarks of LINE-1 RT-mediated retrotransposition (68) . In all cases, the RNA peak identified in plasma corresponded to a smaller segment of the retrotransposed mRNA sequence that could be folded by RNAfold into a stable RNA secondary structure (ΔG ≤ -16.4 kcal/mole; Fig. 8H and I, Fig. S15D ). Comparison of genomic sequences showed the retrotransposition events that inserted the mRNA sequences within the introns were relatively recent, with six occurring in primates and two in placental mammals (Fig. S16 ). In all of these cases (and three additional cases of retrotransposed mRNA sequences described below), the RNA peak identified in plasma corresponded to a segment of the retrotransposed mRNA whose genomic sequence was identical to that in the gene encoding the functional mRNA (verified by searching the peak sequence against the hg19 human reference genome using BLAT and manually checking the sequences in the human genome browser and IGV), making it impossible to determine the origin of the RNA detected in plasma. Analysis of two introns with the largest number of homologous genomic sequences across primates and placental mammals showed that both the retrotransposed mRNA sequence and called peak region were more conserved than the flanking intron sequences, consistent with functional importance (Fig. S17) .

Peaks corresponding to unannotated genomic loci. Only three of the peaks identified by peak calling mapped to unannotated genomic loci (Fig. 3A) . One corresponded to a different putatively structured 3' segment of a retrotransposed ACTB mRNA, with a short poly(A) flanked by two different Alu elements with a third Alu element inserted within the retrotranspsosed mRNA sequence (Fig. S18A) . The second was a 72-nt peak that included multiple copies of a 17nt RNA with discrete 5'-and 3'-ends (Fig. S18B) , and the third was a 158-nt peak consisting of TCCAT(C) 4-6 GTG repeats (Fig. S18C) .

Peaks corresponding to antisense RNAs. Finally, 25 peaks were identified as antisense transcripts of annotated genomic features (Fig. 3A) . Eight of these peaks mapped to introns and the remaining 17 peaks mapped to mRNA exons or pseudogenes. Among the eight antisense peaks mapping to introns, four corresponded to or contained multiple reads for four different discrete short RNAs (<20 nt) complementary to a sequence within the intron (Fig. 9A-D and Supplementary File), with one corresponding to a 13-nt segment of miR-4497, a non-highconfidence predicted miRNA encoded at another locus (Fig. 9A) and another corresponding to a 17-nt segment of a 36-nt tandem repeat (Fig. 9C) . Two other antisense peaks mapping to introns were putatively structured segments (ΔG ≤ -22 kcal/mol) of retrotransposed TMSB4X and FTH1 mRNAs that had been inserted within introns in the antisense orientation relative to the host gene ( Fig. S19A and B) , and the remaining two peaks (ID#973 and 539) contained RNA fragments with heterogenous 5' or 3' ends, with one (ID#539) having the potential to fold into a stable secondary structure (ΔG = -26.7 kcal/mol; Fig. S19C and D) .

Among the 17 other antisense peaks mapping to mRNAs or pseudogenes, 12 mapped to hemoglobin genes (HBA2, HBA1, HBQ1, HBB, etc. ), four to other highly represented genes (RPS29, TMSB4X, INPP5E, and FTL), and one to a pseudogene (Fig. 9E-H and Fig. S20A-C) .

The peak mapping to INPP5E included a discrete short (17-nt) RNA (Fig. 9E ). While many of these peaks may be bona fide antisense transcripts, some included antisense reads that extended across two or more spliced exons (horizontal orange boxes) and/or were mirrored by a partially overlapping antisense DNA peak (red boxes extending vertically across the coverage tracks and read alignments; Fig. 9F -H and Fig. S20 ), the latter suggesting that they contain DNA fragments that are partially protected from DNase I-digestion in an RNA/DNA duplex. Such RNA/DNA duplexes could be remnants of R-loops formed during transcription (69) or cDNAs generated by reverse transcription of a spliced mRNA by an endogenous cellular RT. None of the Reads 1 of these possible antisense DNAs began with an R2 adapter sequence as would be expected for recopying of an initial cDNA by TGIRT-III.

Here, we used TGIRT-seq combined with peak calling to comprehensively profile RNAs present in apheresis-prepared human plasma pooled from multiple healthy individuals.

Extending previous TGIRT-seq analysis of plasma RNAs from a single male individual (20) , we found that human plasma contains largely fragmented mRNAs originating from >19,000 different protein-coding genes, together with abundant full-length mature tRNAs and other structured sncRNAs, as well as specific tRNA fragments, and mature and pre-miRNA, in total representing RNAs originating from >40,000 different human genes. By using peak calling, we found that many of the mRNAs, repeat RNA, and intron RNA fragments that persist in plasma correspond to annotated RBP-binding sites and/or have stable predicted RNA secondary structures or other structural features that may afford protection from plasma nucleases.

Additionally, we identified a family of short putatively structured full-length excised introns RNAs, subsets of which correspond to agotrons and/or mirtrons, as well as putatively structured intron RNA fragments, including a family mapping to conserved structured segments of retrotransposed mRNAs that inserted within long introns. Although not originally a focus of our study, we found that commercial plasma obtained by apheresis from healthy individuals had no detectable reads for pre-screened pathogens (hepatitis B virus, hepatitis C virus, HIV-1, HIV-2, and syphilis), but did contain low levels of reads from bacteria and other RNA and DNA viruses ( Fig. 8) , with possible implications for the use of such plasma for clinical purposes (38) . The ability of TGIRT-seq to simultaneously profile a wide variety of RNA biotypes in human plasma, including structured RNAs that are intractable to retroviral RTs, may be advantageous for identifying optimal combinations of coding and non-coding RNA biomarkers for human diseases.

By using TGIRT-seq with RNA-seq adapters containing a UMI to deconvolute duplicate reads, we confirmed that full-length mature tRNAs are the most abundant class of sncRNAs in human plasma, followed by Y RNAs and 7SL RNA, with tRNA halves, short tRFs, and mature and pre-miRNAs and passenger strands present in lower concentrations. The full-length, mature tRNAs and other sncRNAs present in plasma may be protected from plasma nucleases by stable RNA secondary or tertiary structure, bound proteins, or encapsulation in EVs, the latter suggested by TGIRT-seq analysis of RNAs present in highly purified EVs and exosomes secreted by cultured human cells, which likewise identified full-length mature tRNAs and Y RNAs as the most abundant RNA species (30) . Although the updated TGIRT-seq methods used in this work improved the detection of miRNAs and 5'-and 3'-tRFs, further analysis of these RNAs by RT-qPCR and Fireplex assays indicated that some miRNAs remained substantially underrepresented relative to RNAs >60 nt (manuscript in preparation). Thus, while comprehensive TGIRT-seq of heterogeneously sized RNA preparations can accurately determine the relative abundance of RNAs >60 nt (24, 33) , determining the quantitative relationship between very short RNAs and RNAs >60 nt requires orthogonal methods. Surprisingly, despite the size bias against very short RNAs, we identified 15 novel discrete short RNAs (13-25 nt; Fig.   8 , Fig. 9, Fig. S13, Fig. S15, Fig. S18 and Supplementary File). Only one of these short RNAs (non-high-confidence miR-4497; Fig. 9A ) could be found in miRBase, piRNAdb, tRFdb or cross-linked regions in AGO1-4 PAR-CLIP datasets (62) , raising possibility that they represent other classes of short regulatory RNAs, whose protein-partners and functions, if any, remain to be determined.

The proportion of reads corresponding to mRNAs in plasma prepared by apheresis was relatively low (<5% of mapped reads), but nevertheless sufficient to identify mRNA fragments originating from >19,000 different protein-coding genes, with uniform coverage across mRNA exons for highly represented genes (Fig. 2) . The number of different mRNAs detected in plasma by TGIRT-seq was comparable to that detected in parallel assays of the same plasma by ultra low input SMART-Seq v4, which uses oligo(dT) priming, and higher than that in other studies using different methods that sequence mRNA fragments with or without size selection (4, 13, 14, 41) . Notably, the size distribution of cDNAs generated by ultra low input SMART-seq suggested that a substantial proportion of the mRNA fragments present in plasma have a length >200 nt (36%; Fig. 2A ), which would be lost by incorporating a size-selection step. Although the presence of other more abundant RNA biotypes is a limitation for mRNA detection by TGIRTseq, the proportion of reads corresponding to mRNA could be increased by depleting rRNAs and other abundant non-coding RNAs, by using a maximally efficient RNA extraction kit, and possibly by using plasma prepared from freshly drawn blood to minimize degradation of more labile mRNA fragments, as was done in the previous TGIRT-seq analysis of human plasma RNAs, which gave substantially higher proportions of mRNA reads (20) .

The mRNA fragments present in plasma are likely a mixture of transient RNA degradation intermediates and more persistent RNA fragments that are protected from plasma nucleases by bound proteins, RNA structure, encapsulation in EVs, or a combination of these factors. Both the present and previous TGIRT-seq analyses indicated that many of the mRNA present in plasma have 3' OHs, which are required for efficient 3' RNA-seq adapter addition by TGIRT-templateswitching, particularly under the high-salt reaction conditions used in the previous study (20) . As most cellular enzymes that function in RNA processing or turnover leave 3' OH groups while extracellular RNases, such as RNase 1, leave 3' phosphates (50, 51, 70) , these findings suggest that a high proportion of the mRNA fragments detected in plasma were generated by intracellular RNases that function in mRNA turnover rather than degradation of intact mRNAs that were released into plasma.

To identify those mRNA fragments in plasma that might be associated with bound proteins, we used a peak-calling algorithm and intersected the called RNA peaks with Encode eCLIP annotated RBP-binding site sequences. We thus identified numerous discrete mRNA and lncRNA fragments that are annotated binding sites for ~100 different RBPs, most of which were previously identified as part of the plasma proteome. Although these findings do not directly demonstrate the association of the annotated RBP with the protected RNA fragment, collectively they suggest that many of the mRNA fragments in plasma are protected from plasma nucleases by bound proteins, similar to the inferred protection of plasma DNA fragments by nucleosomes packaging or by the binding of transcription factors or other proteins (35) (36) (37) . The number of protein-protected mRNA fragments that we identified in plasma is likely underestimated as it was based on only two Encode eCLIP datasets and used plasma prepared using the anticoagulant EDTA, which destabilizes some RNPs (71) .

In addition to bound proteins, many of the mRNA fragments that we detected in plasma with or without annotated RBP-binding site sequences were predicted by RNAfold to form stable RNA secondary structures that might contribute to their protection from plasma nucleases. The predicted RNA secondary structures of some mRNA fragments could be classified by Infernal/Rfam as corresponding to known structural elements (e.g., histone 3' stem-loop structures, iron response elements, selenocysteine insertion sequence, and pre-miRNA stem-loop structure; Fig. S10 ), but in most cases, the predicted secondary structure could not be classified by Infernal/Rfam, and its functional significance, if any, remains to be evaluated. Some peaks corresponding to structured mRNA fragments had discrete 5' and 3' ends and could be unannotated sncRNAs (Fig. S13A-D) , while others consisted of or included discrete miRNAsized RNA, possibly unannotated short regulatory RNAs (Fig. S13E-G) . Additionally, some of the mRNA species that we identified as being abundant in plasma, e.g., 5' TOP and aaRSs, were shown previously to be present in EVs and exosomes secreted by HEK-293T cells, a potential vehicle by which these mRNAs might enter plasma (30) .

A surprisingly high proportion of the discrete RNAs identified by peak calling (23.1%) corresponded to RepBase-annotated short tandem repeat and transposable element RNAs (Fig.   3B ). TGIRT enzymes have been found to more completely reverse transcribe and better quantitate short tandem repeat RNAs than do retroviral RTs and can also give full-length reads of Alu and other structured SINE elements RNAs (27, 58) . In all, TGIRT-seq detected discrete RNAs containing 37 different types of short tandem repeats, including telomeric and centromeric repeats. The short tandem repeat RNAs had a narrow size distribution (peak at ~80 nt) with most encompassing a high proportion of the genomic repeat unit. Many but not all of the transposable element and endogenous retroviral RNAs were RNA fragments of a size similar to that of the simple repeat RNAs. The ability of TGIRT-seq to profile repeat RNAs in plasma may be useful for liquid biopsy of RNA repeat expansion diseases, such as myotonic dystrophy and some forms of amyotropic lateral sclerosis.

TGIRT-seq peak calling also identified novel intron RNAs, including a family of short, full-length excised intron RNAs. These full-length excised introns ranged in length from 73-130 nt, are likely linear RNAs, had RNAfold-predicted secondary structures with ΔGs < -18 kcal/mole, and in many cases, corresponded to annotated binding sites for spliceosomal proteins and/or AGO1-4. Twenty three or the 44 full-length excised intron RNAs that we detected in plasma corresponded either to annotated or potential mirtrons, intron pre-miRNAs that are processed by DICER to produce mature-miRNAs (63) (64) (65) and/or to annotated or potential agotrons, structured excised intron RNAs that bind AGO proteins and function directly to repress mRNA translation (60, 61) . Seven of the full-length excised intron RNAs that we detected in plasma were annotated as both an agotron and mirtron, blurring the distinction between these two classes of RNAs. Notably, although we detected mirtron pre-miRNAs in plasma, we did not detect the corresponding mature miRNAs, raising the possibility that mirtron pre-miRNAs are selectively exported or have higher stability in plasma. While extracellular miRNAs have been suggested to function in intercellular communication (72, 73) , mirtron pre-miRNAs and agotrons may be as well or better suited for this role, with the mirtron pre-miRNAs entering cells as precursor RNAs that could be processed by DICER into mature miRNAs, and agotrons bound to AGO proteins entering cells as an RNP that could function directly in regulating gene expression in a miRNA-like manner.

Surprisingly, almost half (21 of 44) of the structured full-length excised introns that we detected in plasma did not correspond to identified agotrons or mirtrons. These introns could be mirtrons or agotrons that have yet to be annotated, or they could have some other functions that remains to be determined. One of the most abundant full-length excised intron RNAs that we detected in plasma, PKD1 intron 29 (ID#333; Fig. 6 ), was highly conserved in both the parent gene and six distant pseudogenes (1 mismatch at a known SNP position), possibly reflecting selection for a sequence-dependent function of this intron. Nevertheless, only one of the fulllength excised intron RNAs that we detected in plasma, a putative mirtron encoding miR-1225, was conserved in sequence across vertebrates and mammals, and only six (not including PKD1 intron 29) were conserved in sequence across primates (PhasCons score ≥0.5; Fig. 7) , suggesting that any sequence-dependent function of these intron RNAs would have been recently acquired.

In addition to full-length excised introns, peak calling identified a variety of intron RNA fragments, with more than half (16 of 31) corresponding to annotated RBP-binding sites and most having stable predicted RNA secondary structures or other structural features (tandem repeats, inverted repeats, or highly AG-or AU-rich sequences) that may afford protection from plasma nucleases. The RBPs potentially associated with these intron RNA fragments were a largely different set than those associated with the full-length excised intron RNAs, and in some cases, could be related to a specific RNA splicing function or structural feature of the intron.

Thus, an intron RNA peak containing complementary segments of a long-inverted repeat corresponded to an annotated binding for the dsRNA-binding protein ILF3, and an intron RNA peak containing tandem AGAA repeats corresponded to annotated binding site for the alternative splicing factor TRA2A, which recognizes this motif (Fig. 8) . Some of the intron RNA peaks that could be folded into stem-loop structures also contained unannotated discrete miRNA-sized RNAs (13-25 nt) that may have been processed out of the stem-loop structure. The latter finding raises the possibility that whatever their origin or original function, if any, structured excised intron RNAs in human cells may be accessible to DICER or other RNase III-like enzymes that cleave double-strand RNAs into discrete small fragments to generate a repository of non-coding short RNAs that could potentially evolve miRNA-like functions.

By analyzing RNA peaks that mapped to both introns and mRNAs, we also identified multiple human introns containing retrotransposed segments of highly expressed mRNAs, conserved, structured portions of which persist in plasma (Fig. 8, Fig. S15, Fig. S16, and Fig.   S17 ). In all cases, the retrotransposition event that inserted the mRNA sequence within the intron appears to have been mediated by a LINE 1 RT and was relatively recent, with most occurring in primates and a few occurring earlier in placental mammals. In all cases, the RNA peak detected in plasma matched a putatively structured portion of the mRNA sequence that was highly conserved within the intron, raising the possibility of functional importance. Additionally, for several highly expressed mRNAs, including three that were progenitors of processed pseudogenes in the human genome (TMSBX4, RPS29, and FTL), we detected antisense reads that could be cDNAs generated by an endogenous RT, possibly precursors of continuing retrotransposition events in human cells ( Fig. 9 and Fig. S20 ). Finally, we note that in addition to their biological and evolutionary interest, full-length excised intron RNAs and intron RNA fragments, like those we detected in plasma, may be uniquely well suited to serve as stable RNA biomarkers, whose expression is linked to that of multiple protein-coding genes. sequence (R1R) and the 6-nt UMI at its 5' end was then ligated to the 3' end of the cDNAs by using Thermostable 5′ AppDNA/RNA Ligase (New England Biolabs) as described (35) . After clean up with a MinElute Reaction Cleanup Kit (Qiagen), ligation products were amplified by PCR using a KAPA Library Amplification Kit (KAPA Biosystems KK2610) with 500 nM of Illumina multiplex primer (a 5′ primer that adds a P5 capture site) and 500 nM of index primer (a 3′ primer that adds an index and a P7 capture site). PCR was done with initial denaturation at 98 °C for 30 sec followed by 12 cycles of 98 °C for 45 sec, 60 °C for 15 sec, and 72 °C for 30 sec, and an additional final incubation at 72 °C for 5 min. The PCR products were purified by using 1.3X Agencourt AMPure XP beads (2 or 3 cycles; Beckman Coulter) and evaluated by using an Sequence data processing. Fastq files were processed by using the TGIRT-map pipeline (74) with modifications for preserving UMI information and more stringent read trimming and filtering, including trimming of partial adapters and trimming each adapter sequence twice for each read. 6-nt UMIs were clipped and appended to the read ID, and only read pairs with average UMI Phred quality score <20 were retained. Adapters, sequencing artifacts (e.g., long homopolymers), and TGIRT-seq byproducts (e.g., primer dimers) were removed with Atropos 

Reads that mapped to protein-coding genes were extracted and quantified by Kallisto (79) using a transcript reference generated by gffread (80) from the hg19 human genome reference sequence and GENCODE transcript annotations release 28. Gene expression profiles for human tissues were downloaded from the Human Protein Atlas (46) . Six platelet RNA-seq datasets from healthy individuals were obtained from NCBI Sequence Read Archive (SRR5907423-SRR5907428) (44) and TPM values were quantified by Kallisto using the hg19 reference sequence. The percentage of bases that mapped to different regions of protein-coding genes was computed by CollectRnaSeqMetrics from Picard (Broad Institute) using a genome BAM file generated from Kallisto. Gene body coverage was computed using RSeQC (81) .

Unmapped reads from datasets for DNase-treated plasma RNA (n=15) or NaOH-treated plasma DNA (n=4) were combined with reads that mapped to E. coli from read-mapping Pass 1 (see above) and analyzed by Kraken2 (82) and additional annotated sncRNAs, including snoRNAs, snRNAs, non-high confidence miRNAs and miscellaneous RNAs (Ensembl; https://useast.ensembl.org/index.html) and then converted to a BED file (84) and deduplicated by UMI and genome coordinates. The counts for each read were readjusted to overcome potential UMI saturation for highly-expressed genes by implementing the algorithm described in Fu et al. (85) , using sequencing tools (https://github.com/wckdouglas/sequencing_tools). BED files from libraries with the same treatments were combined. Reads mapping to human genome blacklist regions (53) were removed. Because we were using human plasma pooled from multiple individuals, misalignments of Mt tRNAs from different Mt DNA haplogroups may produce spurious enrichments in nuclear mitochondrial DNA segments (NUMTs). To remove these false positive enrichments, we also filtered out read pairs that had at least one read that could align to the hg19 mitochondrial genome. The BED file for the combined DNase-treated datasets (n=15) was then split into separate BED files containing only forward strand or reverse strand fragments. BED files from each strand were used as input the for MACS2 (52) 

.was computed for each overlapping record between a peak and a genomic feature. An overlap score of 1 (highest possible overlap score) indicates the peak can be explained by the full-length genomic feature. A best feature annotation for each peak was selected by the highest overlap score. In cases of identical overlap scores, a feature was selected with a priority to RBP-binding sites, repeat regions (RepBase), and NaOH-treated plasma DNA (n=4) datasets. To compensate for repeat sequences with low read counts, which could give misleading strand-specificity, we employed an Empirical Bayes method by using repeat sequences with high total counts (>50 deduplicated reads) regardless of strand from the combined plasma DNA datasets to construct a prior beta distribution of strand specificities for all repeat sequences (Fig. S14C ). This fitted beta distribution represents the probability distribution of detecting 0-100% (+) orientation reads for each repeat sequence in a pure DNA sample. Using the fitted beta distribution hyperparameters and the count data from either the plasma DNA or plasma RNA datasets, we then computed a posterior beta-binomial distributions of the percentages of (+) orientation reads for each repeat sequences, which represents how likely and to what degree the strand specificities were biased towards the (+) or (-) strand orientation. We then compared the posterior distributions for each unique repeat sequences between plasma RNA and DNA to obtain a distribution of differences in the percentage of (+) strand reads in plasma RNA versus plasma DNA (Δ %(+) in RNA vs DNA).

Summary statistics were then extracted from this final distribution, and repeat sequences with Bayes factor >3 for at least 10% more (+) strand fragments in the plasma RNA compared to the plasma DNA samples were identified as significantly enriched in (+)-strand fragments in plasma RNA. The Bayesian statistical testing framework was implemented using pymc3 (89) .

Code availability: All scripts used for data processing are deposited in GitHub:

https://github.com/wckdouglas/cfNA Annotated protein-binding sites are shown as color-coded lines below the coverage plots. For cases in which the peak overlaps multiple RBP-binding sites, the one whose annotated binding site had the largest number of bases overlapping the peak is underlined. Colors other than pink or purple in the read alignments indicate bases that do not match the reference sequence (red, green, blue, and brown indicate thymidine, adenosine, cytidine, and guanosine, respectively). The peak called by MACS2 is delineated by a bracketed line with the peak ID, length, and PhastCons score for 46 vertebrates including humans indicated below (see phylogenetic tree in Fig. 7) . The most stable predicted secondary structure for the peak and its minimum free energy (MFE; ΔG) computed by RNAfold (88) Reads shown in the alignment were down sampled to 100 reads when necessary for display.

NTA/black boxes, non-templated nucleotides added by TGIRT-III. When necessary, reads shown in alignment track were down sampled to maximum 100 for display. 

ID#975, 80 nt PhastCons: 0.00 5'-UCUUCCUCACCCUUACCAAGGA-3' 

Cell-free nucleic acids as biomarkers in cancer patients

Comparative RNA-sequencing analysis of myocardial and circulating small RNAs in human heart failure and their utility as biomarkers

Circulating biomarkers to monitor cancer progression and treatment

Integrative single-cell and cell-free plasma RNA transcriptomics elucidates placental cellular dynamics

The potential of circulating cell free RNA as a biomarker in cancer

Noncoding RNAs in extracellular fluids as cancer biomarkers: the new frontier of liquid biopsies

The circulating transcriptome as a source of cancer liquid biopsy biomarkers

Analysis of cell-free Epstein-Barr virus-associated RNA in the plasma of patients with nasopharyngeal carcinoma

Characterization of extracellular circulating microRNA

Noninvasive in vivo monitoring of tissue-specific global gene expression in humans

Circulating RNAs as new biomarkers for detecting pancreatic cancer

Quantitation of circulating satellite RNAs in pancreatic cancer patients

Detection of circulating extracellular mRNAs by modified small-RNAsequencing analysis

Phospho-RNA-seq: a modified small RNA-seq method that reveals circulating mRNA and lncRNA fragments as potential biomarkers in human plasma

Ribonuclease activity in human plasma

Stability of endogenous and added RNA in blood specimens, serum, and plasma

Whole blood RNA sequencing reveals a unique transcriptomic profile in patients with ARDS following hematopoietic stem cell transplantation

Diverse human extracellular RNAs are widely detected in human plasma

Total extracellular small RNA profiles from plasma, saliva, and urine of healthy subjects

High-throughput sequencing of human plasma RNA by using thermostable group II intron reverse transcriptases

Assessing sample and miRNA profile quality in serum and plasma or other biofluids

Circulating microRNAs as stable blood-based markers for cancer detection

Comprehensive profiling of circulating microRNA via small RNA sequencing of cDNA libraries reveals biomarker potential and limitations

RNA-seq of human reference RNA samples using a thermostable group II intron reverse transcriptase

Thermostable group II intron reverse transcriptase fusion proteins and their use in cDNA synthesis and next-generation RNA sequencing

Template switching mechanism of a group II intron-encoded reverse transcriptase and its implications for biological function and RNA-Seq

Broad and adaptable RNA structure recognition by the human interferon-induced tetratricopeptide repeat protein IFIT5

Base-resolution mapping reveals distinct m 1 A methylome in nuclear-and mitochondrial-encoded transcripts

The m 1 A landscape on cytosolic and mitochondrial mRNA at single-base resolution

Broad role for YBX1 in defining the small noncoding RNA composition of exosomes

DMS-MaPseq for genome-wide or targeted RNA structure probing in vivo

Improved TGIRT-seq methods for comprehensive transcriptome profiling with decreased adapter dimer formation and bias correction

Simultaneous sequencing of coding and non-coding RNA reveals a human transcriptome dominated by a small number of highly expressed non-coding genes

Reducing the structure bias of RNA-Seq reveals a large number of nonannotated non-coding RNA

Facile single-stranded DNA sequencing of human plasma DNA via thermostable group II intron reverse transcriptase template switching

High-resolution characterization of sequence signatures due to non-random cleavage of cell-free DNA

Cell-free DNA comprises an in vivo nucleosome footprint that informs its tissues-of-origin

Effectiveness of convalescent plasma therapy in severe COVID-19 patients

tRFdb: a database for transfer RNA fragments

BCDIN3D regulates tRNAHis 3' fragment processing

Simultaneously monitoring immune response and microbial infections during pregnancy through plasma cfRNA sequencing

The dormant blood microbiome in chronic, inflammatory diseases

Numerous uncharacterized and highly divergent microbes which colonize humans are revealed by circulating cell-free DNA

Granzyme A in human platelets regulates the rynthesis of proinflammatory cytokines by monocytes in aging

A tour through the transcriptional landscape of platelets

Tissue-based map of the human proteome

Proinflammatory activities of S100: proteins S100A8, S100A9, and S100A8/A9 induce neutrophil chemotaxis and adhesion

Translational coregulation of 5′TOP mRNAs by TIA-1 and TIAR

Regulated capture by exosomes of mRNAs for cytoplasmic tRNA synthetases

The many pathways of RNA degradation

Regulation of cytoplasmic mRNA decay

Model-based analysis of ChIP-Seq (MACS)

The ENCODE blacklist: identification of problematic regions of the genome

Robust transcriptome-wide discovery of RNA-binding protein binding sites with enhanced CLIP (eCLIP)

Plasma proteome database as a resource for proteomics research: 2014 update

The major mRNA-associated protein YB-1 is a potent 5' capdependent mRNA stabilizer

Rfam 13.0: shifting to a genome-centric resource for non-coding RNA families

Detection of expanded RNA repeats using thermostable group II intron reverse transcriptase

Excised linear introns regulate growth in yeast

Detecting agotrons in Ago CLIPseq data

Argonaute-associated short introns are a novel class of gene regulators

Transcriptome-wide identification of RNA-binding protein and microRNA target sites by PAR-CLIP

Mammalian mirtron genes

The mirtron pathway generates microRNA-class regulatory RNAs in Drosophila

Intronic microRNA precursors that bypass Drosha processing

Genome-wide searching for pseudouridylation guide snoRNAs: analysis of the Saccharomyces cerevisiae genome

Human Tra2 proteins jointly control a CHEK1 splicing switch among alternative and constitutive target exons

The influence of LINE-1 and SINE retrotransposons on mammalian genomes

Prevalent, dynamic, and conserved R-loop structures associate with specific epigenomic signatures in mammals

The eight human "canonical" ribonucleases: molecular diversity, catalytic properties, and special biological actions of the enzyme proteins

5′ tRNA halves are present as abundant complexes in serum, concentrated in blood cells, and modulated by aging and calorie restriction

Exosome-mediated miRNA delivery promotes liver cancer EMT and metastasis

Exosome-mediated transfer of mRNAs and microRNAs is a novel mechanism of genetic exchange between cells

Limitation of alignment-free tools in total RNA-seq quantification

Atropos: specific, sensitive, and speedy trimming of sequencing reads

Fast gapped-read alignment with Bowtie 2

Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype

ART: a next-generation sequencing read simulator

Near-optimal probabilistic RNA-seq quantification

GFF utilities: GffRead and GffCompare. F1000Res

RSeQC: quality control of RNA-seq experiments

Improved metagenomic analysis with Kraken 2

Pavian: interactive analysis of metagenomics data for microbiome studies and pathogen identification

BEDTools: the swiss-army tool for genome feature analysis

Counting individual DNA molecules by the stochastic attachment of diverse labels

Repbase update, a database of repetitive elements in eukaryotic genomes

Infernal 1.1: 100-fold faster RNA homology searches

ViennaRNA package 2.0

Probabilistic programming in Python using PyMC3

GtRNAdb 2.0: an expanded database of transfer RNA genes identified in complete and draft genomes

PhastWeb: a web interface for evolutionary conservation scoring of multiple sequence alignments using phastCons and phyloP