Ancient microRNA profiles of a 14,300-year-old canid samples confirm taxonomic origin and give glimpses into tissue-specific gene regulation from the Pleistocene.

DNA sequencing is the current key technology for historic or ancient biological samples and has led to many exciting discoveries in the field of paleogenomics. However, functional insights into tissue identity, cellular composition or gene regulation cannot be gained from DNA. Recent analyses have shown that, under favorable conditions, RNA can also be sequenced from ancient samples, enabling studies at the transcriptomic and regulatory level. Analyzing ancient RNA data from a Pleistocene canid, we find hundreds of intact microRNAs that are taxonomically informative, show tissue-specificity and have functionally predictive characteristics. With an extraordinary age of 14,300 years, these microRNA sequences are by far the oldest ever reported. The authenticity of the sequences is further supported by a) the presence of canid / Caniformia-specific sequences that never evolved outside of this clade, b) tissue-specific expression patterns (cartilage, liver and muscle) that resemble those of modern dogs and c) RNA damage patterns that are clearly distinct from those of fresh samples. By performing computational microRNA-target enrichment analyses on the ancient sequences, we predict microRNA functions consistent with their tissue pattern of expression. For instance, we find a liver-specific microRNA that regulates carbohydrate metabolism and starvation responses in canids. In summary, we show that straightforward paleotranscriptomic microRNA analyses can give functional glimpses into tissue identity, cellular composition and gene regulatory activity of ancient samples and biological processes that took place in the Pleistocene, thus holding great promise for deeper insights into gene regulation in extinct animals based on ancient RNA sequencing. . Published by Cold Spring Harbor Laboratory Press for the RNA Society.

INTRODUCTION

Sequencing historic and ancient DNA from up to 780,000-yr-old samples has become a standard approach to infer genetic histories of extinct or extant species (Orlando et al. 2013). While much has been learned from studying these genomes, it remains difficult to infer cellular processes, such as in vivo genome function, from DNA directly. In contrast, the presence and relative expression levels of RNA convey genetic information at the functional level, and can thus give insights into the biological activity, the tissue identity, or even the cellular composition of samples (Newman et al. 2015). The prospect of studying gene-expression, gene-regulation, or species-specific transcripts in extinct species would be a fantastic next step in the paleogenetics field and, at least theoretically, could enable the discovery of extinct genes, typically overlooked by aDNA-based approaches. However, given the putative relative lack of surviving material in historic and ancient samples, due to the release of RNases during decomposition in most tissues (Huynen et al. 2012), “ancient RNA” (aRNA) has rarely been studied and is often disregarded as a source for biological discoveries. On the other hand, preservative conditions for aDNA, such as cold and drought, also inhibit RNases, and many types of RNA molecules are present in hundreds of thousands of copies per cell (see Calabrese et al. 2007 for microRNAs).

With the availability of powerful and sensitive sequencing approaches, a small number of recent studies has capitalized on early observations of very short RNA fragments surviving in some archaeological plant and mummified materials (Rollo 1985; Venanzi and Rollo 1990) and have recovered sequenceable amounts of ancient RNA (aRNA) from historic plant and feces samples. In doing so, they have detected, for instance, viral RNA genomes (Ng et al. 2014; Smith et al. 2014), historical plant microRNAs (Smith et al. 2017) and fragments of plant protein-coding transcripts (Fordyce et al. 2013). While initially only qPCR-based methods were used on human remains, such as the more than 5000-yr-old Tyrolean iceman (Keller et al. 2017), a recent aRNA sequencing study detected disease-related microRNAs from about 1000-yr-old human skeletons (Shaw et al. 2019). While these studies represent a proof of concept for aRNA sequencing, or paleotranscriptomics, they are limited to relatively recent (∼1000-yr-old) samples and, except for one study, restricted to plant seed-material that conveys highly favorable preservative conditions.

A recent sequencing study on three permafrost tissues samples from the ancient “Tumat puppy” (liver, cartilage, and muscle) and two historical canid museum samples (95 and 151-yr-old wolf-skins) showed that aRNA can survive for extended periods in mammalian samples (Smith et al. 2019). With approximately 14,300 yr of age for the Pleistocene permafrost samples, the authors presented the oldest ever sequenced RNA to this date. In a paleotranscriptomic approach, the presence of exon–exon junctions of protein-coding genes and high levels of ribosomal RNA were detected, confirming the data's authenticity. Using two custom bioinformatics methods based on the quantitative comparison of detected short fragments of protein-coding genes between ancient and recent tissue sample profiles, the aRNA profiles could resolve tissue identity for two of the five samples (Smith et al. 2019).

MicroRNAs are ∼22 nt short RNA molecules that are key regulators of protein coding genes with important functions in numerous biological processes including development and disease (Bartel 2018). The so-called “mature microRNAs” and their usually nonfunctional by-products (“star microRNAs”) derive from longer RNA hairpin structures (“precursor microRNAs”) by means of several well-studied processing steps that leave characteristic patterns in sequencing data, which are useful for their annotation (Friedländer et al. 2012; Fromm 2016). Many microRNAs also exhibit tissue-specific expression, making them excellent cell- and tissue markers (Christodoulou et al. 2010; de Rie et al. 2017; McCall et al. 2017), also for difficult samples (Xi et al. 2007; Liu et al. 2009; Courts and Madea 2010; Glynn 2019). Since it is well-established that the number of times that a microRNA is detected by sequencing scales with its molecular abundance in a sample—with certain caveats (Linsen et al. 2009)—it is possible to apply sequencing to quantitatively measure microRNA compositions. Because microRNAs are themselves molecules with regulatory functions, detecting them in the ancient Tumat puppy tissues would give direct insights into gene regulation in the Pleistocene. Because the applied RNA extraction method (adapted from the Ambion miRvana kit) and library preparation strategy (NEBNext Multiplex Small RNA Library Prep Set) clearly favors short RNAs (for details see Smith et al. 2019), we were confident of detecting microRNAs among the sequenced aRNA transcripts from the Tumat puppy (Smith et al. 2019), thus enabling us to expand previous findings beyond protein-coding genes into the gene-regulatory level of ancient samples.

Furthermore, microRNAs are among the most conserved elements in metazoan genomes, giving them great potential as phylogenetic (Sempere et al. 2006; Tarver et al. 2013; Kenny et al. 2015) and taxonomic markers (Fromm et al. 2014). We have recently curated the microRNA complements of 45 Metazoan organisms, including canid, in the curated microRNA gene database MirGeneDB (Fromm et al. 2020) and shown that species-specific microRNAs can be used to trace the organismic origin of samples using our user-friendly software miRTrace (Kang et al. 2018). We applied these two resources to profile and characterize numerous full-length microRNAs from the ancient and historic canid samples. We found that many of the microRNA transcripts are canid-specific or have nucleotide variants that are canid-specific. We also found that the ancient and historical tissues have microRNA profiles that closely resemble their modern day counterparts (cartilage, liver, and muscle). Importantly, the ancient and historical microRNAs have nucleotide substitutions that are clearly distinct from microRNAs from fresh samples, and that are indicative of nucleotide damage over time. The combination of canid-specific sequences, characteristic tissue expression and RNA damage patterns strongly support the authenticity of the sequences. Predicting microRNA targets using only genome sequence, we were able to recapitulate microRNA functions that are consistent with their tissue-specific expression profiles. In summary, using careful statistical analyses and comparisons to recent microRNA tissue-atlas of dog, we show that these microRNAs are informative for ascertaining their authenticity and taxonomic origin, that they predict the tissue identity of all analyzed ancient samples and that they can further provide glimpses into the gene regulation of animals from the Pleistocene.

RESULTS

Detection of microRNAs from ancient and historic samples

Given their documented features, such as short length, frequent association with proteins and high abundance in cells, microRNAs are typically very stable in RNA preparations and thus should be detectable in aRNA samples (Jung et al. 2010). When reanalyzing the smallRNAseq data sets from Smith et al. (2019), we indeed found that of the 111,161,070 sequence read-outs, or “reads” that passed quality controls (Materials and Methods) in the historic and ancient samples, more than ∼0.017% represented microRNA reads (Fig. 1A; Supplemental File 1 miRTrace QC report; Supplemental File 2). Specifically, we found that not only the historic samples (0.023 and 0.045% microRNAs, respectively), but also the much older ancient liver sample (0.0035% microRNA reads) exceeded the percentage of microRNA reads previously reported from the only other NGS data sets on old RNA available (0.003% in 500-yr-old historic barley seed from Smith et al. 2017; and 0.002% in ∼1000-yr-old human bone sample from Shaw et al. 2019, respectively; Fig. 1A; Supplemental Files 3, 4 miRTrace QC reports). The fact that the microRNA content of the Tumat puppy samples and the historic wolf samples was in the same order of magnitude as that of the historic barley and the historic human sample further supports the authenticity of the samples.

FIGURE 1.

Enrichment and damage patterns of historic and ancient RNA samples. (A) microRNA proportion on the total reads and qc passed reads in ancient and historic samples is comparable to the historic barley smallRNA sequencing data from (Smith et al. 2017) and the historic human skeleton smallRNA sequencing data from Shaw et al. (2019). (B) RNA damage patterns increase over time. mapDamage profiles of modern (top), historic (middle), and ancient (bottom) canid tissue samples to canid microRNA loci. Y-axis shows the frequency of observed nucleotide changes and x-axis is the position in the read, respectively. Red lines show C > U differences, while blue lines show G > A differences. ø indicates the mean nucleotide substitution rate.

Ancient and historic samples show characteristic microRNA damage patterns

Nucleotide damage patterns based on deamination events (C to U and G to A, respectively) are a hallmark of ancient DNA samples and can help to authenticate the age of samples (Pääbo 1989; Ginolhac et al. 2011). Smith et al. observed transcriptome-wide RNA damage patterns in the canid data, with prevalent nucleotide substitutions that were not particularly biased toward putative deaminase events (Smith et al. 2019) and we asked if these patterns could also be recovered in microRNA sequences. Indeed, characteristic damage patterns were also observed for microRNAs from historic and ancient samples when compared with microRNAs of modern samples (Fig. 1B). While modern samples (merged liver, muscle, skin) showed virtually no damage patterns (0.2% nucleotide substitutions, Fig. 1B top), historic samples (merged skin 1 and skin 2) bore clear indications of deamination events (2.9% mismatching reads, Fig. 1B middle). The ancient samples (merged ancient liver, cartilage, muscle) showed an additional increase in the number of nucleotide substitutions (3.1% mismatching reads, Fig. 1B bottom). These relative high rates of nucleotide substitutions are clearly distinct from the rates observed in modern samples, and therefore support the authenticity of the ancient microRNA sequences, thereby underlining the usefulness of nucleotide damage analyses for the authentication of RNA extractions.

Detection of hundreds of distinct microRNAs from ancient and historical samples

After detecting comparable proportions of microRNAs on the total read numbers and authentic damage patterns for microRNA loci, we next asked how many distinct microRNAs were detected in each sample and at which level. Across all samples, we detected 74% of all currently known microRNAs in the canid (334 out of the currently 447 annotated canid microRNA genes, Fig. 2A; Fromm et al. 2020), with more than 18,000 sequence reads (Fig. 2B), 13,255 and 4678 from the historic skin samples and 172 (liver), 40 (cartilage), and 24 (muscle) from the ancient samples, respectively (Supplemental File 4; Supplemental Table 1 for summary). The number of detected distinct microRNAs and corresponding reads was highest for the historic samples (258 and 225) and lower for the ancient samples (31 for ancient liver, 13 for ancient cartilage, 14 for ancient muscle) (Fig. 2B,C). Among the ancient samples, the liver showed the best preservation in terms of microRNA content, which is in line with the transcriptome-wide observations reported by Smith et al. (2019), and the experimentally observed postmortem near-sterile conditions in the liver of up to 5 d (Tuomisto et al. 2013).

FIGURE 2.

Abundant numbers of conserved and dog-specific microRNAs are detected in historic and ancient samples. (A) 330 microRNAs are detected in the historic and ancient samples, of those 303 are identical in sequence to human (gray), 11 microRNAs show single nucleotide differences (polymorphic) (green), and 16 microRNAs are absent in human (blue). (B) Total number of microRNA reads detected in each sample (yellow) (logarithmic scale). (C) Detailed numbers of conserved, polymorphic, and canid-specific microRNAs detected in the five ancient and historic samples (colors as in A). (D) Hairpin structures and indication of mature (red) and star (blue) products of six microRNAs specific to canid and currently not known in any other organism.

Canid-specific microRNAs support taxonomic origin of the sequences

Given the extraordinary age of the canid aRNA samples (∼14,300 BP), it is a relevant concern that the detected sequences may not originate from the ancient tissues, but rather from technical sequencing errors or trace levels of contamination from researchers or the environment. Because microRNAs have been shown to be excellent taxonomic markers (Fromm et al. 2014; Kang et al. 2018), we next asked whether the detected microRNAs might be informative and conclusive of a canid origin of the reads and samples, or if microRNAs would indicate any source of contamination from other biological material.

In a first step, our microRNA taxonomic tracing software miRTrace (Kang et al. 2018), which was specifically designed to detect the taxonomic rank of samples based on their distinct microRNA content, and which is very sensitive for contamination events, did not detect any microRNA-based contamination from non-canid eukaryotes such as human or rodents (Supplemental File 1). However, since the majority of microRNAs are highly conserved in their sequence during animal evolution (Fromm et al. 2015; Bartel 2018), the absence of microRNAs specific to putative contaminating organisms cannot be considered as conclusive of the taxonomic origin of the samples, that is, if the microRNA reads are in fact derived from a canid. In additional analyses, beyond the basic taxonomic reference in miRTrace, which currently does not include a canid reference, or the information of conserved microRNAs with species-specific changes in one or several nucleotides, we next investigated if the sequences of the detected microRNAs can be used to unambiguously identify the taxonomic origin of the samples. The accurate annotations of MirGeneDB for 45 metazoan taxa (Fromm et al. 2020) enabled us to make sequence comparisons for each microRNA to human and canid, and test for conservation and variation in the 334 individual canid microRNAs detected.

As expected, we found that the majority of detected microRNAs, in total 303, belong to microRNAs that are fully conserved between human and canid (Fig. 2A, gray bars). Thus, these microRNAs should not be considered informative for taxonomic delimitation within mammalians. Altogether 27 distinct microRNAs supported the canid origin of the data (Supplemental File 4). Among those, three distinct groups of taxonomically informative microRNAs could be distinguished: (i) conserved mammalian microRNAs that showed clear nucleotide differences in the canid in comparison to their orthologs in human (“polymorphic microRNAs”), (ii) conserved mammalian microRNAs that are found in the canid (and other mammals) but are absent from human (“not in human”) and, finally, (iii) microRNAs not (currently) found outside canid (“canid-specific”). Of the 11 observed polymorphic microRNAs (Supplemental File 4) compared with their human orthologous counterparts, six mature microRNAs showed also differences to other organisms (Supplemental Fig. 1A), which, in addition to not being derived from human, also makes an origin from other organisms unlikely. Of the seven conserved microRNAs that are completely missing from the human genome (“not in human,” Supplemental File 4), single nucleotide differences to at least one or several other species orthologous sequences were found in five of them (Supplemental Fig. 1B), making them informative for a distinction with canid too. Finally, we detected nine mature microRNAs from seven distinct microRNA precursors (Cfa-Mir-1838, Cfa-Mir-1840, Cfa-Mir-8859-P1, Cfa-Mir-8859-P2, Cfa-Novel-6, Cfa-Novel-12, Cfa-Novel-23) not found in any other MirGeneDB species except the canid (Canis familiaris), and thus appeared to be candidate canid-specific microRNAs (based on current canid microRNA annotations and the lack of other carnivore species in MirGeneDB), having never evolved in any other animal group (Fig. 2D; Supplemental File 4). This would not only render them candidate biomarkers for the authenticity of the samples (Fig. 2C,D), but also makes them stand out as interesting regulatory molecules with putative canid-specific functions.

Because the canid specificity is currently based on only one carnivore representative in our database, it is possible that these microRNAs already evolved in carnivores and are not canid-specific. Therefore, we used BLAT (Kent 2002) and mapped these seven “canid-specific” microRNA precursors against the 18 currently available carnivore genomes in Ensembl (Ensemble Release 101) (Yates et al. 2020). We found that none of the microRNAs are present in any cat-related carnivore and are, thus, unlikely to have evolved already at the carnivore node (Supplemental Fig. 2, red). Moreover, we found that Mir-1838 and Novel-6, which we detected in ancient liver and historic skin sample 1, respectively, are indeed canid-specific. The other five precursors (Cfa-Mir-1840, Cfa-Mir-8859-P1, Cfa-Mir-8859-P2, Cfa-Novel-12, Cfa-Novel-23) are found in the exact sequence, without any substitutions, in all canids, except for fox (Supplemental Fig. 2, green). Additionally, these five precursors are found, with several substitutions, in most of the other Caniformia species (Supplemental Fig. 2, light green), rendering them as Caniformia-specific microRNAs, which, given the excellent preservation of the clearly canid Tumat puppy and the clear taxonomic origin of the historic wolf skin samples, confirms the canid origin of these microRNAs for our samples.

When comparing the numbers of detected microRNAs belonging to the three distinct groups of taxonomically informative microRNAs between all samples (“polymorphic,” “not in human,” and “canid-specific”), we saw that not only the historic samples, but also the ancient specimens showed representatives of all groups of microRNAs, thus supporting the authenticity of the microRNA reads as being canid in origin (Fig. 2C).

While the presence of canid-specific microRNAs is difficult to explain in terms of technical artifacts, the polymorphic microRNAs could be, at least in principle, contamination events of human origin that resemble the dog sequences due to nucleotide substitutions caused by RNA damage or sequencing errors. Using binomial statistics, we calculate that the probability that 11 independent nucleotide conversions caused by RNA damage or sequencing errors would make human sequences appear canid would be exceedingly low (P < 0.000006, see Materials and Methods). In a separate unbiased analysis, we applied simple sequence matching to query the historical and ancient data for microRNAs that are specific to human, mouse, and dog, but again, we only found matches to the nine dog-specific microRNAs already reported (data not shown). Taken together, we did not find any sign of a technical or biological source of contamination, thus supporting our findings of the taxonomic authenticity of the ancient Tumat puppy samples.

Top expressed microRNAs inform about the cellular and tissue identity of ancient and historic samples

For extant organisms it has on many occasions been shown that a number of microRNAs have clear tissue-specific and sometimes even cell-type-specific expression patterns (Christodoulou et al. 2010; de Rie et al. 2017; McCall et al. 2017). Using custom methods for the coding part of the genome, Smith et al. (2019) demonstrated that total aRNA could be used to confirm the identity of ancient liver and one of the historic skin samples (Smith et al. 2019). We, therefore, next asked whether the tissue identity of the samples, and possibly insights into their cellular composition, could be inferred based on the sequenced microRNAs.

In a first step, because cell-type- or tissue-specific microRNAs usually take up a very large part of the expressed complement (Witwer and Halushka 2016), we checked the top five expressed microRNAs of each sample alone. We identified clear signals of tissue and cell-type-specific microRNAs (Supplemental File 4). Specifically, we identified Mir-205-P1, a skin- and cartilage-specific microRNA, as the microRNA with the highest abundance in both skin samples and ancient cartilage and Mir-203, with similar tissue specificity, in both historic skin samples (Teta et al. 2012). The muscle-specific “myoMir” Mir-133 (Sempere et al. 2004; McCarthy 2008) was found among the top five microRNAs in the ancient muscle sample, and the hepatocyte-specific Mir-122 (Lagos-Quintana et al. 2002; Chang et al. 2004; McCall et al. 2017) was among the top five microRNAs in the ancient liver, again confirming the identity of these tissues using microRNAs as markers.

Interestingly, we also noticed relatively high levels of microRNAs that are not strictly specific to any of the tissues (cartilage, muscle, liver, and skin), but are known to be specific to a range of immune cells (Table 1; McCall et al. 2017). Most notably, the lymphocyte-specific microRNA Mir-155 was detected in all specimens and the highest overall detected microRNA in ancient muscle. Mir-148-P1, which is mast-cell-specific, was detected in all but ancient muscle samples and highest in the ancient liver sample. The high levels of these microRNAs suggest a relatively high number of immune cells in the samples.

TABLE 1.

The relative rank, and actual detected reads in parentheses, of immune-cell-specific microRNAs detected in ancient and historic canid samples

Global microRNA expression of ancient and historic samples resembles modern patterns

Because the top microRNAs already indicated that each of the aRNA samples resembles, at least for one or two of the top expressed microRNAs, the expected microRNA enrichment given their putative biological origin, we next asked if the overall microRNA expression pattern could be informative, too. We therefore compared the canid aRNA samples with more than 66 high-quality smallRNAseq data sets of fresh soft dog tissues from two previously published studies (Supplemental File 5 for miRTrace QC report; Koenig et al. 2016; Penso-Dolfin et al. 2016). In a uniform manifold approximation and projection (UMAP) (McInnes et al. 2018), based on the normalized counts of microRNA sequence reads only, we found that the recent samples clustered well according to their tissue and organ identity (Fig. 3, round spots).

FIGURE 3.

MicroRNA profiles of ancient and historic samples resemble recent tissue expression and are informative for tissue identity. UMAP representation of a total of 71 samples (66 recent, three ancient, two historic) show clustering of major organ groups (red—skeletal muscle and heart, pink—nervous system and thymus, orange—lung, green—kidney and pancreas, blue—skin, brown—gastrointestinal tract, yellow—liver). Historical and ancient samples (squares and diamonds) cluster closely with their tissue atlas equivalents. Note that cartilage is a skin derivative. Tissue abbreviations: (Sk) skin, (Br) brain, (Lu) lung, (Du) duodenum, (Ki) kidney, (Co) colon, (Te) testes, (SN) sciatic nerve, (SM) skeletal muscle, (Il/Je) ileum/jejunum, (Pa) pancreas, (Li) liver, (He) heart, (Pl) plasma, (Th) thymus.

In a second step we projected historical and ancient samples onto this UMAP, based on the previously learned features. Strikingly, historical skin samples are projected closely to the recent skin sample and all ancient samples cluster with related recent samples: Ancient liver is closest to modern day liver samples, ancient muscle is located next to recent skeletal muscle and heart muscle samples, and, finally, ancient cartilage colocates with the recent and historical skin samples (Fig. 3), thus reflecting the shared developmental origin of skin and cartilage. Finally, when plotting a selection of commonly known cell- and tissue-specific microRNAs on top of this UMAP, organ-specific expression patterns could be visualized that further supported the previously shown features of top expressed microRNAs for organ-specific expression (Supplemental Fig. 3).

These findings emphasize that not only individual, highly enriched microRNAs, but also the global microRNA expression “fingerprint,” are excellent tissue- and cell-type markers able to predict the tissue and organ identity of the sequenced samples, even for those being thousands of years old. These findings indicate strong potential for future applications of microRNAs to identify the tissue origin of biological samples of both recent and ancient origin without clear a priori designation, similar to what has been proposed for forensics, parasitology and food quality control (Courts and Madea 2010; Manzano-Román and Siles-Lucas 2012; Kang et al. 2018; Glynn 2019).

Subsampling confirms high degree of microRNA tissue specificity

Ancient and historical samples show overall few microRNA reads. To ensure that the expression profiles observed are reflecting the presented tissue specificity, not low sample complexity, we developed a test to account for low microRNA read numbers specifically in aRNA samples. We compared the aRNA samples to subsamples of comparable size (i.e., read counts of microRNAs) of recent tissue profiles (Fig. 4). Specifically, we tested microRNAs previously reported to be specific to the tissues probed in the ancient and historic tissue samples (muscle, cartilage, liver, and skin, respectively). We found that only subsamples of recent muscle-containing tissue such as heart, skin and skeletal muscle matched or exceeded the relative expression levels of the myoMir Mir-133_3p in the ancient muscle sample. Convincingly, only the liver samples, but not any other tissue, match or exceed the hepatocyte-specific Mir-122_5p levels of the ancient liver sample. While, as expected, only subsamples of recent skin samples match or exceed the Mir-27-P1_3p levels found in historical skin samples, no cartilage sample was available for comparison from recent dog. It was comforting to see, however, that Mir-205-P1 was found in ancient cartilage and in one of the historic skin samples (Fig. 4).

FIGURE 4.

Subsampling of marker microRNAs in reference samples confirms high degree of tissue specificity. Each column represents the result of 1000 subsamplings in which all tissues were up- or down-sampled to match the depth of the respective reference tissue. The length of each bar indicates how often these subsamplings equaled or exceeded the expression level of the indicated microRNA in the reference tissue. For the ubiquitously expressed microRNA Cfa-Let-7-P1a_5p, colon sample #1 was used as reference.

To further test the specificity of our approach, we compared the tissue-specificity findings to the expression of microRNAs either reported to be specific to tissues not included in our ancient and historic samples or known to be ubiquitously expressed in many tissues and cell types. Supporting our findings, microRNA specificity was also confirmed for tissues not included in our samples: Only subsamples originating from them, but not in the ancient or historic samples, match or exceed the levels of the respective selected sample. Further, ubiquitously expressed microRNAs can be found at high levels in many different recent tissues (Fig. 4, shows two examples: Mir-181 and Let-7).

Predicted functions of ancient microRNAs are consistent with their patterns of tissue expression

As microRNA themselves are active regulatory molecules, their presence and differential abundance in different tissues alone can directly provide insights into functional gene activity (Bartel 2018). It is possible to computationally predict the regulatory mRNA targets of microRNAs by sequence complementarity between the microRNA and the target (see Bartel 2009). We then analyzed these predicted target sets for their Gene Ontology (GO) enrichments which represent curated functional annotations for individual genes (Ashburner et al. 2000). Indeed, when we examined the annotated functions of the predicted canid mRNA targets of the tissue-specific microRNAs, we found enrichment for functional annotations corresponding to biological processes clearly related to the respective tissue origin (Supplemental Table 2). For instance, the predicted targets of Mir-203, which is highly abundant in the historic skin samples, are enriched for the GO term “epithelium and blood vessel development,” which is clearly skin-related. Similarly, the targets of Mir-122, which is abundant in ancient liver samples, were enriched for the GO term “carbohydrate metabolism” and the GO term “response to starvation,” representing well-established liver functions. Ancient cartilage-specific Mir-205-P1 targets genes that are enriched for the GO term “platelet-derived growth factor receptor,” which has been described to regulate chondrocyte proliferation (Kieswetter et al. 1997), and the ancient muscle-specific Mir-133 targets genes enriched for GO terms related to muscle contraction (“divalent inorganic cation transmembrane transporter activity”) and neurotransmitters (“GABA receptor binding”). A full list of all specific microRNA and the enrichment of functional annotations within their target gene populations can be found in Supplemental Table 2. The results of this novel approach are very promising for the interpretation of aRNA data and hold great potential for the usefulness of bioinformatics microRNA target prediction for otherwise lost or inaccessible samples such as the case for extinct animals.

DISCUSSION

With an estimated age of 14,300 yr, we present the oldest microRNA transcripts ever sequenced. Using publicly available tools, we report a variety of distinct microRNAs in three ancient and two historic samples that show authenticating microRNA abundances in comparison to previously published historic barley seed and human bone data and characteristic RNA damage patterns. We further find that the microRNA profiles are informative and conclusive of taxonomic origin for all samples. We show that top expressed microRNAs, and also the global microRNA complements of each sample, identify the tissue origin of all samples. Hereby, microRNAs clearly outperformed the previous mRNA-based approach (Smith et al. 2019). Regardless of the goal of the tissue-identification approaches, some microRNAs show cell-type-specific expression patterns and give additional insights into the cellular composition of the samples, including an intriguing signal for the presence of high numbers of immune-cells. Whether these immune-cell-specific microRNAs, or “immuno-Mirs,” are indicative of biological processes such as death and could be used as taphonomical markers, or possibly reflect a higher stability of these molecules, is an exciting new avenue that requires substantial additional work that might be linked to thanatotranscriptomic studies.

Our taxonomic and tissue-identity findings encouraged us to ask more functional questions to get insights or glimpses into the in vivo genome activity of extinct animals. By predicting the targets of the top microRNAs in the samples, we could indeed identify pathways and cellular processes that the samples were likely conducting, and that the microRNAs were regulating, 14,300 yr ago. While we are aware that the technical variation in aRNA samples may be high due to degradation due to their age, the combination of mRNA and microRNA analyses nevertheless showed that global and local RNA expression levels are representative of the given tissue. However, in comparison to analyses of ancient mRNA, microRNAs clearly outperformed in terms of results and ease-of-use. They are straightforward, requiring for basic analyses only simple sequence matching to known reference microRNAs, rather than to full reference genomes or transcriptomes. Furthermore, because microRNAs are captured with protocols that are optimized for highly degraded and fragmented RNA (Smith et al. 2019), microRNA analysis could be used as a complementary paleotranscriptomic approach to aRNA data in general. Importantly, since both conserved and species-specific microRNAs are detected well, and in a few cases even with mature and star reads, it would in principle be possible to reliably predict microRNAs de novo in extinct species, possibly discovering novel extinct microRNA families and genes that are not present in any living animal today (Ambros et al. 2003; Friedländer et al. 2008; Fromm 2016). Since the targets of these microRNAs can be predicted from the genome sequence, as exemplified in our findings above, this implies that we could predict putative functions of these extinct microRNAs that are specific to extinct animals, possibly getting insights into their physiology, ecology, and extinction events thousands of years after their death.

The detection of intact microRNAs in Pleistocene permafrost samples represents a proof of concept and opens up novel opportunities for future studies on the in vivo genome activity in other ancient samples, in particular those from extinct animals such as mammoths, Saber-toothed cats, cave lions and others that the increasingly melting permafrost might provide.

MATERIALS AND METHODS

Data acquisition and quality control

For RNA extraction and Next Generation Sequencing library preparation details see Smith et al. (2019). Illumina smallRNAseq data sets from ancient and historic as well as modern samples (Koenig et al. 2016; Penso-Dolfin et al. 2016) were downloaded from SRA using sra-toolkit 2.9.2 (Leinonen et al. 2010), processed using miRTrace (Kang et al. 2018), and quantified using the “quantifier” function of miRDeep2 (Friedländer et al. 2012) with the MirGeneDB 2.0 dog microRNA complement as reference (Fromm et al. 2020).

RNA damage analyses

Nucleotide damage analyses were performed using the mapDamage package (Jónsson et al. 2013) (mapDamage -i file.bam -r ref.fasta -l 20 -m 20 -y 0.4 –no-stats) after mapping the merged smallRNA sequencing read data to the canid mature microRNAs using bowtie 1 software (Langmead et al. 2009) and allowing for up to three mismatches in a 14 nt long seed alignment (bowtie -n 3 -l 14 –norc -y -a –best -S -f). Specifically, recent skin, liver, and muscle samples, historic skin 1 & 2 and ancient liver, muscle, and cartilage samples were merged, respectively.

Taxonomic determination

The internal miRTrace (Kang et al. 2018) reference of taxon-specific microRNAs was used to check for contaminations and to determine the likely taxonomic origin of the samples. Downstream sequence comparisons to human and 20 other vertebrate microRNA complements are based on MirGeneDB 2.0 annotations (Fromm et al. 2020) and the consistent nomenclature of microRNA gene orthologs, paralogs, and families (Fromm et al. 2015). Alignments were checked using custom scripts and AliView alignment viewer (Larsson 2014). Polymorphic reads were further checked using fastx-toolkit (http://hannonlab.cshl.edu/fastx_toolkit/).

Confirmation of canid specificity of MirGeneDB annotations

The microRNA complement annotations of carnivores in MirGeneDB, but also other microRNA databases such as miRBase (Kozomara et al. 2019) and miRCarta (Backes et al. 2018), currently only consist of dog. To test, therefore, whether the 11 “canid-specific” mature microRNAs (from seven precursors) detected in ancient and historic samples (Cfa-Mir-1838, Cfa-Mir-1840, Cfa-Mir-8859-P1, Cfa-Mir-8859-P2, Cfa-Novel-6, Cfa-Novel-12, and Cfa-Novel-23) were indeed canid-specific and not of deeper phylogenetic origin, we used BLAT (Kent 2002) of the precursor sequences on the Ensembl webserver with default parameters (exception 1 × 10−10 maximum E-value) choosing all 18 available carnivore genomes (see Supplemental Fig. 3 for list) (Yates et al. 2020).

Tissue determination

Small transcriptome-wide comparison of ancient and historic samples with modern samples were conducted using the uniform manifold approximation and projection for dimension reduction (UMAP) approach (McInnes et al. 2018) and were based on microRNAs that were detected in the ancient liver sample—the most complex of the ancient samples with regard to microRNA expression. Modern samples were excluded if they were flagged as low quality based on their mapping statistics or read length distribution by using the miRTrace software. Subsequently, ancient and historical samples were projected onto the UMAP based on the features learned from these tissue atlas samples.

Subsampling analysis

To ensure that the expression profiles observed are reflecting tissue specificity, not low sample complexity, we compared them to subsamples of comparable size of modern tissue profiles. In each case, one sample and one microRNA were selected as reference. All other samples were randomly subsampled 1000 times to the same number of reads as the selected reference sample for direct comparison. We then counted how many times this subsampling resulted in similar or higher levels of the respective microRNA compared to the reference sample.

Binomial testing

Since the 11 microRNAs that have nucleotide differences between dog and human are represented by 347 distinct sequencing reads in our data, we can apply binomial statistics to estimate the probability that these reads have been converted from human to dog versions by sequencing errors. If we assume that the 347 reads have been converted by independent sequencing errors, the probability will be around 1/(3347) = 2.7 × 10−166. This is assuming that 347 sequencing errors will all be of the exact type of nucleotide conversion that will make a human sequence appear as the dog sequence. Importantly, this is a conservative estimate, since it assumes the occurrence of a sequencing error in each read at the exact dog-human informative position, and with current Illumina instruments, sequencing errors are rare. If we rather assume that these reads are the result of 11 nt conversions caused by the long-term storage in ice, and that these 11 molecules have been PCR amplified to 347 reads, then the probability is rather 1/(311) = 0.0000056. Note that these analyses do not make any assumptions about the rate of sequencing errors or conversion events; they only consider the probability that errors or conversions are of the exact type that would make a human sequence look like a dog sequence. Hence, we here underestimate how unlikely it is that our analyses are confounded by sequencing errors or conversion events.

Target prediction and functional analyses

Tissue-specific microRNAs were defined as those that are abundant in the tissues of interest but not abundant in other tissues included in this study. Predicted mRNA targets of these microRNAs were downloaded from TargetScanHuman (Agarwal et al. 2015) using the lift over to dog transcripts. Gene ontology and KEGG pathway enrichment analyses were performed on the top 400 targets of these microRNAs using the R packages topGO (Alexa et al. 2006) and clusterProfiler (Yu et al. 2012). Relevant annotations were selected from the (up to) 25 highest ranking significant terms.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.