Skip to main content
  • Research article
  • Open access
  • Published:

Moose genomes reveal past glacial demography and the origin of modern lineages

Abstract

Background

Numerous megafauna species from northern latitudes went extinct during the Pleistocene/Holocene transition as a result of climate-induced habitat changes. However, several ungulate species managed to successfully track their habitats during this period to eventually flourish and recolonise the holarctic regions. So far, the genomic impacts of these climate fluctuations on ungulates from high latitudes have been little explored. Here, we assemble a de-novo genome for the European moose (Alces alces) and analyse it together with re-sequenced nuclear genomes and ancient and modern mitogenomes from across the moose range in Eurasia and North America.

Results

We found that moose demographic history was greatly influenced by glacial cycles, with demographic responses to the Pleistocene/Holocene transition similar to other temperate ungulates. Our results further support that modern moose lineages trace their origin back to populations that inhabited distinct glacial refugia during the Last Glacial Maximum (LGM). Finally, we found that present day moose in Europe and North America show low to moderate inbreeding levels resulting from post-glacial bottlenecks and founder effects, but no evidence for recent inbreeding resulting from human-induced population declines.

Conclusions

Taken together, our results highlight the dynamic recent evolutionary history of the moose and provide an important resource for further genomic studies.

Background

Glacial and interglacial cycles had a major impact on the evolution of arctic and boreal species [1]. While many large mammals from northern latitudes went extinct toward the end of the Quaternary, several ungulates survived and subsequently recolonized the Holarctic [2]. For these species, rather than leading to extinction, glacial cycles induced range shifts and geographical isolation in refugia which resulted in inter-specific diversification or allopatric speciation [3].

Ungulates inhabiting high-latitude habitats including tundra and boreal forests are ideal models to study biogeographical processes such as vicariance and recolonisation associated with these glacial cycles. These habitats experienced important geographical shifts during glacial cycles [4, 5] and the high mobility of ungulates allowed them to quickly colonise or recolonise unglaciated areas (e.g. [6,7,8]). During the last glaciation, many Eurasian and North American ungulate species persisted south of the ice sheet or in southern refugia from where they recolonised areas following northward glacial retreat [6,7,8,9].

The effects of past climate changes induced contrasting responses among boreal and temperate ungulate species. For example, at the time of the Last Glacial Maximum (LGM), the range of the cold-adapted reindeer (Rangifer tarandus) consisted of both a continuous and large population, extending from Beringia and far into Eurasia that expanded during the Weichselian/Wisconsin glaciation some 115,000 years (115 ka) before present (BP), and of smaller refugial populations south of the ice sheets in Western Europe and in North America [6, 10]. Conversely, the population history of species from temperate climates, such as red deer (Cervus elaphus), fits a classical expansion-contraction model [1, 11], with contraction in isolated southern glacial refugia during cold periods followed by northward expansions during interglacials and in the Holocene [12].

Species inhabiting the taiga, such as moose (Alces alces), present adaptations typical of both cold-adapted and temperate species. The taiga is a boreal biome situated between the tundra and temperate habitats characterised by deciduous forests. Interestingly, variation in antler morphology of moose suggests phenotypic adaptation to open habitats, such as tundra, and to boreal forests [13]. Therefore, moose may have displayed an intermediate demographic response to past climatic changes between that of reindeer and of red deer, perhaps characterised by less severe demographic fluctuations.

Modern moose (Alces alces) first appear in the fossil record some 150–100 ka BP [14,15,16] and there is evidence that moose populations were negatively impacted by climatic changes at the end of the Pleistocene [17]. The glacial refugia moose population in Eurasia seems to have comprised two distinct genetic clades that evolved before the LGM and diversified afterwards [18]. Radiocarbon data further suggest that moose were absent from large parts of northern Europe during the LGM [19]. However, during that period, the distribution of moose extended as far south as the Italian Peninsula, the Balkans, and the Carpathians [17, 18, 20]. During the Holocene, Central and Eastern Europe were recolonised by moose from a glacial refugium in Eastern Europe and Scandinavia was recolonised from the south via a southern land bridge [20,21,22].

A northward shift of Eurasian boreal forests at the end of the LGM not only facilitated the recolonisation of higher latitudes but also allowed moose to colonize North America some 15–10 ka BP [19, 23,24,25], prior to the flooding of the Beringian land bridge some 14–11 ka BP [26]. Consistent with this hypothesis, there is no evidence of moose south of the Wisconsin glaciation ice sheets earlier than 15 ka BP in the fossil record [27]. However, the divergence between Yukon and British Columbia moose was estimated to 25–20 ka BP, and among other lineages to 17–13ka BP, suggesting that this colonisation may actually predate the end of the last glaciation [27].

While both European and North American moose populations display evidence of founder effects and bottlenecks associated with past climate changes, there is also evidence for recent anthropogenic impact on their genetic diversity [21, 28]. In Europe, the recent population history of the species was marked by a range-wide bottleneck dating back some 1800–1200 years BP and a more recent decline in moose numbers was also documented in the 18th to early 20th century [17]. Population declines during the Holocene contributed to shaping the current pattern of genetic diversity, where the Scandinavian population shows the lowest genetic diversity and highest inbreeding [17, 20]. Similarly, the species experienced a recent bottleneck in North America associated with a human-induced decline in the 1800s [27], which could potentially have exposed them to the negative effects of inbreeding.

Here, we generate the first de-novo reference genome for European moose and analyse it together with five previously published nuclear moose genomes to investigate patterns of genome-wide diversity and inbreeding and to test for evidence of founder effects, as well as Late Pleistocene and recent human-driven bottlenecks. We additionally generate five ancient moose mitogenomes and analyse them together with 16 previously published modern mitogenomes to study the phylogenetic relationships and divergence among Eurasian and North American moose lineages. Because of the differences in demographic history between European and North American moose as well as taxonomic uncertainties among lineages [28], this study fills an important gap in our understanding of moose evolution and population history.

Results

Genome assembly, genomic data and mitogenome reconstructions

The highest-quality and final de-novo assembly was generated with ALLPATHS_LG with a size of 2.48 Gb. It comprised 8373 scaffolds with a N50 of 1.7 Mb and had an average scaffold length of 296,974 bp. From 4104 mammalian single-copy orthologs, BUSCO showed the assembly to contain 111 (2.8%) missing, 121 (2.9%) fragmented, and 42 (1.0%) duplicated complete genes. In contrast, the assemblies generated with ABySS and SOAPdenovo had a scaffold N50 of 331.8 Kb and 316.3 Kb, respectively. Furthermore, the BUSCO analysis identified 164 (4%) and 317 (7.8%) missing, 221 (5.4%) and 256 (6.2%) fragmented, as well as 24 (0.6%) and 17 (0.4%) duplicated complete genes for the assemblies generated with ABySS and SOAPdenovo, respectively.

We identified 191 X chromosome-linked scaffolds representing 103 Mb (i.e. ~ 4.1% of the total assembly and ~ 67% of the size of the human X chromosome; Fig. S1). The average genome depth for the five nuclear genomes ranged from 12 to 20 (average = 16.6; Table S1) with the Swedish moose showing the highest depth (20-fold coverage). After filtering for missing and low-quality data, we obtained 3,204,006 high-quality SNPs segregating from the reference.

The final mitochondrial alignment was 16,693 bp long for a total of 21 mitogenomes. The newly sequenced ancient mitogenomes (n = 5) had a coverage ranging from 16.4 to 454.9 (Table S1).

Population structure, divergence estimates and demographic history

To examine the population structure of moose, we performed a Principal Component Analysis (PCA) for the five nuclear genomes and built a phylogeny for the 21 mitogenomes. The PCA based on the autosomal data indicated a clear distinction between Swedish and North American moose (Fig. 1). Also, the Eastern (A. a. americana) and Alaskan moose (A. a. gigas) showed slight genetic distinction from the other North American subspecies. Consistent with the PCA, the mitochondrial phylogeny showed a main European clade which included three previously described European sub-clades and a newly described clade composed of ancient specimens, as well as an Asian/North American clade, divided into two sub-clades (Fig. 2).

Fig. 1
figure 1

Principal Component Analysis (PCA) for five complete moose (Alces alces) genomes. The dataset comprised 3,204,006 high-quality SNPs

Fig. 2
figure 2

Geographical origin of 21 moose (Alces alces) specimens and Bayesian phylogeny based on 16,693 bp mitogenomes. a Circles and triangles depict modern and ancient specimens, respectively. b Asterisks depict samples for which dates were inferred in BEAST. Bayesian posterior probability support (over 0.9) for the branches is given. Divergence times (ka) between main lineages are given as the 95% HPD. Clades are labelled as: EU = European; A = Asian; NA = North American. Eastern, Central and Western European clades are labelled according to Świsłocka et al. [29]. The map was created using the ‘maps’ package (http://keithnewman.co.uk/r/maps-in-r.html) in R [30]

The mtDNA phylogenetic tree supported a divergence between the European and Asian/North American lineage dating back to 71 ka BP (95% Highest Posterior Density (HPD): 119–42 ka BP; Fig. 2b). The divergence between the Asian and North American lineages was estimated at 35 ka BP (95% HPD: 60–20 ka BP), and the divergence between the four European sub-clades at 25 ka BP (95% HPD: 43–16 ka BP; Fig. 2b). The estimated posterior substitution rate for the mitogenomes was 6.7 × 10− 8 substitutions/site/year (95% HPD: 3.53–9.84 × 10− 8).

We inferred changes in effective population size (Ne) of moose using the Pairwise Sequentially Markovian Coalescent (PSMC). Our results supported an increase in Ne for all moose populations during the Eemian interglacial, ca. 130–115 ka BP, followed by a decline during the last glacial period ca. 70 ka BP, and a small recovery at the end of the Pleistocene to a of Ne ~ 2000–8000 (Fig. 3, S2). While the PSMC curves indicated a shared demographic history for all moose until ca. 300 ka BP, the curve for the Swedish moose deviated from the North American ones at that time (Fig. 3).

Fig. 3
figure 3

Past demography for moose (Alces alces) using the PSMC. Each coloured line represents a different individual. The x-axis corresponds to time before present in years on a log scale, assuming an estimated substitution rate of 0.7 × 10− 8 substitutions/site/generation [31] and a generation time of 7 years [32]. The y-axis corresponds to the effective population size Ne

The Bayesian Skyline Plot (BSP) based on 14 European mitogenomes indicated a constant female effective population size (Nef) over the past 30 ka BP with a median of Nef ~ 71,400 assuming a generation time of 7 years (Fig. S3).

Heterozygosity and inbreeding

To compare genome-wide diversity in the five moose nuclear genomes, we estimated heterozygosity and inbreeding coefficients based on the identification of runs of homozygosity (FROH). The average genome-wide heterozygosity for our moose samples was 0.67 heterozygous sites per 1000 bp. The Swedish moose showed higher diversity (0.89) relative to the North American subspecies (average 0.62; Table 1).

Table 1 Genome-wide heterozygosity and inbreeding estimates (FROH) for five moose genomes

We found low to moderate inbreeding for both Swedish and North American moose, with FROH ranging from 8 to 23% of their genome in ROH (≥0.5 Mb; Fig. S4; Table 1). When considering only long ROH (≥2 Mb), mainly arising from recent mating among close relatives, we found low inbreeding coefficients ranging between 0.6 and 3.7%, with A. a. shirasi showing the highest values (Fig. S4; Table 1). The maximum ROH length was 5 Mb (A. a. shirasi 2; Fig. S4).

Discussion

A reference genome for moose

Moose are the largest and heaviest extant cervid species. They play an important ecological role in the boreal and temperate circumpolar forests of Eurasia and North America, with significant impacts on boreal forest regeneration and structure [33], soil fertility [34], and predator abundance (e.g. wolf, Canis lupus [35, 36]). Moose are also socioculturally and economically important in many regions [37]. In Fennoscandia for example, it is one of the most intensely managed species, with up to a third of its total population killed annually [21]. Here, we generated a de-novo genome for the species, with a quality on par to the de-novo genomes of most ruminant genomes currently available (8373 scaffolds and N50 = 1.7 Mb [31]). This reference genome will be an important resource for future studies of the species’ evolutionary history as well as its suite of adaptations to the boreal environment, but it can also be a relevant tool to further refine genetic methods for monitoring and management of moose to mitigate negative impacts on boreal forests [33], maintain healthy hunting stocks [38], and for safeguarding intraspecific diversity in line with intentions of the Convention on Biological Diversity (www.cbd.int [39]).

Moose origins and evolutionary history

Our results indicate a relatively recent divergence between European and Asian/North American moose lineages, dated to between 119 and 42 ka BP, based on modern and ancient complete mitogenomes. These dates are in agreement with an earlier divergence estimate of ca. 150–50 ka BP, based on short mtDNA data from 194 ancient radiocarbon dated samples [19], but at odds with the estimate of ca. 443 ka BP from Świsłocka et al. [29] based on complete mitogenomes. This inconsistency could be due to the fact that both our dataset and the one in Meiri et al. [19] included ancient radiocarbon dated specimens while in Świsłocka et al. [29] inferences are exclusively based on fossil-based calibration of the mutation rate based on modern data alone. Regardless, some of these divergence times have large confidence intervals, sometimes exceeding the fossil appearance of the species. The oldest A. alces remains are dated to 150 ka BP [15] and the latest fossils of earlier forms of Alces (i.e. Cervalces latifrons) date to ca. 186 ka BP [16, 40]. Therefore, further analysis of ancient nuclear genomes from A. alces and of its precursors would be needed to examine the precise timing of appearance of the European and Asian/North American moose lineages and of the origin of the species.

Glacial refugia and modern moose lineages

Previous genetic analyses supported the existence of moose glacial refugia during the LGM in the Alps, the Caucasus, Carpathians, Balkans and northern Italy as well as in western Siberia, the Ural Mountains and Russian plains (see [18] for a review and graphical representation). However, while the taiga biome was much reduced during the LGM in Eurasia, pollen and macrofossil data indicate that the taiga was characterised by a wide yet discontinuous geographical range in the form of isolated patches further east, in southern Ukraine and the Urals, Western Siberia, northern Mongolia and in eastern Siberia over most of the last glacial period [41, 42], potentially indicating the presence of several fragmented refugial moose populations throughout Eurasia. A higher genetic diversity in East Asia, estimated from ancient mtDNA, is also consistent with a LGM refugium in Siberia that could have extended to northern China [19]. Our dataset, combining ancient and modern mitogenomes, allowed us to identify one main European moose lineage composed of up to four clades diverging just before the LGM, ca. 25 ka BP, including three clades described in Świsłocka et al. [29] and one new, previously unknown clade. Sample MH315 from Tatarstan (1227 years BP) is part of the Eastern clade together with two samples from Poland. The Central clade is represented by three samples from eastern Poland. The modern Swedish moose belongs to the Western European clade. Even though limited to one sample, this supports the hypothesis that contemporary Scandinavian moose originate from a refugial population in western/central Europe via a southern colonization that happened after the last glaciation [43]. In fact, our oldest ancient moose sample from Germany, Meng37 (10 ka BP), is also part of this Western clade, indicating that this clade was present in central Europe at the onset of the Holocene. Finally, we describe a new European clade composed solely of three ancient moose samples from Germany dated to 1557–810 years BP (Meng43, Meng44, Qh37). While we caution about the accuracy of these dates, based on a limited number of radiocarbon dated specimens, it is intriguing that this clade is absent among modern moose, thereby suggesting that it disappeared during the Holocene. Such entire clade disappearance has previously been shown in brown bears (Ursus arctos [44]). Conversely, this clade indicates that there could have been more glacial refugia for moose in central Europe than previously identified.

Based on mtDNA data, we estimate a divergence between the Asian and North American moose lineages dating back to ca. 35 ka BP. This estimate predates the colonisation of the North American continent by moose, inferred from short mtDNA sequences and estimated at ca. 15 ka BP [28, 45] but it coincides with the divergence of human lineages in Asia prior to the colonisation of North America [46,47,48]. This suggests the existence of genetic structure in the Asian lineage prior to the North American colonisation via Beringia, where the first Alces fossils are dated to 15 ka BP [23]. Similar to Meiri et al. [19], our mitochondrial phylogenetic tree suggests that all modern North American moose derive from a single migration event. However, we find some genetic differentiation among North American subspecies in the PCA, with Eastern moose (A. a. americana) and Alaskan (A. a. gigas) moose being distinct from the other North American samples which form a tight cluster. We hypothesize that this differentiation could be a consequence of the founder effects during the colonisation of North America. Sequencing of additional ancient and modern genomes from eastern Siberia, Beringia and North America will be essential to solve these questions.

Past demography

Our results indicate that the past 150 ka of moose history was significantly influenced by glacial cycles. Overall, changes in Ne inferred with the PSMC, indicate that moose primarily responded to glacial cycles in the same way as temperate species such as red deer (C. elaphus), with expansions during interglacials and contractions in cryptic and isolated refugia during glacials [12]. For example, we found evidence for a demographic expansion possibly coinciding with the onset of the Eemian interglacial, which was followed by a decline ca. 70 ka BP. Interestingly, the lack of moose fossil remains from Europe during the LGM is consistent with either a demographic decline [19] or a shift in geographical distribution. Thus, it is possible that moose demographic history largely reflects historical shifts in boreal forests and the taiga during cold periods, including the last glaciation [19]. Following the decline, moose populations experienced a slight demographic recovery, some 15–10 ka BP at the Pleistocene/Holocene transition. Importantly, Meiri et al.’s [19] observation of a temporal gap in radiocarbon dates across Asia, but only two base pairs difference between two pre- and post-LGM moose from northern Siberia supports the hypothesis of habitat tracking in response to ice sheet movements. The moose population may thus have tracked its habitat as it shifted south during the LGM and then expanded north after the ice receded, instead of populations in the north going extinct and subsequently being replaced as climate conditions improved after the LGM. However, additional genomic data from diverse Eurasian moose populations would be required to further test this hypothesis.

There were striking differences between the demographic trajectories of Swedish and North American moose. While the Swedish moose experienced a gradual decline from the Eemian up to the Pleistocene/Holocene transition, the North American moose genomes showed signatures of increase in Ne during the Eemian and declined later, ca. 50 ka BP. The different trajectories between the two moose could be artefactual. However, since moose had not colonised North America at the time [28, 49], this difference could be an indication of different population histories between the European and Asian populations. For example, European moose seem to have experienced a stronger range contraction than the Asian moose during the last glaciation [19]. We also found evidence for population expansions in both European and North American moose at the end of the LGM, ca. 15 ka BP, which likely reflects the northward advances as the species tracked its habitat, similar to previous advances during warming periods, some 59 ka BP and 14 ka BP, as shown by Hundertmark et al. [28].

We caution that PSMC is known to lack power in estimating Ne in very recent times (i.e. < 10 ka BP) due to fewer coalescent events [50]. We bypassed this limitation by reconstructing female demographic changes using Bayesian inference based on the mitogenome data in European moose. We found a constant Nef over the past 30 ka BP, which is consistent with previous analyses of modern European mitochondrial data [18, 20]. Therefore, our results support a scenario where European moose experienced a spatial rather than a demographic expansion, with only limited founder effects, as they shifted their range northwards at the end of the LGM. Finally, while the confidence intervals for the Nef estimates were wide, we note a disconnect between a mitochondrial Nef of ~ 71,000 and a nuclear Ne of 2000-8000. This difference could be due to variance in male reproductive success in a highly polygynous species such as moose, with a mating system characterized by intrasexual competition, and where successful males maintain harems [51, 52].

Moose genomic diversity

We found that the majority of ROH were short (< 2 Mb) in all moose genomes, which indicates that most of the observed inbreeding is due to background relatedness [53, 54], potentially resulting from glacial or post-glacial bottleneck events [17, 20, 28, 43]. For example, the three American samples showed the highest background inbreeding and lowest diversity, consistent with a single or several founder effects when colonizing North America at the end of the last glaciation [27]. Interestingly, the Alaskan moose, A. a. gigas, had the lowest inbreeding and highest diversity of all North American moose samples. This could indicate that moose colonised Alaska from a large refugial population, possibly in Beringia and that they subsequently experienced serial founder effects in their colonisation of the rest of North America. In fact, both the Swedish and Alaskan moose displayed the highest levels of heterozygosity compared to the other moose specimens, suggesting the existence of large refugial populations in Europe and Beringia during the LGM [17, 20, 55].

We only detected low levels of inbreeding arising from long ROH (≥2 Mb), commonly caused by recent mating with relatives during bottlenecks [53, 54]. This result is at odds with previous studies which identified limited mitochondrial and nuclear diversity in moose populations from Europe [21, 28] and North America [27, 56]. Moreover, evidence for a severe reduction in effective population size to less than 3% of their former size (down to Ne ~ 400), lasting for hundreds of generations and possibly dating back as early as the fifteenth century, has been previously reported in Swedish moose populations [21]. Similarly, substantial reductions in Ne in moose populations from Lapland and northern Finland have been associated with 18th century bottlenecks [57]. Interestingly, Norwegian moose populations do not seem to have experienced such population reductions [58, 59]. It is thus possible that recent admixture with some of these populations may have prevented an increase in inbreeding in the Swedish population. Recent moose sightings in southeastern Germany and the expansion of populations from Poland further suggest that the species has potential for long range dispersal [60]. Altogether, our results suggest that the recent demographic history of European moose may be more complicated than that of a single bottleneck scenario.

While inbreeding was generally low, it is however worth noting that the two Rocky Mountain moose (A. a. shirasi) specimens had the largest inbreeding coefficients and the longest ROH (~ 4–5 Mb). This is consistent with these southernmost North American populations being founded from relatively few individuals after the declines of the late 1800s and with low gene flow from northern populations until the middle of the 20th century [27].

Conclusion

In this study, we present a reference genome for the European moose, which will serve as an important resource for monitoring as well as for future studies on the evolutionary history and population genomics of the species. Through analysis of several moose genomes, we provide a glimpse into the demography and population history of the species in Europe and North America. Our results indicate that current moose lineages trace back their origin to several refugial populations during the LGM. Throughout their history, moose have experienced similar population demographic fluctuations as temperate ungulates (e.g. red deer, C. elaphus [12]). However, their overall spatial response to the Pleistocene/Holocene transition are more consistent with a range shift, reflecting historical shifts in the taiga.

Methods

Sampling and data collection

We used stored muscle tissue from a young female moose from Sweden (Province of Gävleborg, Central Sweden; approximate RT90 coordinates 1570800/6855900) to generate a de-novo assembly (see below; Table S1). We also obtained bone samples from five ancient moose specimens from Russia (n = 1) and Germany (n = 4; Table S1). Radiocarbon dating of one German sample (Meng37) was performed at the Curt-Engelhorn-Centre for Archaeometry gGmbH (CEZA) laboratory at the Reiss-Engelhorn-Museen in Mannheim (sample identifiers starting with “MAMS”) using a MICADAS AMS system [61]. Conventional 14C ages were calibrated to calendar ages with the IntCal13 data set [62] and the SwissCal software (L. Wacker, ETH Zurich). Radiocarbon dating for the Russian sample (MH315) was performed at the Oxford Radiocarbon Accelerator Unit (sample identifiers starting with “OxA”). The AMS-date was turned into calendar years using the IntCal13 calibration curve in Oxcal v4.2 [63].

We then supplemented our newly generated data with published nuclear and mitochondrial genome data for four North American moose (ENA project number: PRJNA325061; Table S1) and 11 mitogenomes for Eurasian and North American moose (Table S1). We therefore used two different datasets: a genomic dataset comprising the nuclear genomes from five modern moose samples (one newly sequenced), and a mitogenomic dataset composed of 21 samples in total - 16 modern and five ancient ones (six newly sequenced; Fig. 2a, Table S1).

DNA extraction and sequencing

In order to generate a de-novo assembly, we extracted total DNA from muscle tissue from the female moose using a Thermo Scientific KingFisher Duo magnetic particle processor (ThermoFisher Scientific) with the KingFisher Cell and Tissue DNA Kit. We then prepared two paired-end libraries with 180 bp inserts using a TruSeq DNA kit (Illumina, CA, USA) according to the manufacturer’s specification, but automated on an Agilent NGS workstation (Agilent, CA, USA). We also constructed two mate-pair libraries according to the Nextera gel-plus protocol (Illumina, CA, USA) using a Blue Pippin (Sage Science Inc., MA, USA) for size selecting target fragments: one 2.5–4.5 kb mate-pair library (MPS) and one 5–7 kb mate-pair library (MPL) using the 0.75% 1–10 kb Gel Cassette with Marker S1. All libraries were indexed to enable de-multiplexing after sequencing. These libraries were sequenced on an Illumina HiSeqX v.2.5 instrument (HiSeq Control Software 3.3.76/RTA 2.7.6) with a 2 × 150 bp setup, where the 180 bp library was sequenced on 1.5 lanes. The mate-pair libraries were multiplexed in equimolar ratios and sequenced on 0.5 lanes (in total) at SciLifeLab (Stockholm, Sweden). BCL to FASTQ conversion was performed using BCL2FASTQ v1.8.3 from the CASAVA software suite.

DNA from the ancient sample from Russia (MH315; Table S1) was extracted using the silica-column protocol described in Yang et al. [64] and modified in Gamba et al. [65]. A USER treated genomic library [66] was then built following Meyer & Kircher [67] and subjected to mitogenome capture as described in Maricic et al. [68], using deer-specific baits and as described in Heino et al. [69]. The resulting library was purified, quantified on a 2100 Bioanalyzer (Agilent) and sequenced on a HiSeq lane with paired-end 2 × 126 bp setup.

DNA was extracted from the petrous bones of four additional ancient samples from Germany (Meng37, Meng43, Meng44, Qh37; Table S1) following Dabney et al. [70]. Single stranded libraries were prepared according to Gansauge et al. [71] and pre-treated with 0.5 μl of USER enzyme (New England Biolabs, Ipswich, MA, USA) for 15 min at 37 °C to remove uracil residues resulting from post-mortem damage (modified from Meyer et al. [72]). The resulting libraries were then amplified and indexed using two double-unique p5-p7 tailed primers to generate dual-indexed library molecules. The number of PCR cycles was determined in advance from qPCR results, as described in Gansauge & Meyer [73]. Amplified libraries were then pooled in equimolar amounts according to their concentration and length distribution determined using Qubit 2.0 and 2200 TapeStation (Agilent Technologies), respectively. Finally, libraries were sequenced on an Illumina NovaSeq6000 S4 platform using 2 × 100 bp paired-end setup at SciLifeLab (Stockholm, Sweden).

De-novo reference genome assembly

The raw sequencing reads were processed using Trimmomatic v0.32 [74] to remove low quality sequences and adaptor read-through. Trimmed reads were assembled using three different methods: ALLPATHS-LG r.52488 [75] with the option “HAPLOIDIFY = True”, ABySS v1.3.5 [76] and SOAPdenovo 2.04-r240 [77]. For assembly evaluation we used QUAST v4.5.4 [78] and BUSCO v3.0.2 [79] with the “mammalia_odb9” dataset. The assembly showing the highest contiguity and best BUSCO scores (Genbank project number PRJNA668262; Table S1) was used for mapping of resequenced data and downstream analyses.

We used the genomic dataset to identify the X chromosome-linked scaffolds of the moose reference genome using their relative genomic coverage as in Malde et al. [80]. Briefly, we estimated which scaffolds larger than 10 kb displayed ca. half of the coverage than the genomic average in the three male moose samples but ca. the average in the two female samples (Fig. S1).

Bioinformatic procedures

In order to generate complete nuclear moose genomes, we mapped paired-end short read data for the newly generated (ENA Project number: PRJNA40679; Table S1) and the four downloaded moose samples (ENA Project number: PRJNA325061; Table S1) to our de-novo moose assembly. To avoid biases, we processed both the newly generated and downloaded data using the same bioinformatic pipeline. Briefly, forward and reverse reads were trimmed to remove Illumina adapter sequences using Trimmomatic v0.32 with default settings [74] and then mapped to the reference genome using BWA mem v0.7.13 [81]. We then used SAMtools v1.8 [82] for coordinate sorting, indexing, and removing duplicates from the alignments. Next, we re-aligned reads around indels using GATK IndelRealigner v3.4.0 [83]. For all downstream analyses, we only retained reads with mapping quality ≥30. We estimated the depth of coverage for the five genomes using Qualimap v2.2.1 [84].

Next, we used bcftools v1.8 [82, 85] to call variants for each moose genome. We retained SNPs with a minimum depth of coverage of 1/3 of the average coverage and maximum depth of 10 times the average coverage of each individual genome. We also retained SNPs with base quality ≥30 and excluded those within 5 bp of indels. For all downstream analyses, we excluded all scaffolds (n = 191) identified as linked to X chromosomes. We also hard masked all repeat regions across the genome using BEDtools v2.27.1 [86]. Finally, we merged all five individuals into one single vcf file keeping only those positions genotyped in all samples using PLINK v1.9 [87].

For the mitochondrial dataset, we mapped the raw sequencing data from the ancient moose specimens sequenced here (Table S1) to the mitogenome of an Eurasian moose (Genbank accession number: MF784604.1). We used BWA aln v0.7.8 [81] with settings adjusted for ancient DNA studies as described in Pečnerová et al. [88]. After removing PCR duplicates with a python script which takes into account both the start and end of the alignments (github.com/pontussk/samremovedup), we generated consensus sequences from the BAM files in ANGSD v0.921 [89] using the majority rule (−doFasta 2) and a minimum depth threshold of 3. To avoid biases, in the four ancient moose samples subjected to partial USER treatment (Meng37, Meng43, Meng44, Qh37), which leaves post-mortem damage at the fragments termini, we trimmed 1 bp from both sides of each sequence before mapping.

The mitogenomes from the five modern moose for which shotgun sequencing data was available (newly sequenced and downloaded) were generated as above, but with the following modifications: mapping was performed with BWA mem, PCR duplicates were excluded using SAMtools, and the minimum depth threshold was set to 100. The final mitochondrial alignment of modern and ancient newly sequenced samples and downloaded sequences (n = 21) was generated in MAFFT v7.407 [90].

Population structure and phylogenetic analyses

We first used the R package SNPRelate [91] to perform a Principal Component Analysis (PCA) based on the genetic covariance matrix estimated from the genotypes using our filtered SNP genomic dataset.

Secondly, we used BEAUti/BEAST v1.8.4 [92] to build a phylogeny for the 21 mitogenomes. In order to estimate the molecular age of three undated samples (Meng43, Meng44 and Qh37), we created a ‘taxa’ partition of these undated samples. We set the median value of the dates for the two radiocarbon dated samples as calibrated tip dates (Table S1) and used the ‘sampling with individual priors’ for the undated samples partition (Meng43, Meng44, Qh37). We determined the best fitting evolutionary model for the mitogenome dataset with jModelTest v2.1.9 [93] and selected a model of HKY + G using the Bayesian Inference Criterion as previously used by DeCesare et al. [27]. We used a constant size model with a strict molecular clock and set the clock rate to a normal distribution with a mean value at 9 × 10− 9 substitutions/site/year based on Zurano et al. [94] and a standard deviation of 0.01. The ages of the three ancient samples without radiocarbon dates were determined using a uniform prior ranging from 0 to 100 ka BP. We then ran all models using BEAST v1.8.4 for 50 million generations with sampling every 1000 generations. To ensure that convergence had occurred, we visualized all output log files in Tracer v1.7.1 [95] and then combined all trees in a single tree with LogCombiner v1.8.4. We used TreeAnnotator v1.10.4 [96] to remove the first 10% of runs as burn-in from the tree files. Finally, we visualised and built the phylogenies in Figtree v1.4.4 (github.com/rambaut/figtree).

Demography and divergence estimates

To reconstruct past changes in effective population size (Ne) of moose, we used the Pairwise Sequentially Markovian Coalescent (PSMC v0.6.5 [50]) model and applied it to the five moose nuclear genomes. This model identifies historical recombination events across a diploid genome and infers the time to the most recent common ancestor (TMRCA) between independent segments of the genome. Assuming that pairwise sequence divergence is proportional to the time of the coalescent event, regions of low heterozygosity correspond to recent coalescent events while regions of high heterozygosity correspond to more ancient coalescent events. Because the rate of coalescence is inversely proportional to Ne, it can then be used to estimate temporal changes in Ne.

After excluding all sites in the scaffolds (n = 191) identified as linked to the X chromosome, we generated consensus sequences for all autosomes of the five modern genomes using the SAMtools mpileup [82] command and the vcf2fq command from vcfutils.pl. We then excluded sites with base quality and mapping quality below 30, and depth below 1/3 of the average coverage, which was estimated for each genome independently. Since PSMC’s model is highly sensitive to false heterozygous sites, for this analysis we used a more strict threshold of maximum depth excluding positions with more than two times the average coverage estimated for each genome. We set N (the number of iterations) = 25, t (Tmax) = 15 and p (atomic time interval) = 64 (4 + 25*2 + 4 + 6, for each of which parameters are estimated with 28 free interval parameters) for the inference of TMRCA between each chromosome from each individual genome. We scaled population parameters assuming a generation time of 7 years [32] and a substitution rate of 0.7 × 10− 8 substitutions/site/generation, inferred from a rate of 1 × 10− 9 substitutions/site/year from Chen et al. [31].

Given that PSMC results are less accurate at times more recent than 10 ka BP [50], we also estimated fluctuations in female effective population size (Nef) using the mitogenomes from 14 European moose using BEAUti/BEAST v1.8.4 [92]. We used the same parameters as described above for the phylogenetic tree but used a Bayesian Skyline Plot (BSP) model. We then visualised all output log files with Tracer v1.7.1 [95] to ensure that convergence had occurred and to generate the BSP.

Heterozygosity and inbreeding

In order to estimate the individual genomic diversity we used mlRho v2.7 [97]. mlRho allows to estimate the population mutation rate (θ), which approximates the per-site expected heterozygosity under the infinite sites model [97]. We excluded bases with quality below 30, reads with mapping quality below 30 and positions with root-mean-squared mapping quality below 30 from modern bam files and excluded all scaffolds identified as linked to X chromosomes. We also filtered out sites with depth lower than 1/3 of and higher than 10 times the average coverage. This allowed us to avoid biases associated with erroneous mapping to the reference genome and with false heterozygous sites caused by variable coverage across the genome which can result from structural variation.

We then quantified individual inbreeding coefficients by estimating the number and lengths of runs of homozygosity (ROH). ROH are long genomic tracts without heterozygous sites that can be used to inform about past and recent inbreeding levels [53, 54]. We first converted the filtered multi-individual vcf file comprising the five nuclear genomes into a ped file and identified all ROH in autosomal scaffolds. We then used PLINK v1.9 [87] to identify ROH with the following settings: a sliding window size of 1000 SNPs (homozyg-window-snp 1000); a maximum of 3 heterozygous sites per window (homozyg-window-het 3); a minimum of 5% of overlapping windows that must be called homozygous to define any given SNP as ‘in a homozygous segment’ (homozyg-window-threshold 0.05); a minimum of 10 SNP per window (homozyg-snp 10); a minimum of 500 kb coverage (homozyg-kb 500); a minimum SNP density of one SNP per 50 kb (homozyg-density 50); a maximum distance between two neighbouring SNPs of 1000 kb (homozyg-gap 1000); a maximum of 999 heterozygous sites within ROH of 999 (homozyg-het 999). Finally, we calculated individual inbreeding coefficients (FROH) as the overall proportion of their genome contained in ROH by summing the size of all ROH per individual divided by the total genome size (i.e. autosomes only). We considered two ROH size cut-offs at ≥0.5 Mb and ≥ 2 Mb, representing background relatedness from matings between distant relatives and recent inbreeding events, respectively [53, 54].

Availability of data and materials

de-novo reference genome for European Moose (A. alces): Genbank accession project number PRJNA668262 (Table S1).

Raw paired-end trimmed fastq data for one complete nuclear genome and five mitogenomes: ENA (https://www.ebi.ac.uk/ena) project accession number PRJEB40679 (Table S1).

Abbreviations

LGM:

Last glacial maximum

BP:

Before present

PCA:

Principal component analysis

HPD:

Highest posterior density

PSMC:

Pairwise sequentially markovian coalescent

TMRCA:

Time to the most recent common ancestor

BSP:

Bayesian skyline plot

N e :

Effective population size

N ef :

Female effective population size

ROH:

Runs of Homozygosity

FROH :

Inbreeding coefficients based on the identification of Runs of Homozygosity

References

  1. Hewitt G. Some genetic consequences of ice ages, and their role in divergence and speciation. Biol J Linn Soc. 1996;58:247–76. https://doi.org/10.1006/bijl.1996.0035.

    Article  Google Scholar 

  2. Geist V. Deer of the world: their evolution, behaviour, and ecology. Stackpole Books; 1998.

  3. Avise JC, Walker D. Pleistocene phylogeographic effects on avian populations and the speciation process. Proc Biol Sci. 1998;265:457–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Edwards ME, Brubaker LB, Lozhkin AV, Anderson PM. Structurally Novel Biomes: A response to past warming in BERINGIA. Ecology. 2005;86:1696–703. https://doi.org/10.1890/03-0787.

    Article  Google Scholar 

  5. Bigelow NH. Climate change and Arctic ecosystems: 1. Vegetation changes north of 55°N between the last glacial maximum, mid-Holocene, and present. J Geophys Res. 2003;108. https://doi.org/10.1029/2002jd002558.

  6. Flagstad O, Røed KH. Refugial origins of reindeer (Rangifer tarandus L.) inferred from mitochondrial DNA sequences. Evolution. 2003;57:658–70.

    Article  CAS  PubMed  Google Scholar 

  7. Latch EK, Heffelfinger JR, Fike JA, Rhodes OE Jr. Species-wide phylogeography of north American mule deer (Odocoileus hemionus): cryptic glacial refugia and postglacial recolonization. Mol Ecol. 2009;18:1730–45.

    Article  PubMed  Google Scholar 

  8. Sage RD, Wolff JO. Pleistocene glaciations, fluctuating ranges, and low genetic variability in a large mammal (Ovis dalli). Evolution. 1986;40:1092. https://doi.org/10.2307/2408767.

    Article  PubMed  Google Scholar 

  9. Stewart JR, Lister AM. Cryptic northern refugia and the origins of the modern biota. Trends Ecol Evol. 2001;16:608–13. https://doi.org/10.1016/s0169-5347(01)02338–2.

  10. Kurtén B. Pleistocene mammals of Europe. Routledge; 2017.

  11. Hewitt GM. Post-glacial re-colonization of European biota. Biol J Linn Soc. 1999;68:87–112. https://doi.org/10.1111/j.1095-8312.1999.tb01160.x.

    Article  Google Scholar 

  12. Meiri M, Lister AM, Higham TFG, Stewart JR, Straus LG, Obermaier H, et al. Late-glacial recolonization and phylogeography of European red deer (Cervus elaphus L.). Mol Ecol. 2013;22:4711–22. https://doi.org/10.1111/mec.12420.

    Article  PubMed  Google Scholar 

  13. Stewart KM, Pierce BM, Hundertmark KJ, Gasaway WC. Geographical variation in antler morphology of Alaskan moose: putative effects of habitat and genetics. Alces. 2002;38:155–65.

    Google Scholar 

  14. Boeskorov GG. A review of the systematics of Pliocene and Pleistocene moose, part 1. Cranium. 2005;22:26–55.

    Google Scholar 

  15. Breda M, Marchetti M. Systematical and biochronological review of Plio-Pleistocene Alceini (Cervidae; Mammalia) from Eurasia. Quat Sci Rev. 2005;24:775–805. https://doi.org/10.1016/j.quascirev.2004.05.005.

    Article  Google Scholar 

  16. Lister AM. Evolution of mammoths and moose: the Holarctic perspective. Morphol Change Quatern Mammals North Am. 1993:178–204. https://doi.org/10.1017/cbo9780511565052.009.

  17. Niedziałkowska M, Hundertmark KJ, Jędrzejewska B, Sidorovich VE, Zalewska H, Veeroja R, et al. The contemporary genetic pattern of European moose is shaped by postglacial recolonization, bottlenecks, and the geographical barrier of the Baltic Sea. Biol J Linn Soc. 2016;117:879–94. https://doi.org/10.1111/bij.12713.

    Article  Google Scholar 

  18. Niedziałkowska M. Phylogeography of European moose (Alces alces) based on contemporary mtDNA data and archaeological records. Mamm Biol. 2017;84:35–43. https://doi.org/10.1016/j.mambio.2017.01.004.

    Article  Google Scholar 

  19. Meiri M, Lister A, Kosintsev P, Zazula GD, Barnes I. Population dynamics and range shifts of moose (Alces alces) during the Late Quaternary. J Biogeography. 2020;June:1–12.

  20. Niedziałkowska M, Hundertmark KJ, Jędrzejewska B, Niedziałkowski K, Sidorovich VE, Górny M, et al. Spatial structure in European moose (Alces alces): genetic data reveal a complex population history. J Biogeogr. 2014;41:2173–84. https://doi.org/10.1111/jbi.12362.

    Article  Google Scholar 

  21. Wennerström L, Ryman N, Tison J-L, Hasslow A, Dalén L, Laikre L. Genetic landscape with sharp discontinuities shaped by complex demographic history in moose (Alces alces). J Mammal. 2016; x:gyv146.

  22. Björck S. A review of the history of the Baltic Sea, 13.0-8.0 ka BP. Quat Int. 1995;27:19–40. https://doi.org/10.1016/1040-6182(94)00057-c.

    Article  Google Scholar 

  23. Guthrie RD. New carbon dates link climatic change with human colonization and Pleistocene extinctions. Nature. 2006;441:207–9.

    Article  CAS  PubMed  Google Scholar 

  24. Udina IG, Danilkin AA, Boeskorov GG. Genetic diversity of moose (Alces alces L.) in Eurasia. Genetika. 2002;38:1125–32.

    CAS  PubMed  Google Scholar 

  25. Guthrie RD. Mammalian evolution in response to the Pleistocene-Holocene transition and the break-up of the mammoth steppe: two case studies. Acta Zool Cracov. 1995;38:139–54.

    Google Scholar 

  26. Elias SA, Short SK, Hans Nelson C, Birks HH. Life and times of the Bering land bridge. Nature. 1996;382:60–3. https://doi.org/10.1038/382060a0.

    Article  CAS  Google Scholar 

  27. DeCesare NJ, Weckworth BV, Pilgrim KL, Walker ABD, Bergman EJ, Colson KE, et al. Phylogeography of moose in western North America. J Mammal. 2020;101:10–23. https://doi.org/10.1093/jmammal/gyz163.

    Article  Google Scholar 

  28. Hundertmark KJ, Shields GF, Udina IG, Bowyer RT, Danilkin AA, Schwartz CC. Mitochondrial phylogeography of moose (Alces alces): late pleistocene divergence and population expansion. Mol Phylogenet Evol. 2002;22:375–87.

    Article  CAS  PubMed  Google Scholar 

  29. Świsłocka M, Matosiuk M, Ratkiewicz M, Borkowska A, Czajkowska M, Mackiewicz P. Phylogeny and diversity of moose (Alces alces, Cervidae, Mammalia) revealed by complete mitochondrial genomes. Hystrix It J Mamm. 2020;31. https://doi.org/10.4404/hystrix-00252-2019.

  30. R Development Core Team. R: A Language and Environment for Statistical Computing. R Found Stat Comput Vienna Austria, 0:ISBN 3–900051–07-0. 2016.

  31. Chen L, Qiu Q, Jiang Y, Wang K, Lin Z, Li Z, et al. Large-scale ruminant genome sequencing provides insights into their evolution and distinct traits. Science. 2019;364. https://doi.org/10.1126/science.aav6202.

  32. Gaillard JM. Are moose only a large deer?: some life history considerations. Alces. 2007;43:1–12.

    Google Scholar 

  33. Edenius L, Bergman M, Ericsson G, Danell K. The role of moose as a disturbance factor in managed boreal forests. Silva Fennica. 2002:36. https://doi.org/10.14214/sf.550.

  34. Pastor J, Dewey B, Naiman RJ, McInnes PF, Cohen Y. Moose browsing and soil fertility in the boreal forests of isle Royale National Park. Ecology. 1993;74:467–80. https://doi.org/10.2307/1939308.

    Article  Google Scholar 

  35. Mech LD, David Mech L, Fieberg J, Barber-Meyer S. An historical overview and update of wolf-moose interactions in northeastern Minnesota. Wildl Soc Bull. 2018;42:40–7. https://doi.org/10.1002/wsb.844.

    Article  Google Scholar 

  36. Schwartz CC, Swenson JE, Miller SD. Large carnivores, moose, and humans: a changing paradigm of predator management in the 21st century. Alces. 2003;39:41–63.

    Google Scholar 

  37. Timmermann HR, Rodgers AR. Others. Moose: competing and complementary values. Alces. 2005;41:85–120.

    Google Scholar 

  38. Mattsson L. Moose management and the economic value of hunting: towards bioeconomic analysis. Scand J For Res. 2008;5:575–81. https://doi.org/10.1080/02827589009382640.

    Article  Google Scholar 

  39. Laikre L, Larsson LC, Palmé A, Charlier J, Josefsson M, Ryman N. Potentials for monitoring gene level biodiversity: using Sweden as an example. Biodivers Conserv. 2008;17:893–910. https://doi.org/10.1007/s10531-008-9335-2.

    Article  Google Scholar 

  40. Lister AM. The impact of quaternary ice ages on mammalian evolution. Philos Trans R Soc Lond Ser B Biol Sci. 2004;359:221–41.

    Article  Google Scholar 

  41. Markova AK, Smirnov NG, Kozharinov AV, Kazantseva NE, Simakova AN, Kitaev LM. Late Pleistocene distribution and diversity of mammals in northern Eurasia. Paleontol I Evolucia. 1995;28–29:5–143.

    Google Scholar 

  42. Tarasov PE, Volkova VS, Webb T III, Guiot J, Andreev AA, Bezusko LG, et al. Last glacial maximum biomes reconstructed from pollen and plant macrofossil data from northern Eurasia. J Biogeogr. 2000;27:609–20.

    Article  Google Scholar 

  43. Hundertmark KJ, Bowyer RT. Genetics, evolution, and phylogeography of moose. Alces. 2004;40:103–23.

    Google Scholar 

  44. Calvignac S, Hughes S, Tougard C, Michaux J, Thevenot M, Philippe M, et al. Ancient DNA evidence for the loss of a highly divergent brown bear clade during historical times. Mol Ecol. 2008;17:1962–70.

    Article  CAS  PubMed  Google Scholar 

  45. Cronin MA. Intraspecific variation in mitochondrial DNA of north American Cervids. J Mammal. 1992;73:70–82. https://doi.org/10.2307/1381867.

    Article  Google Scholar 

  46. Goebel T, Waters MR, O’Rourke DH. The late Pleistocene dispersal of modern humans in the Americas. Science. 2008;319:1497–502.

    Article  CAS  PubMed  Google Scholar 

  47. Kitchen A, Miyamoto MM, Mulligan CJ. A three-stage colonization model for the peopling of the Americas. PLoS One. 2008;3:e1596.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Mulligan CJ, Kitchen A, Miyamoto MM. Updated three-stage model for the peopling of the Americas. PLoS One. 2008;3:e3199.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Hundertmark KJ, Bowyer TR, Shields GF, Schwartz CC. Mitochondrial Phylogeography Of Moose (Alces Alces) In North America. J Mammal. 2003;84:718–28. https://doi.org/10.1644/1545-1542(2003)084<0718:MPOMAA>2.0.CO;2.

  50. Li H, Durbin R. Inference of human population history from individual whole-genome sequences. Nature. 2011;475:493–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Mysterud A, Solberg EJ, Yoccoz NG. Ageing and reproductive effort in male moose under variable levels of Intrasexual competition. J Anim Ecol. 2005;74:742–54.

    Article  Google Scholar 

  52. Bowyer RT, Rachlow JL, Stewart KM, Van Ballenberghe V. Vocalizations by Alaskan moose: female incitation of male aggression. Behav Ecol Sociobiol. 2011;65:2251–60. https://doi.org/10.1007/s00265-011-1234-y.

  53. Thompson EA. Identity by descent: variation in meiosis, across genomes, and in populations. Genetics. 2013;194:301–26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Pemberton TJ, Absher D, Feldman MW, Myers RM, Rosenberg NA, Li JZ. Genomic patterns of homozygosity in worldwide human populations. Am J Hum Genet. 2012;91:275–92.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Shafer ABA, Cullingham CI, Côté SD, Coltman DW. Of glaciers and refugia: a decade of study sheds new light on the phylogeography of northwestern North America. Mol Ecol. 2010;19:4589–621.

    Article  PubMed  Google Scholar 

  56. DeCesare NJ, Smucker TD, Garrott RA, Gude JA. Moose status and management in Montana. Alces: A Journal Devoted to the Biology and Management of Moose. 2014;50:35-51.

  57. Kangas V-M, Kvist L, Kholodova M, Nygrén T, Danilov P, Panchenko D, et al. Evidence of post-glacial secondary contact and subsequent anthropogenic influence on the genetic composition of Fennoscandian moose (Alces alces). J Biogeogr. 2015;42:2197–208.

    Article  Google Scholar 

  58. Haanes H, Røed KH, Solberg EJ, Herfindal I, Sæther B-E. Genetic discontinuities in a continuously distributed and highly mobile ungulate, the Norwegian moose. Conserv Genet. 2011;12:1131–43. https://doi.org/10.1007/s10592-011-0214-0.

    Article  Google Scholar 

  59. Kangas V-M, Kvist L, Laaksonen S, Nygrén T, Aspi J. Present genetic structure revealed by microsatellites reflects recent history of the Finnish moose (Alces alces). Eur J Wildl Res. 2013;59:613–27. https://doi.org/10.1007/s10344-013-0712-0.

    Article  Google Scholar 

  60. Schönfeld F. Presence of moose (Alces alces) in southeastern Germany. Eur J Wildl Res. 2009;55:449–53. https://doi.org/10.1007/s10344-009-0272-5.

    Article  Google Scholar 

  61. Kromer B, Lindauer S, Synal H-A, Wacker L. MAMS – a new AMS facility at the Curt-Engelhorn-Centre for Archaeometry, Mannheim, Germany. Nucl Instrum Methods Phys Res, Sect B. 2013;294:11–3. https://doi.org/10.1016/j.nimb.2012.01.015.

    Article  CAS  Google Scholar 

  62. Reimer PJ, Bard E, Bayliss A, Warren Beck J, Blackwell PG, Ramsey CB, et al. IntCal13 and Marine13 radiocarbon age calibration curves 0–50,000 years cal BP. Radiocarbon. 2013;55:1869–87.

    Article  CAS  Google Scholar 

  63. Ramsey CB. Radiocarbon calibration and analysis of stratigraphy: the OxCal program. Radiocarbon. 1995;37:425–30.

    Article  Google Scholar 

  64. Yang DY, Eng B, Waye JS, Dudar JC, Saunders SR. Improved DNA extraction from ancient bones using silica-based spin columns. Am J Phys Anthropol. 1998;105:539–43.

    Article  CAS  PubMed  Google Scholar 

  65. Gamba C, Hanghøj K, Gaunitz C, Alfarhan AH, Alquraishi SA, Al-Rasheid KAS, et al. Comparing the performance of three ancient DNA extraction methods for high-throughput sequencing. Mol Ecol Resour. 2016;16:459–69.

    Article  CAS  PubMed  Google Scholar 

  66. Briggs AW, Stenzel U, Meyer M, Krause J, Kircher M, Pääbo S. Removal of deaminated cytosines and detection of in vivo methylation in ancient DNA. Nucleic Acids Res. 2010;38:e87.

    Article  CAS  PubMed  Google Scholar 

  67. Meyer M, Kircher M. Illumina Sequencing Library Preparation for Highly Multiplexed Target Capture and Sequencing. Cold Spring Harbor Protocols 2010.6 (2010): pdb-prot5448. https://doi.org/10.1101/pdb.prot5448.

  68. Maricic T, Whitten M, Pääbo S. Multiplexed DNA sequence capture of mitochondrial genomes using PCR products. PLoS One. 2010;5:e14004.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Heino MT, Askeyev IV, Shaymuratova DN, Askeyev OV, Askeyev AO, van der Valk T, et al. 4000-year-old reindeer mitogenomes from the Volga-Kama region reveal continuity among the forest reindeer in northeastern part of European Russia. In: Chapter 5. Archaeоlogy of the eurasian steppes no 4 2019, stone age and chalcolithic. Kazan: Tatarstan Academy of Sciences; 2019. p. 179–90.

    Google Scholar 

  70. Dabney J, Knapp M, Glocke I, Gansauge M-T, Weihmann A, Nickel B, et al. Complete mitochondrial genome sequence of a middle Pleistocene cave bear reconstructed from ultrashort DNA fragments. Proc Natl Acad Sci U S A. 2013;110:15758–63.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Gansauge M-T, Gerber T, Glocke I, Korlevic P, Lippik L, Nagel S, et al. Single-stranded DNA library preparation from highly degraded DNA using T4 DNA ligase. Nucleic Acids Res. 2017;45:e79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  72. Meyer M, Kircher M, Gansauge M-T, Li H, Racimo F, Mallick S, et al. A high-coverage genome sequence from an archaic Denisovan individual. Science. 2012;338:222–6. https://doi.org/10.1126/science.1224344.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Gansauge M-T, Meyer M. Single-stranded DNA library preparation for the sequencing of ancient or damaged DNA. Nat Protoc. 2013;8:737–48. https://doi.org/10.1038/nprot.2013.038.

    Article  CAS  PubMed  Google Scholar 

  74. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  75. Gnerre S, Maccallum I, Przybylski D, Ribeiro FJ, Burton JN, Walker BJ, et al. High-quality draft assemblies of mammalian genomes from massively parallel sequence data. Proc Natl Acad Sci U S A. 2011;108:1513–8.

    Article  CAS  PubMed  Google Scholar 

  76. Simpson JT, Wong K, Jackman SD, Schein JE, Jones SJM, Birol I. ABySS: a parallel assembler for short read sequence data. Genome Res. 2009;19:1117–23.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Luo R, Liu B, Xie Y, Li Z, Huang W, Yuan J, et al. SOAPdenovo2: an empirically improved memory-efficient short-read de novo assembler. Gigascience. 2012;1:18.

    Article  PubMed  PubMed Central  Google Scholar 

  78. Gurevich A, Saveliev V, Vyahhi N, Tesler G. QUAST: quality assessment tool for genome assemblies. Bioinformatics. 2013;29:1072–5.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31:3210–2.

    Article  CAS  PubMed  Google Scholar 

  80. Malde K, Skern R, Glover KA. Using sequencing coverage statistics to identify sex chromosomes in minke whales. arXiv [q-bio.GN]. 2019. http://arxiv.org/abs/1902.06654.

  81. Li H, Durbin R. Fast and accurate long-read alignment with burrows-wheeler transform. Bioinformatics. 2010;26:589–95.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  82. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. Bioinformatics. 2009;25:2078–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  83. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297–303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  84. Okonechnikov K, Conesa A, García-Alcalde F. Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data. Bioinformatics. 2016;32:292–4.

    CAS  PubMed  Google Scholar 

  85. Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011;27:2987–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  86. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  87. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  88. Pečnerová P, Palkopoulou E, Wheat CW, Skoglund P, Vartanyan S, Tikhonov A, et al. Mitogenome evolution in the last surviving woolly mammoth population reveals neutral and functional consequences of small population size. Evol Lett. 2017;1:292–303.

    Article  PubMed  PubMed Central  Google Scholar 

  89. Korneliussen TS, Albrechtsen A, Nielsen R. ANGSD: analysis of next generation sequencing data. BMC Bioinform. 2014;15:356.

    Article  Google Scholar 

  90. Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30:772–80.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  91. Zheng X, Levine D, Shen J, Gogarten SM, Laurie C, Weir BS. A high-performance computing toolset for relatedness and principal component analysis of SNP data. Bioinformatics. 2012;28:3326–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  92. Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu C-H, Xie D, et al. BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol. 2014;10:e1003537.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  93. Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9:772.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  94. Zurano JP, Magalhães FM, Asato AE, Silva G, Bidau CJ, Mesquita DO, et al. Cetartiodactyla: updating a time-calibrated molecular phylogeny. Mol Phylogenet Evol. 2019;133:256–62.

    Article  PubMed  Google Scholar 

  95. Rambaut A, Drummond AJ, Xie D, Baele G, Suchard MA. Posterior summarization in Bayesian Phylogenetics using tracer 1.7. Syst Biol. 2018;67:901–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  96. Rambaut A, Drummond AJ. TreeAnnotator v1. 7.0. Available as part of the BEAST package at http://beast.bio.ed.ac.uk. 2013.

  97. Haubold B, Pfaffelhuber P, Lynch M. mlRho - a program for estimating the population mutation and recombination rates from shotgun-sequenced diploid genomes. Mol Ecol. 2010;19(Suppl 1):277–84.

    Article  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Bodil Cronholm for laboratory assistance. We also thank Prof. Dr. Eberhard “Dino” Frey (State Museum of Natural History in Karlsruhe) and Frank Menger (private collection in Groß-Rohrheim) for providing access to four ancient moose specimens (Meng37, Meng43, Meng44, Qh37).  We acknowledge support from the Uppsala Multidisciplinary Centre for Advanced Computational Science for assistance with massively parallel sequencing and access to the UPPMAX computational infrastructure. We also thank the Wallenberg Advanced Bioinformatics Infrastructure (WABI; Science for Life Laboratory, Stockholm) for Long-Term Bioinformatics support from the National Bioinformatics Infrastructure Sweden (NBIS). Part of the work was carried out with the support of the Centre for Material Analysis, University of Oulu, Finland.

Authors’s contributions

N.D. and D.D.d.M. conceived the study. A.H., F.A., and M.T.H performed the laboratory work. L.D., K.L. and R.A.O. generated the de-novo reference genome. M.H., M.T.H., J.A., I.V.A, O.V.A., D.N.S. and A.O.A. provided ancient samples. N.R. and L.L. provided the modern sample used for de-novo genome assembly. M.T.H, F.A., D.D., R.F., S.L., W.R., M.H. performed the radiocarbon dating of ancient samples. N.D., D.D.d.M., T.v.d.V. and F.A. analysed the data. N.D. wrote the main body of the manuscript while all authors contributed to writing the manuscript. All authors have read and approved the manuscript.

Funding

M.T.H. acknowledges funding from Emil Aaltonen Foundation. D.D.d.M. and N.D. are funded by the Carl Tryggers Foundation (CTS 17:109 and 19:257). L. D acknowledges support from FORMAS (grant 2018–01640) and the Swedish Research Council (grant 2017–04647). We also acknowledge support from the Uppsala Multidisciplinary Centre for Advanced Computational Science for assistance with massively parallel sequencing and access to the UPPMAX computational infrastructure. Parts of this study (F.A., D.D., R.F., S.L., W.R., M.H.) were funded by the Klaus Tschira Stiftung Heidelberg within the project “Eiszeitfenster Oberrheingraben” at the Reiss-Engelhorn-Museen Mannheim. Open Access funding provided by Stockholm University.

Author information

Authors and Affiliations

Authors

Corresponding authors

Correspondence to Nicolas Dussex or David Díez-del-Molino.

Ethics declarations

Ethics approval and consent to participate

The moose sample used for the assembly generation was obtained from a frozen tissue bank collection from 1980 maintained by L.L. and N.R. at Stockholm University. No Ethics approval was required.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Additional file 2:

Figure S1. Identification of scaffolds linked to chromosome X using sequencing coverage and scaffold length. Green dots correspond to scaffolds assigned as linked to chromosome X. The red lines indicate the median coverage for all scaffolds. The blue lines represent half the median coverage which corresponds to the coverage of scaffolds linked to chromosome X. Figure S2. Past demography for moose (Alces alces) using PSMC. Thin lines represent 100 bootstrap runs. The x-axis corresponds to time before present in years on a log scale, assuming an estimated substitution rate of 0.7 × 10− 8 substitutions/site/generation [30] and a generation time of 7 years [31]. The y-axis corresponds to the effective population size Ne. Figure S3. Past demography for moose (Alces alces) using a Bayesian Skyline Plot (BSP). Demographic reconstruction was inferred in BEAST using 14 European 16,693 bp mitogenomes. Timing of events was estimated assuming a mean rate of 9 × 10− 9 substitutions/site/year based on Zurano et al. [93] and a standard deviation of 0.01. The x axis is in calendar years before present and y axis represents changes in effective population size (shown as the product of Nef and generation time T). The black line corresponds to the median estimate and the blue lines show the 95% highest posterior density intervals. Figure S4. Distribution of runs of homozygosity (ROH) in moose (Alces alces). ROH ≥500 kb are shown.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Dussex, N., Alberti, F., Heino, M.T. et al. Moose genomes reveal past glacial demography and the origin of modern lineages. BMC Genomics 21, 854 (2020). https://doi.org/10.1186/s12864-020-07208-3

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-020-07208-3

Keywords