Abstract
Bird-mediated seed dispersal is crucial for the regeneration and viability of ecosystems, often resulting in complex mutualistic species networks. Yet, how this mutualism drives the evolution of seed dispersing birds is still poorly understood. In the present study we combine whole genome re-sequencing analyses and morphometric data to assess the evolutionary processes that shaped the diversification of the Eurasian nutcracker (Nucifraga), a seed disperser known for its mutualism with pines (Pinus). Our results show that the divergence and phylogeographic patterns of nutcrackers resemble those of other non-mutualistic passerine birds and suggest that their early diversification was shaped by similar biogeographic and climatic processes. The limited variation in foraging traits indicates that local adaptation to pines likely played a minor role. Our study shows that close mutualistic relationships between bird and plant species might not necessarily act as a primary driver of evolution and diversification in resource-specialized birds.
Similar content being viewed by others
Introduction
Seed dispersal by birds is an essential ecological process that contributes to the regeneration and viability of plant communities1,2,3. The resulting mutualistic relationship between seed-dispersing birds and their food plants potentially leads to resource specialization and can subsequently act as a driving force in speciation2,4,5,6,7. Nutcrackers (genus Nucifraga) exemplify a case of mutualism and are known for their preference for caching pine seeds. Seeds that are not retrieved have the potential to germinate, thus making nutcrackers a crucial vector in pine seed dispersal8,9,10. Although nutcrackers have been reported to collect and cache the seeds of spruce (Picea) and hazel (Corylus)11,12 as well, they primarily prefer the seeds of white pines (Pinus, subgenus Strobus)13,14. This pine subgenus is characterized by wingless, nutrient-rich nuts and relies on animals to disperse their seeds as the pine cones of the majority of these species do not open at maturity13,15,16,17. The mutualistic relationship between nutcrackers and white pines has been extensively studied between the Nearctic Clark’s nutcracker (Nucifraga columbiana) and North American pines13,18,19 and has similarly been described for Eurasian nutcrackers and European and Asian pines10,20,21,22,23,24,25. Furthermore, white pines are distributed across a large geographic range26 and have deep divergences27,28. Yet their role in the diversification of nutcrackers remains unexplored.
Nutcrackers fulfil an essential role in the maintenance and viability of pine ecosystems by actively dispersing their seeds. Yet, despite their ecological importance, little is known about the evolution of nutcrackers and the processes that shaped their diversification. The variation in beak morphology previously reported between subspecies indicates that nutcrackers adapted to different pine species throughout their range and supports the hypothesis that their evolution was driven by its mutualism with pines8,21,29,30. Subsequently, beak morphology has been a primary trait on which subspecies have been distinguished, such as the “slender-billed” (N. c. macrorhynchos) and “thick-billed nutcrackers” (N. c. caryocatactes)29,30. However, the taxonomic classification of the Eurasian nutcrackers has been plagued by inconsistency and until now has solely been based on phenotypic descriptions between (sub)species, acknowledging either one (N. caryocatactes29,31), two (N. caryocatactes and N. multipunctata32,33,34,35) or three (N. caryocatactes, N. multipunctata, N. hemispila36) Eurasian nutcracker species. In the current paper, we follow the most recent classification of Madge et al. and Gill et al.34,35 which implies two species (N. caryocatactes and N. multipunctata) with eight defined subspecies in N. caryocatactes. Additionally, Madge et al.35 recognize two distinct groups within N. caryocatactes (northern and southern group), without taxonomic rank.
Whether the diversification of nutcrackers was indeed driven by their link with pines remains unclear, and understanding their evolution requires a combined effort of distribution, morphometric and genomic data. So far, a genomic assessment of the nutcracker species complex remains absent and the few studies that analysed the diversification of nutcrackers were based on mitochondrial markers limited to populations in the Northern Palearctic region37,38. Hence, a comprehensive genome-wide sampling covering the complete nutcracker range is needed to gain a deeper understanding of the driving forces behind the diversification of Eurasian nutcrackers.
In the current study, we combined both whole-genome re-sequencing and morphometric data, covering the complete Eurasian nutcracker distribution range, to examine what processes shaped the evolution of the nutcracker species complex. We show that beak morphology only marginally contributed to the overall phenotypic variation between (sub)species and that the phylogeographic pattern of Eurasian nutcrackers resembles those of other Sino-Himalayan and Palearctic non-mutualistic passerine birds. We therefore conclude that the early diversification of nutcrackers was most likely driven by environmental and climatic forces, instead of a strong mutualism with pines.
Results
Genetic, phylogenomic and morphometric comparison between the Eurasian nutcracker species
Principal component analysis (PCA) based on 110,979 linkage-pruned single nucleotide polymorphisms (SNPs) identified three distinct genetic clusters among the 31 re-sequenced Eurasian nutcracker samples (Fig. 1a, Fig. 1b and Supplementary Fig. 1). The first principal component explained 35.7% of the variation and separated the northern N. caryocatactes group from both N. multipunctata and the southern N. caryocatactes group, whereas the second principal component explained 13.4% of the variation and separated N. multipunctata from the other two groups. Admixture analyses based on the same SNP dataset with ancestry components (K) ranging from 2 to 4 supported this pattern and identified K = 3 as the optimal number of clusters (Fig. 1c and Supplementary Fig. 2).
a Species range according to Madge et al.35 and Gill et al.34 with sampling locations (black dots) of N. caryocatactes (northern) (N = 17), N. caryocatactes (southern) (N = 8) and N. multipunctata (N = 3). Subspecies distributions are depicted in different colours according to the legend. b Principal component analysis based on 110,979 SNPs separated the data in three distinct clusters. c Population structuring based on admixture analysis for K ranging from 2 to 4. Cross-validation indicates K = 3 as the optimal number of clusters. d Species tree under the Multiple Species Coalescence model constructed from 192 nuclear genomic alignments of 1 Mbp. e Linear discriminant analysis of phenotypic variation among Eurasian nutcrackers based on an independent morphometric dataset of 90 museum specimens. Bird images by Jon Baldur Hlidberg, courtesy of the authors.
Maximum likelihood (ML) phylogenetic inference was performed based on whole mitochondrial (mt) genome alignments, revealing three distinct evolutionary lineages corresponding to a northern N. carycocatactes clade, southern N. caryocatactes clade and N. multipunctata clade. However, the internal node received a low support value and therefore reflects a poorly resolved tree (Supplementary Fig. 3). Maximum likelihood phylogenetic inference and multispecies coalescent (MSC) species tree reconstruction of the nuclear genome based on 192 non-overlapping 1 Million base pairs (Mbp) alignments similarly revealed three distinct genetic lineages and supported a topology in which N. caryocatactes (northern) is a sister group to N. caryocatactes (southern) and N. multipunctata (Fig. 1d, Supplementary Figs. 4 and 5). However, with a normalized quartet score of 71.65%, the MSC species tree indicates conflicting phylogenetic signals within the nuclear genome.
To assess potential reticulation as a cause of the discordance between the mitochondrial and nuclear data and low node support value, we subsequently constructed a phylogenetic network based on the 192 alignments of 1 Mbp, which consistently revealed genomic conflict for the deepest divergences between nutcrackers for thresholds ranging between 0.01 and 0.30 thus depicting reticulate evolution for alternative branches that were supported by up to 30% of the 192 alignments (Fig. 2). One third (32.2%) of the alignments supported the split of N. multipunctata from N. caryocatactes (northern) and N. caryocatactes (southern), whereas 67.7% of the 1 Mbp alignments supported the split of N. caryocatactes (northern) from N. caryocatactes (southern) and N. multipunctata.
Phylogenetic network of the Eurasian nutcracker species complex based on 192 alignments of 1 Mbp of the nuclear genome. The numbers at the branches indicate the percentage of genomic alignments supporting the split. Some branch lengths are shortened for visualization purposes. Bird images by Jon Baldur Hlidberg, courtesy of the authors.
Divergence time estimates of the Eurasian nutcracker species complex were based on the mtDNA alignment of 59 individuals covering 31 taxa (Supplementary Fig. 6 and Supplementary Table 1). The nutcrackers and crows diverged in the late Miocene at ~8.00 (5.74−10.20) million years ago (Ma), whereas the Nearctic Clark’s nutcracker (N. columbiana) and Eurasian nutcrackers diverged during the late Miocene/early Pliocene at around ~4.72 (3.2–6.23) Ma. Within the Eurasian nutcracker species complex, the deepest split occurred during the late Pliocene/early Pleistocene at ~2.32 (1.59–3.11) Ma and the youngest split at ~1.87 (1.17–2.60) Ma (N. caryocatactes northern and N. multipunctata). The within-species divergence occurred only recently at ~340–50 kilo years ago (ka) (Supplementary Fig. 7).
Phenotypic variation was examined based on an independent morphometric dataset of 90 museum specimens. First, a Principal Coordinate Analysis (PCoA) was performed based on 12 phenotypic traits to assess the morphometric variation within the data without any prior group assumptions (Supplementary Figs. 8 and 9). Results indicate two clusters, separating the northern N. caryocatactes group from the southern N. caryocatactes group and N. multipunctata. This variation was primarily driven by wing length and white extent on the tail feathers. In contrast, linear discriminant analysis (LDA) based on a priori defined groups revealed the presence of three distinct clusters (Fig. 1e, Supplementary Fig. 10). The first linear discriminant accounted for 76.52% of the phenotypic variation and was similar to the PCoA primarily characterized by variation in white extent on the tail feathers, separating the northern N. caryocatactes group from the southern N. caryocatactes group and N. multipunctata (Supplementary Table 2). However, the second linear discriminant explained 8.63% of the variation, separating N. multipunctata from both the northern and southern N. caryocatactes groups, and was characterized by variation in tail length. Comparison of beak traits revealed only a weak yet significant difference in bill morphology (Welch’s t-test p ≤ 0.05), with shorter and higher bills in the southern N. caryocatactes nutcrackers compared to the northern N. caryocatactes nutcrackers (Supplementary Fig. 11).
Genetic, phylogenetic and morphometric comparison between Eurasian nutcracker subspecies
We performed both PCA and admixture analyses for the northern and southern N. caryocatactes groups separately to examine genetic differentiation at the subspecies level (Figs. 3 and 4). PCA based on 142,799 SNPs identified three clusters within the northern N. caryocatactes group (Fig. 3b). The first principal component explained 11.8% of the variation and separated the European subspecies (N. c. caryocatactes) from the other subspecies. The second principal component explained 8.6% of the variation and separated the Central Asian subspecies (N. c. rothschildi). No clear genetic differentiation was found between the Siberian (N. c. macrorhynchos) and Japanese (N. c. japonica) subspecies. The admixture analyses reflected a similar pattern as the PCA for K ranging from 2 to 4. However, cross-validation tests suggest K = 2 as the most likely number of genetic clusters, only assigning N. c. caryocatactes to its own genetic group (Supplementary Figs. 2 and 12). Both the PCA and admixture analyses identified one individual as a likely F1-hybrid between N. c. caryocatactes and N. c. macrorhynchos.
a Distribution of the four subspecies within the northern N. caryocatactes group. Two specimens sampled in the range of N. c. caryocatactes were taxonomically and genetically identified as N. c. macrorhynchos. b PCA based on 142,799 SNPs separated the northern N. caryocatactes nutcrackers in three clusters. The specimen located in the middle of N. c. caryocatactes and N. c. macrorhynchos/N. c. japonica reflects a potential F1-hybrid. c Linear discriminant analyses based on morphometric measurements of 63 northern N. caryocatactes museum specimens. d Reconstruction of historical effective population size through time based on PSMC. N. c. rothschildi was excluded due to limited genome coverage.
a Distribution of the four subspecies within the southern N. caryocatactes group. b Within-species genetic variation based on a PCA of 101,611 SNPs. c Linear discriminant analyses based on morphometric measurements of 23 southern N. caryocatactes museum specimens. N. c. interdicta had to be excluded due to lack of trait data. d Reconstruction of historical effective population size through time based on PSMC. N. c. hemispila was excluded from the analyses due to insufficient genome coverage.
PCoA of morphometric traits revealed no clear variation between groups (Supplementary Figs. 8 and 9) and a LDA based on a priori defined groups revealed only marginal phenotypic differentiation among subspecies, as indicated by the overlap in clusters (Fig. 3c). The main phenotypic traits that contributed to the linear discriminants consisted of wing length, white extent on the tail feathers and variation in tail length (Supplementary Table 3). Trait-specific analyses of beak morphology revealed significant variation between subspecies, with N. c. caryocatactes having wider bills and N. c. japonica shorter bills compared to the other northern N. caryocatactes subspecies (Supplementary Fig. 13).
Reconstruction of the demographic history using pairwise sequentially Markovian coalescent (PSMC) suggests that the Siberian (N. c. macrorhynchos) and Japanese (N. c. japonica) nutcrackers underwent a sharp increase in effective population size (Ne) around 50–100 ka ago, followed by a rapid decrease around 30–40 ka ago (Fig. 3d, Supplementary Figure 14). In contrast, the European population (N. c. caryocatactes) had a low and continuously decreasing Ne throughout the past 400,000 years. No demographic estimates were included for N. c. rothschildi due to insufficient genome coverage.
PCA based on 101,611 SNPs identified three clusters within the southern N. caryocatactes group. The first principal component explained 22.5% of the variation and separated the Taiwanese subspecies (N. c. owstoni) from the other subspecies (Fig. 4b). The second principle component explained 16.6% of the variation and separated the subspecies of North-West/Central Himalayas (N. c. hemispila) from the other subspecies. The two subspecies in South and North China showed strong genetic similarity (N. c. macella and N. c. interdicta). Admixture analyses with K ranging from 2 to 4 identified distinct genetic groups for each of the subspecies at K = 4 (Supplementary Fig. 12), though cross-validation tests suggest K = 2 as the most likely number of clusters (Supplementary Fig. 2).
Similar to the northern group, the PCoA of the morphometric traits revealed no clear variation between the southern groups (Supplementary Figs. 8 and 9). However, LDA identified N. c. macella and N. c. hemispila as distinct subspecies (Fig. 4c). The first linear discriminant explained 76.1% of the variance and was mainly driven by differences in bill length and the distal distance of the wings (distance between the tip of outermost primary and tip of the wing), whereas the second linear discriminant explained 23.9% of the variance and was primarily driven by variation in wing length, separating N. c. owstoni from the other two subspecies. However, none of the traits differed significantly between subspecies (Supplementary Table 4). Measurements for N. c. interdicta were incomplete and had to be excluded from this analysis. Separate analyses of beak morphology revealed no significant variation in beak traits between the southern N. caryocatactes subspecies (Supplementary Fig. 13).
Analyses of the historical effective population size revealed a continuous decline in Ne of the Taiwanese subspecies (N. c. owstoni) from 400 ka ago until present (Fig. 4d and Supplementary Fig. 14). Both the southern and northern Chinese subspecies (N. c. macella and N. c. interdicta) underwent a decrease in Ne around 400–100 ka ago. However, N. c. macella underwent a rapid increase and decrease in Ne around 100–60 ka ago, whereas N. c. interdicta continued to decrease in the same period. Both populations seem to have stabilized afterwards. No demographic estimates were included for N. c. hemispila due to insufficient coverage in both samples.
Discussion
In the current study, we used a whole-genome sequencing approach in combination with morphometric analyses to assess the evolutionary and phylogeographic history of the Eurasian nutcracker, a food caching corvid known for its mutualism with pines. The results presented here provide clear evidence that the Eurasian nutcracker species complex is characterized by three evolutionary lineages, corresponding to a widespread northern Palearctic clade (N. caryocatactes (northern)), a widespread Sino-Himalayan clade (N. caryocatactes (southern)) and a narrow-range endemic Western Himalayan clade (N. multipunctata). Our data supports three species-level taxa of Eurasian nutcrackers in accordance with previous classifications recognizing the southern group as its own species N. hemispila36, warranting a taxonomic revision. Divergence time estimates placed the deepest divergences of the Eurasian nutcrackers at 2.32 Ma, followed by a second split around 1.87 Ma (Supplementary Figs. 6 and 7). However, our analyses could not fully resolve the phylogenetic relationship between the three lineages, most likely due to the short time span between the two speciation events which resulted in high levels of incomplete lineage sorting39,40. Past introgressive gene flow may complicate the phylogenetic reconstruction of nutcrackers, but given the unresolved topology, the hypothesis is with current methods not possible to examine. The evolutionary history of the Eurasian nutcracker species complex therefore seems to be characterized by a rapid speciation event that occurred at the beginning of the Pleistocene. Rather than a bifurcating phylogeny, the phylogenetic history of the Eurasian nutcracker species complex is instead better depicted as a reticulate network (Fig. 2), a pattern that seems to become common for many avian species with the increasing number of whole-genome studies40,41,42,43.
The spatial diversification pattern of nutcrackers is rather different from those found in other corvids44,45,46. Many Palearctic corvid species do not have any close relatives in the (Sino-)Himalayas37,47, whereas several Sino-Himalayan corvid species represent an Indomalayan diversification without any closer affiliation to the Palearctic48. Instead, the divergence times and phylogeographic distribution of nutcrackers seem to closely match those of other montane passerines of the temperate forest belt in the Sino-Himalayas, particularly for those clades encompassing endemics of the Western Himalayan subalpine conifer forests. Of the randomly distributed 192 non-overlapping 1 Mbp alignments of the nuclear genome, 67.7% suggested a phylogenetic relationship characterized by close relatedness between the Western Himalayan and Sino-Himalayan nutcrackers which formed a sister group to the northern Palearctic nutcrackers (Figs. 1d, 2, Supplementary Figs. 15 and 16), a phylogeographic pattern which is common among Sino-Himalayan passerine birds49,50,51,52,53 (Supplementary Fig. 17). The alternative phylogeographic pattern in which the northern Palearctic and Sino-Himalayan nutcrackers form a sister group to the Western-Himalayan nutcrackers was less frequently supported by the genomic analyses (Fig. 2, Supplementary Figs. 15 and 16). This phylogeographic pattern has only previously been found in leaf warblers (e.g. Phylloscopus subviridis)54 and therefore seems to be less common among Sino-Himalayan birds. Overall, the similarity in phylogeographic patterns and divergence times between the nutcrackers and Sino-Himalayan passerines, suggests that similar environmental and climatic forces have shaped the diversification of these species. Future studies utilizing ecological niche modelling would be important to better understand these phylogeographic patterns on a broader scale.
The northern Palearctic region is characterized by an East-West (E-W) phylogeographic pattern for passerine birds, caused by strong Pleistocene climatic changes and has been observed in several corvid species37,45,46,52,55,56,57,58,59. Similarly, our genome-wide analyses showed a clear E-W disjunction between the European nutcrackers and the other northern Palearctic nutcrackers, in contrast to previous studies which showed no clear E-W divergence based on the mitochondrial control region37,38. Additionally, one hybrid individual identified in our dataset originated from Sweden and suggests ongoing gene flow between the European and Siberian nutcrackers, likely linked to frequent mass irruptive migrations of Siberian nutcrackers into Western Europe during times of food scarcity11,12,31,60. In Sweden, local colonies of Siberian nutcrackers (N. c. macrorhynchos) became established after the invasion event of 197761 which would explain the occurrence of admixed individuals in that area. Interestingly, both the Japanese (N. c. japonica) and Siberian nutcrackers (N. c. macrorhynchos) showed strong genetic similarity and overlap in their historical effective population size, despite their large distribution range, suggesting either a panmictic ancestral population in shared refugia37,62,63 or extensive gene flow and migration during the late Pleistocene via the existence of land bridges between Japan and continental Eurasia64.
Our analyses could similarly identify another clear E-W phylogeographic disjunction in the Sino-Himalayan nutcrackers, separating the North-West/Central Himalayan nutcrackers (N. c. hemispila) from those of South-West China, originating around 300–100 ka. The divergence seems to have co-occurred with cyclic contractions and expansions of conifer forests during the late Pleistocene65 and reflects a phylogeographic pattern found in other montane passerine birds of the Sino-Himalayas51,66. Furthermore, the Taiwanese subspecies (N. c. owstoni) was genetically distinct compared to the other Sino-Himalayan nutcrackers and together with the continuously decreasing effective population size over the last 300–400 ka, reflects a loss of genetic diversity typically found in island populations as a result of drift, inbreeding and genetic isolation67.
Despite previous reports of phenotypic variation in beak morphology between nutcracker subspecies, it is still unclear whether the mutualism between nutcrackers and pines resulted in co-evolutionary processes in which local adaptation to pine cone morphology was the driving factor of beak variation. For example, crossbills (Loxia) illustrate a typical co-evolving arms-race with lodgepole pines (Pinus contorta latifolia) where divergent selection on beak morphology contributes to reproductive isolation of populations resulting in different crossbill species6,68,69. Moreover, only limited genome-wide differentiation was found between crossbill ecotypes, present only at a small number of loci that arose over as little as 6000 years6,7. If nutcrackers represent a similar case of co-evolution, we would expect that phenotypic differentiation between subspecies was primarily driven by variation in beak morphology. Our analyses show that variation in beak morphology between nutcracker species only contributed marginally to the phenotypic variation between species. Phenotypic diversification of the three main Eurasian nutcracker lineages was primarily driven by variation in flight-related traits (e.g. tail length) as well as in the white extent on the tail feathers, thus suggesting selection pressures unrelated to food resources as a driver for phenotypic diversification.
At the subspecies level, the phenotypic variation of northern Palearctic nutcrackers was similarly driven by variation in flight-related traits and white extent on the tail feathers. No significant morphological variation was found between the Sino-Himalayan subspecies, though the low number of samples and lack of N. c. interdicta in the dataset requires more extensive sampling to fully understand the phenotypic variation within the Sino-Himalayan nutcrackers. Yet, these initial results seem to support the idea that variation in beak morphology only marginally contributed to the phenotypic differentiation between species and subspecies and that flight-related traits are driving diversification. An extended in-depth analysis focussing only on beak traits revealed significant variation between several northern Palearctic subspecies. For example, the “Thick-billed” nutcracker (N. c. caryocatactes) had significantly wider bills compared to the other subspecies, which has previously been linked to their preference for both hazelnuts and pine seeds29,30,70. Interestingly, the Japanese nutcrackers had a significantly shorter beak compared to the other northern Palearctic subspecies, despite the lack of genome-wide differentiation with Siberian nutcrackers (Fig. 3b and Supplementary Fig. 13). As was shown in the case of Loxia crossbills, as well as in Darwin finches71,72, diversification of ecotypes can occur in short timescales without genome-wide diversification. It is therefore possible that the Siberian and Japanese nutcrackers only recently became adapted to different food sources. Indeed, the Japanese nutcrackers have been linked with Japanese stone pine (Pinus pumila)23,24 while the Siberian nutcracker feeds preferentially on Siberian pine (Pinus sibirica)73, which typically has longer seeds compared to Japanese stone pines (10–14 mm vs 7–10 mm, respectively)17,74. Additionally, some nutcracker subspecies occur in regions not occupied by pine species and have most likely adapted to different food sources, such as Asian spruce (Picea schrenkiana) in the case of N. c. rothschildi (Supplementary Fig. 1 d). Combined with the known preference of N. c. caryocatactes for both pine and hazel seeds, nutcrackers therefore seem to be characterized by a more opportunistic feeding strategy and support the hypothesis that local adaptation to pines played only a minor role in the diversification of these birds. However, subsequent studies are needed to directly link the variation in nutcracker beaks to the pine cone and pine seed morphology. Furthermore, more extensive genomic sampling is needed in order to link beak traits to specific genomic regions71,72 to get a deeper understanding of the underlying genetic basis of beak variation in nutcrackers.
Although nutcrackers are characterized by a strong mutualism with pines, our results show that resource specialization and subsequent local adaptation resulted in only minor phenotypic differentiation between nutcracker lineages. The strongest overall differences between species lie in flight-related traits rather than beak characters. Furthermore, the divergence times and phylogeographic patterns of nutcrackers were highly congruent with those of other non-mutualistic passerine birds, suggesting that similar evolutionary and climatic forces shaped the evolution of these species. Divergence time estimates could link these to climatic fluctuations during the Pleistocene period. Our study therefore suggests that a mutualistic relationship between plants and animals does not necessarily need to coincide with strong selective pressures on speciation and can instead be characterized by a high degree of evolutionary flexibility where other processes such as climatic and environmental factors are more determinant in their diversification.
Materials and methods
De novo assembly
As a basis for the genomic analyses, a reference nutcracker genome was assembled using 10X Genomics linked read technology. A blood sample from a northern N. caryocatactes specimen (Voucher ID: NHMO-BI-40098/1-B Natural History Museum of Oslo) was sent to Beijing Genomic Institute (BGI) Hong Kong for DNA extraction and library preparation. Sequencing was subsequently conducted on a single Illumina HiSeq X Ten lane (2 × 151 bp), yielding a total of 88.9 Gbp of raw sequence data.
Quality control of raw read data was performed with FastQC v0.11.775 and subsequent de novo assembly of the genome was performed with the 10X genomics Supernova v.1.0.0 assembly pipeline76. The quality of the assembly was assessed with QUAST v. 6.4.377 and gene completeness was assessed with BUSCO v.3.0.278 using the Aves gene set (available at busco.ezlab.org). Synteny differences between the Nucifraga assembly and the collared flycatcher (Ficedula albicollis) assembly (acc. No. GCA_000738735.2) were analysed with JupiterPlot79 (Supplementary Fig. 18). For downstream analyses, scaffolds smaller than 100 Kbp were filtered out, with the exception of those that contained BUSCO genes. The final assembly had a total length of 1.14 Gbp, with a contig N50 of 217.12 Kbp and scaffold N50 of 12.61 Mbp (Supplementary Table 5). BUSCO assessment identified 93.4% complete, 4.0% fragmented and 2.6% of missing genes of the Aves set (n = 4915).
RepeatModeler v.1.0.1180 was used to identify de novo repeats in the nutcracker genome. The resulting de novo repeat database was merged with the custom annotated repeat library used for studying crows81. The custom repeat library was subsequently used to mask repetitive sequences in the nutcracker genome with RepeatMasker v.4.1.082. Repetitive regions comprised 9.23% of the genome (Supplementary Table 6).
Re-sequencing, mapping and genotype calling
In the current paper, we have assumed the classification of Madge et al.35 and Gill et al.34 that acknowledge two Eurasian nutcracker species (N. caryocatactes and N. multipunctata) with eight defined subspecies in N. caryocatactes. Additionally, Madge et al.35 recognize two distinct groups within N. caryocatactes (northern and southern group), though they do not officially acknowledge them as independent species. We therefore sampled a total of 31 individuals covering each species and including at least a minimum of two individuals per subspecies (Fig. 1a). Additionally, one individual from the closely related Clark’s nutcracker (Nucifraga columbiana) was sequenced to use as an outgroup for the phylogenomic analyses. Samples consisted of either blood, tissue or toe pads, originating from varying natural history collections (Supplementary Table 7). DNA from tissue and blood samples was extracted with the Macherey-Nagel NucleoSpin® Tissue kit according to the manufacturer’s protocol. DNA from museum samples was extracted from 3 × 3 mm toe pads83. Sequencing libraries were prepared with 300 bp insert-sizes and sequenced with the Illumina Hiseq X technology at BGI (HongKong) or at the National Genomics Institute in Stockholm (SciLife). Additionally, we retrieved publicly available short-read sequencing data for the closely related hooded crow (Corvus cornix) (Acc. SRR1266946 -SRR1266949) as an additional outgroup.
Reads were processed using a custom-designed, clean-up workflow (https://github.com/mozesblom). Cleaning steps included quality control with FastQC v0.11.7, deduplication with Super-Deduper v2.0 (https://github.com/dstreett/Super-Deduper), trimming with Trimmomatic v.0.38 set at SLIDINGWINDOW:4:20 and MINLEN:3084, merging of overlapping reads with PEAR v0.9.685 for reads that overlapped for ≥20 base pairs and removal of low complexity reads with a cut-of threshold of ≥0.5. The resulting reads were mapped to the de novo nutcracker genome with the default settings of BWA mem v.0.7.1786 (Supplementary Table 8), and duplicates were marked with Picard v.2.18.21 (http://broadinstitute.github.io/picard/). Subsequently, we used freebayes v1.3.187 to perform variant calling with the following settings:-report-monomorphic –min-mapping-quality 20 -C4 and -F 0.3. For the whole-genome alignments, indels were removed and consensus sequences were created from the VCF files using custom perl scripts (https://github.com/mobilegenome/phylogenomics/), masking out ambiguous sites with N’s and removing all positions with missing data.
Phylogenetic analyses
We reconstructed the mitochondrial genomes for all individuals by mapping the cleaned reads to the mitochondrial genome of the Clark’s nutcracker (Nucifraga columbiana) (KF509923) using BWA mem, followed by variant calling with freebayes87, using the same settings as described above for the nuclear genome. Bcftools’ consensus function was then used to call consensus sequences. Mitochondrial genomes were aligned with MAFFT v7.407 (settings L-INS-i) and we included previously published mitochondrial genomes of the Magpie (Pica pica—HQ915867) and hooded crow (Corvus cornix—CM002877). Subsequent Maximum Likelihood tree inference was performed with IQ-TREE v.1.6.11 using the optimal substitution model identified by ModelFinder88 implemented in IQ-TREE89.
Phylogenomic inference of the nuclear genome proceeded using a method similar to a non-overlapping sliding window approach90. First, consensus sequences were generated for each individual, followed by the removal of heterozygous sites and N’s for all individuals using bedtools v2.28.0. Sequences were subsequently aligned per scaffold with custom python scripts (https://github.com/mobilegenome/phylogenomics/). To identify the minimum nucleotide length of the sliding windows that carry sufficient phylogenetic information to reject alternative topologies, increasingly longer windows were analysed by using an approximately unbiased (AU) test91 as implemented in IQ-TREE v.1.6.11 (Supplementary Figs. 15 and 16). The AU test was performed with 10,000 replicates for window lengths of 50 Kbp to 2 Mbp in steps of 50 Kbp, with 100 randomly sampled alignments per window length. Of the three potential nutcracker topologies that were tested with the AU test, one could be significantly rejected when windows had a length of ≥1 Mbp. Increasing the window length further could not significantly accept or reject the other two potential topologies (Supplementary Fig. 16). Subsequently, the scaffold alignments were split into 1 Mbp non-overlapping windows and only those windows that contained sufficient phylogenetic content (>2000 informative sites) were retained, resulting in a total of 192 1 Mbp alignments (192 Mbp, representing 19.6% of the repeat-filtered genome).
Maximum Likelihood trees were inferred for each 1 Mbp alignment in IQ-TREE using 1,000 ultrafast bootstrap replicates and the implemented automatic ultrafast model selection. From the resulting trees, we constructed a multispecies coalescent species tree using ASTRAL-III v5.6.192. The trees were rooted with Clark’s nutcracker and hooded crow as outgroups. Lastly, a consensus network was generated based on the 192 inferred ML trees. The “Consensus Network” method was used as implemented in SplitsTree v.4.15.193 with a minimum proportion of trees supporting the splits ranging between 0.05 and 0.4.
Allele-based analyses
VCF files were filtered for biallelic positions and SNPs were removed when coverage was less than 8× or more then 2.5 times the mean coverage for that particular individual. To account for the bias resulting from mapping to an in-group reference genome and the biases introduced by combining both museum and fresh samples in our study, we applied a strict filtering approach in which at maximum one individual was allowed to have missing data at each site. The final number of SNPs within the total dataset was reduced from 11,602,883 to 837,481 SNPs after filtering. Two allele frequency-based analyses were implemented to examine the genetic variation within the data. First, a PCA was performed between and within species with Plink version 1.994. To account for linkage disequilibrium, sites within a 50-SNP stepping window with a correlation coefficient higher than 0.1 were omitted. Secondly, ADMIXTURE v1.3.095 was used to estimate the population structure within the data. The VCF-file was converted to plink’s PED format with VCFtools and admixture analyses were run for K = 2 to K = 4. ADMIXTURE’s cross-validation procedure was used to determine the optimal number of clusters.
Divergence time estimation
Divergence times were estimated using BEAST2 (version 2.7.0)96. The bmodeltest97 with a transition-transversion split was used to find a fitting substitution model. The mitochondrial protein-coding genes were extracted and aligned using MAFFT v7.407 (settings L-INS-i)98 and manually inspected. The ND6 gene was excluded from the dataset due to being encoded on the opposite strand. The final alignment of 10,851 nucleotides long included 59 individuals, covering 31 species (Supplementary Table 9). Five fossil-based calibration points were selected from the literature99,100 and used to estimate the divergence times (Supplementary Table 1). Zero offset values were based on the minimum fossil ages and the log mean and standard deviation was adjusted to achieve the optimal match of the normal prior distribution for the fossil age interval101. Additionally, we applied a normal prior to the root age to avoid the occurrence of implausibly old root ages. All calibrated nodes were set to be monophyletic. Four independent Markov chain Monte Carlo (MCMC) were run for 100,000,000 generations sampling every 10,000. Trace files were examined with Tracer v1.7.1102 to ensure chain convergence and appropriate effective sample sizes (ESS > 200). The tree log files from each iteration were combined with LogCombiner with a 50% burn-in and maximum clade credibility (MCC) tree was summarized in TreeAnnotator with a 20% burn-in. Both LogCombiner and TreeAnnotator are part of the BEAST2 package.
Population size history
The demographic history was assessed with PSMC v.0.6.5103. Consensus sequences in FASTQ format were generated for each individual that had a coverage of ≥17X104 using the ‘mpileup’ command in BCFtools v.1.9 (https://samtools.github.io/bcftools/) and the included script ‘vcfutils.pl’. Sites with a read depth below 8 or above twice the sample’s median depth and a minimum mapping quality below 20 were removed. PSMC was run for 25 iterations with the upper limit of the TMRCA set to -t 5, the initial h/q value to -r 1, and 34 atomic time intervals (4 + 30*2 + 4 + 6 + 10), as is recommended for avian genomes104. We performed 100 bootstrap replicates by randomly sampling with replacement 1 Mbp blocks from the consensus sequence. Results were scaled using a mutation rate (m) of 3.18 × 10−9 substitutions per generation81 and a generation time of 7 years105 (Supplementary Fig. 14).
Morphometric analyses
To assess the phenotypic differentiation among nutcrackers, we measured 12 morphological traits on nutcracker specimens from museum collections. We used morphological traits associated with nutcracker foraging and diet, flight ability, manoeuvrability, bipedal locomotion and plumage colour (detailed description of all measures in Supplementary Table 10). Traits were measured on a total number of 118 specimens from the Zoological Research Museum Alexander Koenig, Bonn, Germany, the Natural History Museum, Berlin, Germany, Senckenberg Natural History Collections, Dresden, Germany and the Senckenberg Forschungsinstitut und Naturmuseum, Frankfurt, Germany. First, phenotypic variation was assessed without a priori defined groups by implementing a PCoA as part of the ape package (v. 5.5)106 in R v. 3.5.2107. Subsequently, a LDA was performed to assess the morphometric variation between groups defined by our genomic analysis. Specimens with missing data were excluded from the analyses, resulting in a total sample size of N = 90. Similarly, a LDA was performed to assess morphometric variation at the subspecies level for both the northern N. caryocatactes (N = 63) and southern N. caryocatactes (N = 23) nutcrackers. LDA were performed with the ggord package108 in R v. 3.5.2107. Traits that contributed primarily to the linear discriminants were further examined. First, a Shapiro-Wilk test was performed to test whether the traits adhered to a normal distribution. Subsequently, a Welch’s t-test with Bonferroni correction was performed to compare means between groups. A similar approach was implemented to examine the variation in beak traits between groups.
Vegetation maps of the main food resources linked to nutcrackers26 were plotted with R108 (Supplementary Fig. 1 and Supplementary Table 11). In total the distribution of 26 plant species that have been linked to, or occur in the same geographical region as the nutcrackers were selected and plotted. This included two species of spruce (Picea), eight species of hazel (Corylus) as well as 16 species of pine (Pinus) that occur in the Palearctic region. All 11 species of white pines (Pinus subgenus Strobus Section Quinquefoliae subsection Strobus) occurring in the Palearctic region were selected.
Statistics and reproducibility
Information about the statistical analyses is given in the methods section, results section and the supplementary file.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.
Data availability
The data produced in the current study are available at the NCBI Sequence Read Archive under BioProject PRJNA682958 accession numbers SAMN17014695–SAMN17014724. The raw 10X genomics sequencing data used for the nutcracker reference assembly have been deposited at NCBI under accession number SAMN17050349. All supporting data are uploaded to Dryad, https://doi.org/10.5061/dryad.stqjq2c57109.
Code availability
All software versions and scripts needed to reproduce the results are specified in the Methods.
References
Howe, H. & Smallwood, J. Ecology of seed dispersal. Annu. Rev. Ecol. Syst. 13, 201–228 (1982).
Vander Wall, S. B. & Beck, M. J. A comparison of frugivory and scatter-hoarding seed-dispersal syndromes. Botanical Rev. 78, 10–31 (2012).
Eriksson, O. Evolution of angiosperm seed disperser mutualisms: the timing of origins and their consequences for coevolutionary interactions between angiosperms and frugivores. Biol. Rev. 91, 168–186 (2016).
Snow, D. W. Coevolution of birds and flowers*. Curtis’s Botanical Mag. 11, 198–206 (1994).
Grant, P. & Grant, B. 40 Years of Evolution: Darwin’s Finches on Daphne Major Island. (Princeton Univ. Press, 2014)
Parchman, T. L., Buerkle, C. A., Soria‐Carrasco, V. & Benkman, C. W. Genome divergence and diversification within a geographic mosaic of coevolution. Mol. Ecol. 25, 5705–5718 (2016).
Parchman, T. L. et al. Resource stability and geographic isolation are associated with genome divergence in western Palearctic crossbills. J. Evol. Biol. 31, 1715–1731 (2018).
Lanner, R. M. Adaptation of whitebark pine for seed dispersal by Clark’s Nutcracker. Can. J. Forest Res. 12, 391–402 (1982).
Berggren, L. Occurrence of Swiss stone pine (Pinus cembra) in a boreal forest landscape in relation to regeneration characteristics and habitat preferences by nutcrackers (Nucifraga caryocatactes). Master thesis, Institutionen for skoglig zooekologi. (2005).
Neuschulz, E. L., Mueller, T., Bollmann, K., Gugerli, F. & Böhning-Gaese, K. Seed perishability determines the caching behavior of a food-hoarding bird. J. Anim. Ecol. 84, 71–78 (2015).
Haffer, J. in Handbuch der Vogel Mitteleuropas. (eds. Glutz von Blotzheim, U. N., Bauer, K.) vol. 13/III: Corvidae—Sturnidae. (AULA-Verlag, 1993).
Cramp, S. & Perrins, C. M. Handbook of the Birds of Europe, the Middle East and North Africa. Vol. VIII: Crows to Finches. (Oxford Univ. Press, 1994).
Tomback, D. F. & Linhart, Y. B. The evolution of bird-dispersed pines. Evol. Ecol. 4, 185–219 (1990).
Gernandt, D. S., López, G. G., García, S. O. & Liston, A. Phylogeny and classification of Pinus. Taxon 54, 29–42 (2005).
Linhart, Y. B. & Tomback, D. F. Seed dispersal by nutcrackers causes multi-trunk growth form in pines. Oecologia 67, 107–110 (1985).
Schütt, P., Schuck, H. J. & Stimm, B. Lexikon der Forstbotanik. (Ecomed Verlag, 1992).
Farjón, A. A. Handbook of the World’s Conifers. Vol II. (Brill, 2010).
Vander Wall, S. B. Food Hoarding in Animals, 1st edn. (University of Chicago Press, 1990).
Tomback, D. F. Dispersal of whitebark pine seeds by Clark’s nutcracker: a mutualism hypothesis. J. Anim. Ecol. 51, 451–467 (1982).
Reimers, N. F. Birds of the Cedar-pine forests of south-central Siberia and their role in the life of the cedar-pine. Tr. Biol. Inst. Sibirsk. Otd. Akad. Nauk. USSR 5, 121–166 (1959).
Mattes, H. Die Lebensgemeinschaft von Tannenhäher, Nucifraga caryocatactes (L.) und Arve, Pinus cembra L., und ihre forstliche Bedeutung in der oberen Gebirgswaldstufe (Flück-Wirth, 1982).
Hutchins, H. E., Hutchins, S. A. & Liu, B. W. The role of birds and mammals in Korean pine (Pinus koraiensis) regeneration dynamics. Oecologia 107, 120–130 (1996).
Kajimoto, T., Onodera, H., Ikeda, S., Daimaru, H. & Seki, T. Seedling establishment of subalpine stone pine (Pinus pumila) by nutcracker (Nucifraga) seed dispersal on Mt. Yumori, northern Japan. Arc. Alp. Res. 30, 408–417 (1998).
Hayashida, M. Seed dispersal of Japanese stone pine by the Eurasian Nutcracker. Ornithol. Sci. 2, 33–40 (2003).
Neuschulz, E. L., Merges, D., Bollmann, K., Gugerli, F. & Böhning‐Gaese, K. Biotic interactions and seed deposition rather than abiotic factors determine recruitment at elevational range limits of an alpine tree. J. Ecol. 106, 948–959 (2017).
Scherges, A. M. Ecomorphological adaptations to feeding and foraging behaviour in Eurasian Nucifraga nutcrackers. Master thesis, Rheinische Friedrich-Wilhelms-Universität Bonn. (2018).
Hao, Z. Z., Liu, Y. Y., Nazaire, M., Wei, X. X. & Wang, X. Q. Molecular phylogenetics and evolutionary history of sect. Quinquefoliae (Pinus): implications for Northern Hemisphere biogeography. Mol. Phylogenet. Evol. 87, 65–79 (2015).
Saladin, B. et al. Fossils matter: improved estimates of divergence times in Pinus reveal older diversification. BMC Evol. Biol. 17, 95 (2017).
Kleinschmidt, O. Corvus Nucifraga. (Berajah, Zoographica infinita, 1909–1911).
Mattes, H. in Engadin: Studien zu seiner Ökologie und Funktion im Arvenwald. 1–78 (Ferdinand Schöningh, 1978).
Lanner, R. M. Made for each other: a symbiosis of birds and pines. (Oxford Univ. Press, 1996).
Howard, R., Moore, A., Dickinson, E. & Remsen, J. The Howard and Moore complete checklist of the birds of the world. (Aves Press, 2003)
Clements, J. F. et al. The eBird/Clements Checklist of Birds of the World: v2019 (2019). https://www.birds.cornell.edu/clementschecklist/download/.
Gill, F., Donsker, D., Rasmussen, P., Eds. IOC World Bird List (v. 10.1). https://doi.org/10.14344/IOC.ML.10.1. (2020).
Madge, S., del Hoyo, J. Christie, D. A. Collar, N. & Kirwan G. M. In Birds of the World (eds. Billerman, S. M., Keeney, B. K., Rodewald, P. G., & Schulenberg, T. S.) (Cornell Lab of Ornithology, 2020).
del Hoyo, J. et al. HBW and BirdLife International Illustrated Checklist of the Birds of the World. Volume 2: Passerines. (Lynx Edicions and BirdLife International, 2016).
Haring, E., Gamauf, A. & Kryukov, A. Phylogeographic patterns in widespread corvid birds. Mol. Phylogenet. Evol. 45, 840–862 (2007).
Dohms, K. M. & Burg, T. M. Limited geographic genetic structure detected in a widespread Palearctic corvid Nucifraga caryocatactes. PeerJ 2, e371 (2014).
Ottenburghs, J. et al. A tree of geese: a phylogenomic perspective on the evolutionary history of True Geese. Mol. Phylogenet. Evol. 101, 303–313 (2016).
Suh, A. The phylogenomic forest of bird trees contains a hard polytomy at the root of Neoaves. Zool. Scr. 45, 50–62 (2016).
Toews, D. P. et al. Plumage genes and little else distinguish the genomes of hybridizing warblers. Curr. Biol. 26, 2313–2318 (2016).
Ottenburghs, J. et al. A history of hybrids? Genomic patterns of introgression in the True Geese. BMC Evol. Biol. 17, 201 (2017).
Campagna, L. et al. Repeated divergent selection on pigmentation genes in a rapid finch radiation. Sci. Adv. 3, e1602404 (2017).
Kryukov, A. P. et al. Deep phylogeographic breaks in magpie Pica pica across the Holarctic: concordance with bioacoustics and phenotypes. Zool. Sci. 34, 185–200 (2017).
Aoki, D. et al. Quaternary‑related genetic differentiation and parallel population dynamics of the Eurasian Jay (Garrulus glandarius) in the circum‑Japan Sea region. J. Ornithol. 159, 1087–1097 (2018).
Song, G. et al. Complete taxon sampling of the avian genus Pica (magpies) reveals ancient relictual populations and synchronous Late‐Pleistocene demographic expansion across the Northern Hemisphere. J. Avian Biol. 49, jav–01612 (2018).
Kryukov, A. P. Phylogeography and hybridization of corvid birds in the Palearctic Region. Vavilov J. Genet. Breed. 23, 232–238 (2019).
Fernando, S. W., Peterson, A. T. & Li, S. Reconstructing the geographic origin of the New World jays. Neotrop. Biodivers. 3, 80–92 (2017).
Töpfer, T. et al. A molecular phylogeny of bullfinches Pyrrhula Brisson, 1760 (Aves: Fringillidae). Mol. Phylogenet. Evol. 58, 271–282 (2011).
Tietze, D. T., Martens, J., Sun, Y. ‐H., Severinghaus, L. & Päckert, M. Song evolution in the coal tit Parus ater. J. Avian Biol. 42, 214–230 (2011).
Päckert, M. et al. Horizontal and elevational phylogeographic patterns of Himalayan and Southeast Asian forest passerines (Aves: Passeriformes). J. Biogeogr. 39, 556–573 (2012).
Päckert, M. et al. A revised phylogeny of nuthatches (Aves, Passeriformes, Sitta) reveals insight in intra- and interspecific diversification patterns in the Palearctic. Vertebr. Zool. 70, 241–262 (2020).
Päckert, M. et al. Out of Tibet and out of Himalayas: centers of origin and diversification of passerine birds in the Qinghai-Tibet Plateau region. Ecol. Evol. 10, 9283–9300 (2020).
Alström, P. et al. Complete species-level phylogeny of the leaf warbler (Aves: Phylloscopidae) radiation. Mol. Phylogenet. Evol. 126, 141–152 (2018).
Fok, K. W., Wade, C. M. & Parkin, D. T. Inferring the phylogeny of disjunct populations of the azure-winged magpie Cyanopica cyanus from mitochondrial control region sequences. Proc. R. Soc. B 269, 1671–1679 (2002).
Päckert, M. et al. The great tit (Parus major)—a misclassified ring species. Biol. J. Linn. Soc. Lond. 86, 153–174 (2005).
Zhang, R. et al. Comparative phylogeography of two widespread magpies: importance of habitat preference and breeding behavior on genetic structure in China. Mol. Phylogenet. Evol. 65, 562–572 (2012).
Tritsch, C. et al. Improved sampling at the subspecies level solves a taxonomic dilemma—a case study of two enigmatic Chinese tit species (Aves, Passeriformes, Paridae, Poecile). Mol. Phylogenet. Evol. 107, 538–550 (2017).
Salinas, P. et al. Genetic diversity, differentiation and historical origin of the isolated population of rooks Corvus frugilegus in Iberia. J. Avian Biol. 52, e02689 (2021).
Hollyer, J. N. The invasion of Nutcrackers in autumn 1968. Brit. Birds 63, 353–373 (1968).
Newton, I. The Migration Ecology of Birds. (Academic Press, 2008).
Nazarenko, A. A. On faunistic cycles (extinction–expansion–extinction) with special reference to the East Palearctic dendrophilous avifauna. Zh. Obshei. Biol. 43, 823–835 (1982).
Kahlke, H. D. Die Eiszeit. (Urania Verlag, 1994).
Nazarenko, A. A. Ornithofaunistical exchange through the South and the North Asia on the eastern periphery of the continent—the last glacial-interglacial cycle. Zh. Obshei Biol. 51, 89–106 (1990).
Manish, K. & Pandit, M. K. Geophysical upheavals and evolutionary diversification of plant species in the Himalaya. PeerJ 6, e5919 (2018).
Päckert, M. et al. Evolutionary history of passerine birds (Aves: Passeriformes) from the Qinghai–Tibetan plateau: from a pre-Quarternary perspective to an integrative biodiversity assessment. J. Ornithol. 156, 355–365 (2015).
Frankham, R. Do island populations have less genetic variation than mainland populations? Heredity 78, 311–327 (1997).
Benkman, C. W., Parchman, T. L., Favis, A. & Siepielski, A. M. Reciprocal selection causes a coevolutionary arms race between crossbills and lodgepole pine. Am. Nat. 162, 182–194 (2003).
Thompson, J. N. The Geographic Mosaic of Coevolution (University of Chicago Press, 2005).
Löhrl, H. Der Tannenhäher (Nucifraga carycatactes) beim Sammeln und Knacken von Nüßchen. Anz. Orn. Ges. Bayern 9, 187–196 (1970).
Lamichhaney, S. et al. Evolution of Darwin’s finches and their beaks revealed by genome sequencing. Nature 518, 371–375 (2015).
Han, F. et al. Gene flow, ancient polymorphism, and ecological adaptation shape the genomic landscape of divergence among Darwin’s finches. Genome Res. 27, 1004–1015 (2017).
del Hoyo, J. et al. HBW and BirdLife International Illustrated Checklist of the Birds of the World. Volume 2: Passerines. (Lynx Edicions and BirdLife International, 2020).
Ruguzova, G., Pokhylchenko, O., Ivanova, I. & Yaremchuk, Y. Some features of Pinus pumila (Pall.) Reg. and Pinus armandii Franch. Seed formation in the conditions of introduction. Agric. Forest. 62, 203–212 (2016).
Andrews, S. FastQC: a quality control tool for high throughput sequence data. http://www.bioinformatics.babraham.ac.uk/projects/fastqc. (2010).
Weisenfeld, N. I., Kumar, V., Shah, P., Church, D. M. & Jaffe, D. B. Direct determination of diploid genome sequences. Genome Res. 27, 757–767 (2017).
Gurevich, A. et al. QUAST: quality assessment tool for genome assemblies. Bioinformatics 29, 1072–1075 (2013).
Seppey, M., Manni, M. & Zdobnov, E. M. BUSCO: assessing genome assembly and annotation completeness. Methods Mol. Biol. 1962, 227–245 (2019).
Chu, J. Jupiter Plot: A Circos-based tool to visualize genome assembly consistency (Version 1.0). https://zenodoorg/record/1241235XA92q2hKiUk. (2017).
Smit, A. & Hubley, R. RepeatModeler Open-1.0. http://www.repeatmasker.org. (2015).
Vijay, N. et al. Evolution of heterogeneous genome differentiation across multiple contact zones in a crow species complex. Nat. Commun. 7, 13195 (2016).
Smit, A. et al. RepeatMasker Open-4.0. http://www.repeatmasker.org. (2013).
Mundy, N., Unitt, P. & Woodruff, D. Skin from feet of museum specimens as a non-destructive source of DNA for avian genotyping. Auk 114, 126–129 (1997).
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).
Zhang, J., Kobert, K., Flouri, T. & Stamatakis, A. PEAR: a fast and accurate Illumina Paired-End reAd mergeR. Bioinformatics 30, 614–620 (2014).
Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics 25, 1754–1760 (2009).
Garrison, E. & Marth G. Haplotype-based variant detection from short-read sequencing. arXiv preprint arXiv:1207.3907 [q-bio.GN] (2012).
Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A. & Jermiin, L. S. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat. Methods 14, 587–589 (2017).
Nguyen, L. T., Schmidt, H., Haeseler, A. & Minh, B. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268–274 (2015).
Kumar, V. et al. The evolutionary history of bears is characterized by gene flow across species. Sci. Rep. 7, 46487 (2017).
Shimodaira, H. An approximately unbiased test of phylogenetic tree selection. Syst. Biol. 51, 492–508 (2002).
Zhang, C., Rabiee, M., Sayyari, E. & Mirarab, S. ASTRAL-III: polynomial time species tree reconstruction from partially resolved gene trees. BMC Bioinformatics 19, 153 (2018).
Huson, D. H. & Bryant, D. Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 23, 254–267 (2006).
Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007).
Alexander, D. H. & Lange, K. Enhancements to the ADMIXTURE algorithm for individual ancestry estimation. BMC Bioinformatics 12, 246 (2011).
Bouckaert, R. et al. BEAST 2.5: An advanced software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 15, e1006650 (2019).
Bouckaert, R. R. & Drummond, A. J. bModelTest: Bayesian phylogenetic site model averaging and model comparison. BMC Evol. Biol. 17, 42 (2017).
Nakamura, T., Yamada, K. D., Tomii, K. & Katoh, K. Parallelization of MAFFT for large-scale multiple sequence alignments. Bioinformatics 34, 2490–2492 (2018).
Oliveros, C. H. et al. Earth history and the passerine superradiation. Proc. Natl Acad. Sci. USA 116, 7916–7925 (2019).
Moyle, R. G. et al. Tectonic collision and uplift of Wallacea triggered the global songbird radiation. Nat. Commun. 7, 12709 (2016).
Benton, M. J., Donoghue, P. C. J. & Asher, R. In The Timetree of Life (eds. Hedges, S. B. and Kumar, S.) 35–86 (Oxford Univ. Press, 2009).
Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior Summarization in Bayesian Phylogenetics Using Tracer 1.7. Syst. Biol. 67, 901–904 (2018).
Li, H. & Durbin, R. Inference of human population history from individual whole-genome sequences. Nature 475, 493–496 (2011).
Nadachowska-Brzyska, K., Burri, R., Smeds, L. & Ellegren, H. PSMC analysis of effective population sizes in molecular ecology and its application to black-and-white Ficedula flycatchers. Mol. Ecol. 25, 1058–1072 (2016).
BirdLife International Species factsheet: Nucifraga caryocatactes. http://www.birdlife.org (2020).
Paradis, E. & Schliep, K. ape 5.0: an environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528 (2019).
R Core Team. R: A language and environment for statistical computing (R Foundation for Statistical Computing, 2020).
Beck, M. W. ggord: Ordination Plots with ggplot2. R package version 1.0.0. https://zenodo.org/badge/latestdoi/35334615 (2017).
Nilsson Janke, Maria et al. Supporting data: Speciation and population divergence in a mutualistic seed dispersing bird, Dryad, Dataset, https://doi.org/10.5061/dryad.stqjq2c57 (2022).
Acknowledgements
We thank Dr. Stephan Weigl, Biologiezentrum Linz, Dr. Jocelyn Hudon, Royal Alberta Museum, Canada, Dr. Jan Bolding Kristensen, Natural History Museum of Denmark, Dr. Shigeki Asai, Yamashina Institute for Ornithology, Japan, Dr. Lars Erik Johannessen, Natural History Museum, Norway, Dr. Ulf Johansson, Natural History Museum, Sweden and Dr. Sylke Frahnert of the Museum für Naturkunde in Berlin for providing samples from their collections. Additionally, we express our gratitude to Dr. Mozes Blom, Dr. Fritjof Lammers and Dr. Tilman Schell for their support and feedback regarding bioinformatics. We thank Alice Scherges and Prof. Susanne Fritz for the vegetation maps. We thank Jon Baldur Hlidberg (www.fauna.is) for the bird images. We thank the Erika und Walter Datz-Stiftung for providing additional financial support for the project. Travel funding for China was provided by DFG PA1818/3-1 to M.P. and by Feldbausch-Stiftung and Wagner-Stiftung to J.M. for field studies in Asia (both Johannes Gutenberg-University of Mainz, Germany). E.L.N. was supported by the German Research Foundation (NE1863/2-1, NE1863/2-2) and M.I. was supported by the Swedish Research Council (2019-03900). Part of the work done by Y.A.R. was carried out under the state program “Taxonomic and Biochronological Analysis of the Animal World as a Basis for Studying and Preserving the Structure of Biological Diversity” (number 121032300105-0). Additionally, we acknowledge support from the National Genomics Infrastructure in Stockholm funded by Science for Life Laboratory, the Knut and Alice Wallenberg Foundation and the Swedish Research Council and SNIC/Uppsala Multidisciplinary Center for Advanced Computational Science for assistance with massively parallel sequencing and access to the UPPMAX computational infrastructure. Lastly, the present study is a result of the Centre for Translational Biodiversity Genomics (LOEWE-TBG) and was supported through the programme ‘LOEWE – Landes-Offensive zur Entwicklung Wissenschaftlich-ökonomischer Exzellenz’ of Hesse’s Ministry of Higher Education, Research and the Arts.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Contributions
M.N. and J.d.R. designed the study. M.I., M.N. and J.d.R. performed laboratory procedures and sequencing. J.d.R. led the bioinformatic processing, data analyses and interpretation with contributions from M.N., E.L.N., M.S. and M.P. E.L.N., M.S. and T.T. provided the morphometric dataset. M.P. and J.M. assisted in the biogeographic interpretation of the data. M.P., M.I., J.M., A.P.K., Y.A.R., Y.S., E.L.N. and A.J. provided samples and funding. M.N. and J.d.R. wrote the manuscript with input from all co-authors. All authors gave final approval for publication.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Peer review
Peer review information
Communications Biology thanks Jente Ottenburghs and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Luciano Matzkin and Luke R. Grinham.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
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 license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license 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 license, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
de Raad, J., Päckert, M., Irestedt, M. et al. Speciation and population divergence in a mutualistic seed dispersing bird. Commun Biol 5, 429 (2022). https://doi.org/10.1038/s42003-022-03364-2
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s42003-022-03364-2
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.
