Environmental genomics of Late Pleistocene black bears and giant short-faced bears

Abstract

Analysis of ancient environmental DNA (eDNA) has revolutionized our ability to describe biological communities in space and time,^1-3 by allowing for parallel sequencing of DNA from all trophic levels.^4-8 However, because environmental samples contain sparse and fragmented data from multiple individuals, and often contain closely related species,⁹ the field of ancient eDNA has so far been limited to organellar genomes in its contribution to population and phylogenetic studies.^5,6,10,11 This is in contrast to data from fossils^12,13 where full-genome studies are routine, despite these being rare and their destruction for sequencing undesirable.^14-16 Here, we report the retrieval of three low-coverage (0.03×) environmental genomes from American black bear (Ursus americanus) and a 0.04× environmental genome of the extinct giant short-faced bear (Arctodus simus) from cave sediment samples from northern Mexico dated to 16-14 thousand calibrated years before present (cal kyr BP), which we contextualize with a new high-coverage (26×) and two lower-coverage giant short-faced bear genomes obtained from fossils recovered from Yukon Territory, Canada, which date to ∼22-50 cal kyr BP. We show that the Late Pleistocene black bear population in Mexico is ancestrally related to the present-day Eastern American black bear population, and that the extinct giant short-faced bears present in Mexico were deeply divergent from the earlier Beringian population. Our findings demonstrate the ability to separately analyze genomic-scale DNA sequences of closely related species co-preserved in environmental samples, which brings the use of ancient eDNA into the era of population genomics and phylogenetics.

Keywords: American black bear; ancient environmental genomics; environmental DNA; genomics; giant short-faced bear; paleoontology; phylogeny; population genetics.

Graphical abstract

Figure 1. American black bear phylogeny

(A) Map showing the black bear samples used. (B) Principal component analysis using smartpca, which accounts for the high amount of missing data in the Mexican samples by projecting the ancient samples onto a PCA created from the modern samples. (C) Genetic Hamming distance of UE1212 to each of the modern samples on biallelic SNPs, scaled to account for missing data, mapped to a color scale, and plotted on a phylomap using a neighbor-joining tree of the modern samples (results for UE1210 and UE1605 are shown in Figures S3A and S3B). (D) Inferred admixture graph, using two polar bear genomes (STAR Methods) as an outgroup in our admixture analysis. All data were parsed and plotted using admixtools2. We determined seven best-fitting graphs with highly similar topologies and many shared characteristics. The best of these is shown here, with a score of 4.922, and with a worst excess f4 residual of −2.182 for the configuration (East,Kenai;Mexican,Polar), and the rest are shown in Figure S4. See also Figures S1–S4.

Figure 2. Working model of American black bear phylogeography

(A) Pre-LGM–LGM conditions, with the ice sheet extending at ~21.5 kyr BP, and the hypothesized refugia to which the pre-LGM black bear population was suppressed. (B) Post-LGM conditions, with gray arrows indicating the northward recolonization of ice-free areas. See also Figures S3 and S4 and Table S2.

Figure 3. Photographs and descriptions of the three specimens used to generate the giant short-faced bear (Arctodus simus) genomes

(A–C) YG 24.1, a complete cranium of a giant short-faced bear (A. simus) that was collected from Pleistocene age permafrost sediments exposed at a placer gold mine along Ophir Creek, Yukon Territory, Canada. All measurements on the cranium demonstrate this is a very small individual compared to other specimens of this species that have been described from the region. All sutures appear to be fused and all adult teeth are present and fully erupted, demonstrating this individual was an adult. As a high degree of sexual dimorphism has been demonstrated for A. simus, it is likely this small cranium represents an adult female. (D–F) YG 76.4, a complete radius bone from a giant short-faced bear (A. simus) that was collected from Pleistocene age permafrost sediments exposed at a placer gold mine along Hester Creek, Yukon Territory, Canada. The very large size of this radius precludes it from being any other large Pleistocene carnivore known from the region. The bone exhibits a high degree of bone exostosis on major muscle attachments, suggesting this represents an older adult male. This radius (YG 76.4) articulates with specimen YG 129.1, a complete right ulna, which was collected at the same locality. (G) YG 546.562, a small fragment of a right femur diaphysis collected from Pleistocene age permafrost sediments exposed at a placer gold mine along Canyon Creek, Yukon Territory, Canada. The thick cortical bone wall and curvature of the diaphysis clearly compare well with those of a giant short-faced bear (A. simus). See Table S3 for measurements and radiocarbon ages. See also Figures S1 and S2 and Data S1.

Figure 4. Giant short-faced bear genomic and population estimates

(A) Biallelic transversion SNPs in UE1605, partitioned by read mapping (uniquely to the black bear mitochondrion, uniquely to Andean bear, or shared) and placed onto a mitochondrial Ursid tree. Lines above the black backbone lines of the tree indicate SNPs mapping uniquely to Andean bear; lines below the tree indicate mapping uniquely to black bear. The (+1) indicates a single supporting SNP in the black bear mapping leading to the Andean bear clade. (B) Phylogenetic tree and divergence times of the eight extant bear species and the extinct giant short-faced bear, as inferred from analysis of nuclear genomes. Branch lengths represent time before present (mya). The mean age of each node is shown, with 95% credibility intervals in parentheses and depicted as blue bars around each node. (C) PSMC plot for YG 546.562. See Figures S1 and S2 and Data S1.