Diverse human extracellular RNAs are widely detected in human plasma

Abstract

There is growing appreciation for the importance of non-protein-coding genes in development and disease. Although much is known about microRNAs, limitations in bioinformatic analyses of RNA sequencing have precluded broad assessment of other forms of small-RNAs in humans. By analysing sequencing data from plasma-derived RNA from 40 individuals, here we identified over a thousand human extracellular RNAs including microRNAs, piwi-interacting RNA (piRNA), and small nucleolar RNAs. Using a targeted quantitative PCR with reverse transcription approach in an additional 2,763 individuals, we characterized almost 500 of the most abundant extracellular transcripts including microRNAs, piRNAs and small nucleolar RNAs. The presence in plasma of many non-microRNA small-RNAs was confirmed in an independent cohort. We present comprehensive data to demonstrate the broad and consistent detection of diverse classes of circulating non-cellular small-RNAs from a large population.

Figure 1. Sequencing data analysis using the Genboree-sequencing pipeline.

Small-RNAseq reads were processed and quantified using the exceRpt tool available on the Genboree Workbench. ExceRpt incorporates several modifications to existing analysis methods used to assess cytosolic miRNAs that specifically address experimental issues pertinent to exRNA profiling, such as variable contamination of ribosomal RNAs, and the presence of endogenous non-miRNA small-RNAs. Briefly, the software processes each sample independently through a casade of read-alignment steps designed to remove likely contaminants and endogenous sequences before aligning to exogenous miRNAs.

Figure 2. Overlapping expression of piRNAs and snoRNAs between post-MI and FHS Cohorts.

Using the same detection threshold as in our discovery phase (average RPM≥1 K combined), the replication sample detected 111 piRNAs and 164 snoRNAs with statistically significant correlation between the RPM values for snoRNAs (Spearmans ρ=0.32, P=0.02) and piRNAs (Spearman's ρ=0.53, P=0.01) between the FHS and BIDMC cohorts.

Figure 3. Distribution of plasma exRNA expression in 2,763 Framingham Heart Study participants.

From the RNAseq data we developed a list of observed exRNAs to target using high-throughput RT-qPCR (Fluidigm BioMark system). The expression of 331 human miRNAs, 97 human piRNAs, and 43 human snoRNAs were measured (Cq<23). Because piRNAs have a 3′ modification (2′-O-methylation), we utilized a modified RT-qPCR approach. The abundance and distribution of each type of small-RNA in the FHS cohort is illustrated in Figure 3 and shows that 50% of the FHS cohort expressed at least 11 snoRNAs and 20 piRNAs in plasma. These results demonstrate considerable variability in terms of number of exRNAs detected in plasma in each participant as well as the overall percent detected.

Figure 4. Measurement of plasma miRNAs and piRNAs by RT-qPCR in 2,763 FHS participants.

The prevalence and quantity of expression for each miRNA (a) and piRNA (b) are measured. RNAs are sorted from least to most expressed with the N and per cent expressed (that is, Cq<23). The quantification cycle (Cq) axis is sorted in reverse order so that expression increases from bottom to top (a) or from left to right (b) along the axis. The mean (±s.d.) Cq values are depicted with the horizontal capped lines and small squares (a,b) or dots (b) for each RNA. Light horizontal grey gridlines stratify RNAs by per cent of expression at <10% and >90%. Red squares or dots are used for statistically significant associations for age (P<0.05 after adjusting for multiplicity using the FDR) and navy squares are used for RNAs not meeting statistical significance. All supporting values included age and sex of participants are included in Supplementary Tables 1–3 and Supplementary Data 1–3.

Figure 5. Measurement of plasma snoRNAs and association with age and sex by RT-qPCR in 2,763 FHS participants.

(a) The prevalence and quantity of expression for each snoRNA are measured. RNAs are sorted from least to most expressed with the N and per cent expressed (that is, Cq<23). The quantification cycle (Cq) axis is sorted in reverse order so that expression increases from left to right along the axis. The mean (±s.d.) Cq values are depicted with the horizontal capped lines and small squares for each RNA. Light horizontal grey gridlines stratify RNAs by per cent of expression at <10% and >90%. (b) Range of fold-change values associated with age (median split at 66.5 years), sex and cardiovascular risk using the Framingham Risk Score. Red squares or boxes are used for statistically significant associations (P<0.05 after adjusting for multiplicity using the FDR) and navy squares and/or boxes are used for RNAs not meeting statistical significance. There were no statistically significant associations with sex for any of the RNAs assessed. In contrast, 39% (n=133) of miRNAs, 13% (n=13) of piRNAs and 21% (n=9) of snoRNAs were associated with age. As shown in Figure 4b and 5a, all piRNAs and snoRNAs and all but three of the miRNAs showed increased expression in subjects <66.5 years of age. All supporting values included age and sex of participants are included in Supplementary Tables 1–3 and Supplementary Data 1–3.

Figure 6. MicroRNA target prediction and pathway enrichment analysis.

(a) Venn diagram showing the overlaps of the predicted sets of unique mRNAs targeted by miRNAs for each factor. (b) Significantly (Benjamini Hochberg corrected P value<0.05) enriched KEGG pathways for target genes for factors 1–4. The number of miRNA involved in targeting at least one gene in a given pathway are shown in parentheses.