Ancient RNA from Late Pleistocene permafrost and historical canids shows tissue-specific transcriptome survival

Abstract

While sequencing ancient DNA (aDNA) from archaeological material is now commonplace, very few attempts to sequence ancient transcriptomes have been made, even from typically stable deposition environments such as permafrost. This is presumably due to assumptions that RNA completely degrades relatively quickly, particularly when dealing with autolytic, nuclease-rich mammalian tissues. However, given the recent successes in sequencing ancient RNA (aRNA) from various sources including plants and animals, we suspect that these assumptions may be incorrect or exaggerated. To challenge the underlying dogma, we generated shotgun RNA data from sources that might normally be dismissed for such study. Here, we present aRNA data generated from two historical wolf skins, and permafrost-preserved liver tissue of a 14,300-year-old Pleistocene canid. Not only is the latter the oldest RNA ever to be sequenced, but it also shows evidence of biologically relevant tissue specificity and close similarity to equivalent data derived from modern-day control tissue. Other hallmarks of RNA sequencing (RNA-seq) data such as exon-exon junction presence and high endogenous ribosomal RNA (rRNA) content confirms our data's authenticity. By performing independent technical library replicates using two high-throughput sequencing platforms, we show not only that aRNA can survive for extended periods in mammalian tissues but also that it has potential for tissue identification. aRNA also has possible further potential, such as identifying in vivo genome activity and adaptation, when sequenced using this technology.

Fig 1. Regressions of ancient liver and historical skin samples, Method 1: Relationships between 95th percentile of expressed genes in each control tissue sample (x-axis) and each ancient sample or control samples from other tissues (y-axis).

Black points in graphs comparing ancient samples are the relationships between the control tissue and the equivalent ancient tissue. Red points overlaid show the relationship between the control tissue and other ancient tissues specified in the graph subtitle. Yellow lines are least squares linear regression fit for black points. Green lines are least squares linear regression fit for red points. Filled lines indicate a significant linear regression. Dashed lines indicate a nonsignificant linear regression. (A) BGISEQ-500 data, de-duplicated; (B) HiSeq-2500 data, de-duplicated; (C) BGISEQ-500 data, duplicates retained; (D) HiSeq-2500 data, duplicates retained. The underlying data for this figure are derived from Varistran output, summaries of which can be found in S2 Data and S3 Data.

Fig 2. Comparison of ancient and control tissues using Method 2.

Coverage scores (y-axis) were calculated based on the mean coverage of reads to each named gene in the CanFam3.1 transcriptome, followed by filtering to the 95th percentile of all genes represented. Each gene was then assigned a most-associated tissue based on data from an Affymetrix array derived from 10 canine tissues (x-axis). Each tissue was then assigned a cumulative score based on the coverage scores of each gene in the 95th percentile. Orange bars represent modern control tissues and blue bars represent ancient/historical tissues. (A) Historical Skin 2 versus control skin. (B) Ancient Tumat liver versus control liver. The underlying data for this figure can be found in S7 Data under Tissue_summary tabs, and are derived from primary data found in S4 Data under the Skin2_HS and Tumat_liver_HS tabs and from S5 Data under the ctrl_skin and crtl_liver tabs.

Fig 3. Regressions of all samples, Method 2: Relationships between the 95th percentile of expressed genes in ancient tissues (x-axis) versus control samples (y-axis).

Values are calculated based on per-tissue scores (see Methods), having removed duplicate reads from mapping data. Black data points and trend line refer to BGISEQ-500 data, while orange data points and trendline refer to Illumina HiSeq-2500 data. (A) Skin 1, (B) Skin 2, (C) Tumat cartilage, (D) Tumat liver, and (E) Tumat muscle. The underlying data for this figure can be found in S7 Data, ‘Regressions’ tab.