Genome sequencing and assembly
Genome assemblies ranged in size from 799.9 Mbp in Melanocharis versteri to 1053.5 Mbp in Sericornis nouhuysi. The number of scaffolds ranged from 14,086 scaffolds in Melipotes ater to 87,957 scaffolds in Ficedula hyperythra and N50 ranged between ca. 40 Kbp to and 25 Mbp. Benchmarking Universal Single-Copy Orthologs (BUSCO) analyses of genome completeness ranged from a high proportion of complete BUSCOs in Melipotes ater, 86.8% to only 66.7% complete BUSCOs in Rhipidura albolimbata. For most species, the proportions of complete BUSCOs were 75–80%. Overall, the proportion of missing BUSCOs was low, ranging from 6.6% in Melipotes ater to 15.2% in Rhipidura albolimbata (see Supplementary Table 1 for all genome assembly statistics and Supplementary Fig. 1 for the number of SNP variants per species).
Kinship analyses of individuals within populations
Sampling of closely related individuals can dramatically bias estimates of population structure and demographics. Two Pachycephala schlegelii individuals (A117 and A118) showed a pairwise kinship coefficient of 0.144, indicative of being half-siblings. The two individuals were collected at the same locality on the same date. Similarly, two Ifrita kowaldi individuals (D116 and D117) showed a pairwise kinship coefficient of 0.135, also suggestive of being half-siblings. In this case, the individuals were collected on the same sampling locality on two consecutive days. To not bias downstream demographic analyses, one of the P. schlegelii (A118) and one of the I. kowaldi (D117) individuals were excluded from all subsequent analyses. For all other species, no closely related individuals were identified.
Genetic differentiation
Estimated levels of differentiation between populations were initially based on three approaches; (i) calculation of FST (the fixation index), which quantifies the degree of genetic differentiation between populations based on the variation in allele frequencies, ranging between 0 (complete mixing of individuals) and 1 (complete separation of populations) (Fig. 1), (ii) Standardized pairwise FST used to conduct a Principal Component Analysis (PCA) in order to visualize population structure (Supplementary Fig. 1) and (iii) Admixture analysis as implemented in STRUCTURE (a clustering algorithm that infers the most likely number of groups [K]), in which individuals are grouped into clusters according to the proportion of their ancestry components (Supplementary Fig. 1). As a preliminary analysis, we calculated FST and constructed PCA plots for the four congeneric (incl. Sericornis/Aethomyias [until recently placed in the genus Sericornis]) species pairs in our study (Supplementary Fig. 2), which were aligned using the same reference genome. This was done to ascertain that no samples had been misidentified and to gauge levels of differentiation between distinct species. All species were genetically well separated and FST values ranged from 0.08 for the two Ptiloprora species to 0.20 for the two Ficedula species.
For five out of six species from Buru/Seram, genetic differentiation (FST) was high between islands (Fig. 1), and comparable to differentiation between named congeneric species in this study (e.g. Ptiloprora and Melipotes); Ceyx lepidus (FST = 0.16), Thapsinillas affinis (FST = 0.15), Ficedula buruensis (FST = 0.13) and Pachycephala macrorhyncha (FST = 0.09). In contrast, differentiation in Ficedula hyperythra was consistent with population-level differentiation (FST = 0.04). In all cases, individuals from Buru and Seram were clearly differentiated in the PCA and STRUCTURE plots (Supplementary Fig. 1A). For Ceyx lepidus, Ficedula buruensis and Pachycephala macrorhyncha, samples were collected at multiple elevations and we therefore calculated genetic differentiation between elevations (Buru: 1097 m versus 1435 m and Seram: 1000 m versus 1300 m) to determine any potential parapatric differentiation along the gradients. In all possible comparisons, FST values did not differ significantly from 0. Moreover, PCA plots showed that samples did not cluster according to elevation (Supplementary Fig. 3A).
Three of the thirteen New Guinean population pairs occurring in Mount Wilhelm and Huon showed relatively high genetic divergences: Melipotes fumigatus/ater (FST = 0.08), Paramythia montium (FST = 0.09) and Ifrita kowaldi (FST = 0.07) (Fig. 1) with populations clearly separated (Supplementary Fig. 1). By contrast, the two lowland species Toxorhamphus novaeguineae and Melilestes megarhynchus showed little genetic differentiation, FST = 0.00. For the remaining species, genetic differentiation between Mount Wilhelm and Huon ranged between FST = 0.01–0.05. Despite this moderate level of genetic differentiation, the populations of Mount Wilhelm and Huon could be clearly distinguished in the PCA plots. In all cases STRUCTURE suggested a scenario with K = 2 with some mixing of individuals, except for Rhipidura albolimbata, in which K = 1 was suggested.
For five bird species we included an additional population from Mount Scratchley, which is also situated in the Central Range but ~400 km to the southeast of Mount Wilhelm. Genetic differentiation of this population from the other two populations was comparable with that between Mount Wilhelm and Huon. The highest genetic differentiation was found in Paramythia montium (FST = 0.10 both between Mount Wilhelm and Mount Scratchley and between Huon and Mount Scratchley). In the case of Peneothello sigillata, the Mount Scratchley population appeared genetically well-differentiated from both the populations of Mount Wilhelm (FST = 0.06) and Huon (FST = 0.07). In both cases, STRUCTURE suggested a scenario of K = 3, with individual assignments matching the three geographically circumscribed populations. For Pachycephala schlegelii, genetic differentiation was relatively high between Huon and Mount Scratchley (FST = 0.05), but low between Mount Wilhelm and Mount Scratchley (FST = 0.01). Accordingly, STRUCTURE suggested a scenario with K = 2 groups. For the remaining two species Sericornis nouhuysi showed some differentiation (FST = 0.03) between Mount Wilhelm and Huon and Aethomyias papuensis showed minor differentiation (FST = 0.02 between Mount Scratchley and Huon (Supplementary Table 2), but for both species, STRUCTURE suggested a scenario of K = 2 with considerable mixing of individuals between populations.
Samples from Mount Wilhelm were collected at elevations ranging from 1700 to 3700 m, again allowing us to test for differences within populations on a single slope, a finding that would be consistent with incipient parapatric speciation. No species showed significant differences in FST when comparing individuals from different elevations, and concordantly there was little clustering of individuals by elevation in the PCA plots. Even when individuals were collected as far as 2000 elevational meters apart (as in the case of Origma robusta), genetic differentiation was low (FST = 0.01). In Huon, all samples were collected at the same elevation, except for Ifrita kowaldi, for which genetic differentiation of FST = 0.03 was found between individuals collected at 2300 m and 2950 m (Supplementary Fig. 3B, Supplementary Table 2). These analyses however, suffer from very small sample sizes that hinder a thorough analysis of parapatric speciation events. Furthermore, we note that divergence with gene flow may not manifest as a genome-wide phenomenon (at least, not until the taxa are so differentiated that gene flow has ceased). Instead, it may proceed via selection acting to create small ‘islands of differentiation’ within the genome against a background of negligible differentiation22,23. Such analyses require large sample sizes and are therefore not possible herein.
Correlations between genetic divergence and elevation
If lineages colonize mountains from the lowlands, followed by range contraction and differentiation in the highlands, we would expect a signature of larger genetic differentiation (FST) between populations inhabiting higher elevations. We found no relationship between genetic differentiation (FST) and the altitudinal floor (the lowest elevation at which a species/population occurs) for the five Moluccan species, but for all New Guinean taxa with the exception of Melipotes fumigatus/ater we found a significant positive correlation (r = 0.83, p < 0.001). This relationship remained significant even when excluding the two lowland taxa (Toxorhamphus novaeguineae and Melilestes megarhynchus, r = 0.70, p = 0.022).
Demographic inferences using Pairwise Sequentially Markovian Coalescent (PSMC) analyses
The demographic history of the eighteen bird species in our study was first inferred by analyzing the whole-genome sequence of one single individual per species using Pairwise Sequentially Markovian Coalescent (PSMC)24 (Supplementary Fig. 4). PSMC relies on the distribution of heterozygous sites across the genome and infers the distribution of the time since the most recent common ancestor (TMRCA) between each pair of alleles at all loci across the whole genome of one single individual. This provides information about how effective population sizes change over time.
The inferred demographic histories of the eighteen bird species encompassed the time period from ca. 10 Mya to 20 Kya (note that most of the species included herein are younger than 10 My and as such past demographic dynamics reflect an ancestral form, likely with a very different distribution, Supplementary Fig. 4). Initial population sizes for most species were ca. 300,000–400,000 individuals, with the lowest initial population size being Ifrita kowaldi with ca. 200,000 individuals, and the highest being Toxorhamphus novaeguineae and Peneothello sigillata with ca. 1,000,000 individuals.
Inferred demographic fluctuations over time show that most species had an initial period of relative demographic stability during the Miocene (5.3–23 Mya) and Pliocene (2.6–5.3 Mya). During this period, most species showed a rather flat curve. The most dramatic demographic changes for the majority of species occurred in the Pleistocene (beginning around 2.6 Mya) but the timing of population increases and decreases vary between species. All PSMC results were consistent when bootstrapping the data but clearly become more dubious toward the present (Supplementary Fig. 4).
PSMC plots generated for the re-sequenced individuals (Figs. 2–4) allowed us to compare the demographic histories of the different populations—those on Mount Wilhelm, Huon, and Mount Scratchley in New Guinea; and on Buru and Seram in the Moluccas. These plots closely mimicked the demographic changes inferred from the de novo assembled reference sequences (Supplementary Fig. 4). Individuals from the same population showed largely similar demographic trajectories and only differed for the most recent time periods. For three New Guinean species (Melanocharis versteri, Toxorhamphus novaeguineae and Melilestes megarhynchus), populations from Mount Wilhelm and Huon shared the same demographic history over time. However, for all other New Guinean species, populations from Mount Wilhelm and Huon showed contrasting demographic patterns from the mid-late Pleistocene towards the present: populations from Huon experienced an important demographic decline, dropping to a historical low by the Late Pleistocene, while populations from Wilhelm remained stable and did not experience the decline observed in Huon. This contrasting demographic pattern can be mostly clearly observed in Origma robusta, Ifrita kowaldi, Rhipidura albolimbata, Ptiloprora guisei and Melipotes fumigatus/ater. Similarly, populations of the species occurring on Mount Wilhelm, Huon and Mount Scratchley showed comparable demographic histories up until the mid-late Pleistocene, at which point each population began to behave slightly differently, with Huon always showing the lowest population sizes.
The brown curve is the PSMC estimate for the de novo sequence data (for PSMC estimates for 100 bootstrapped sequences see Supplementary Fig. 4). The PSMC estimates for the resequencing data are presented for Ficedula buruensis, Ficedula hyperythra, Pachycephala macrorhyncha, Thapsinillas affinis and Ceyx lepidus) with populations from Buru (red) and Seram (blue). The most recent 20 Ky have been removed from the PSMC plots and the current Ne estimate indicated at T = 10 Kya is based on the demographic analyses in fastsimcoal2. The shaded area indicates the Pleistocene. Illustrations of the focal bird species from del Hoyo et al.68.
The brown curve is the PSMC estimate for the de novo sequence data (for PSMC estimates for 100 bootstrapped sequences see Supplementary Fig. 4). The PSMC estimates for the resequencing data are presented for Melipotes fumigatus/ater, Ifrita kowaldi, Origma robusta, Ptiloprora guisei, Rhipidura albolimbata and Melanocharis versteri with populations from Mount Wilhelm (purple) and Huon (blue). The most recent 20 Ky have been removed from the PSMC plots and the current Ne estimate indicated at T = 10 Kya is based on the demographic analyses in fastsimcoal2. The shaded area indicates the Pleistocene. Illustrations of the focal bird species from del Hoyo et al.68.
The brown curve is the PSMC estimate for the de novo sequence data (for PSMC estimates for 100 bootstrapped sequences see Supplementary Fig. 4). The PSMC estimates for the resequencing data are presented for Melilestes megarhynchus, Toxorhamphus novaeguineae including a third population from the lowlands north of the Ramu/Markham River (light green) and New Guinean species (Aethomyias papuensis, Peneothello sigillata, Pachycephala schlegelii, Paramythia montium and Sericornis nouhuysi) including a third montane population from Mount Scratchley (dark green) in the south-eastern central range. Populations from Mount Wilhelm (purple) and Huon (blue). The most recent 20 Ky have been removed from the PSMC plots and the current Ne estimate indicated at T = 10 Kya is based on the demographic analyses in fastsimcoal2. The shaded area indicates the Pleistocene. Illustrations of the focal bird species from del Hoyo et al.68.
Of the species occurring in the Moluccas, Ficedula hyperythra and Thapsinillas affinis showed a common demographic history on both Buru and Seram. The two Ceyx lepidus populations also appear to have a similar demographic history up until 1 Mya, when the effective population size spiked on Buru, but not on Seram (Fig. 2).
For the three species with no de novo genome (Aethomyias papuensis, Pachycephala macrorhyncha and Ficedula buruensis), we conducted PSMC analyses with re-sequenced data mapped to the following non-conspecific reference genomes: Sericornis nouhuysi (for the closely related Aethomyias papuensis), Pachycephala schlegelii (for the closely related Pachycephala macrorhyncha) and Ficedula hyperythra (for the closely related Ficedula buruensis). While the PSMC plots for these three species should be interpreted with some caution, we observed the same large population expansions during the Pleistocene, but the expansion occurred in the late- rather than the mid-Pleistocene.
Overall, the PSMC plots obtained from resequencing data (Figs. 2–4) were largely concordant with the PSMC plots obtained from de novo sequencing (Supplementary Fig. 4). However, plots from resequencing data differed in showing variable trajectories from 20 to 30,000 years ago to the present, suggesting that recent demographic changes and current population sizes inferred with PSMC are uncertain and should be interpreted with caution. In contrast, the demographic histories of the species between further back in time appeared robust, with good agreement between the PSMC plots derived from the de novo and resequencing data.
Demographic inference using explicit model testing
The demographic history of the populations of pairs of the eighteen bird species in our study was further investigated using the model-testing approach implemented in fastsimcoal225. We present a summary of the results in Table 1 and the results from all models tested in Supplementary Data 1, including current and ancestral population sizes, migration rates and divergence times and summarize the major findings below.
Models including migration always performed better than models without migration. Migration rates between Buru and Seram were low for the five species compared, ranging between proportions of 10−6 and 10−7 per generation. Migration rates between Mount Wilhelm and Huon were more variable, and while low migration rates (ca. 10−6–10−7) were suggested for those species with high genetic differentiation (i.e. Melipotes fumigatus/ater and Ifrita kowaldi), higher migration rates (ca. 10−5) were suggested for the lowland species Toxorhamphus novaeguineae and Melilestes megarhynchus and for those species with low genetic differentiation (e.g. Melanocharis versteri, Rhipidura albolimbata and Ptiloprora guisei). For a few species, migration rates were close to zero (10−11–10−18), but we found no clear directionality in migration rates between Mount Wilhelm and Huon.
Ancestral population sizes were similar across species, with an average of ca. 325,000 individuals, ranging from ca. 100,000 in Thapsinillas affinis to ca. 600,000 individuals in Ficedula buruensis (see Supplementary Data 1 for confidence intervals). Current population sizes for the New Guinean species were generally higher in Mount Wilhelm than Huon, with an average population size of ca. 400,000 individuals per species for the former, and an average of ca. 200,000 individuals for the latter. Current population sizes ranged from ca. 75,000 individuals in Aethomyias papuensis to ca. 900,000 in Rhipidura albolimbata in Mount Wilhelm, and from ca. 35,000 individuals in Paramythia montium to ca. 700,000 in M. megarhynchus in Huon. Current population sizes in Mount Scratchley were also generally lower than in Mount Wilhelm, but similar to the values found in Huon, ranging from ca. 25,000 individuals in Peneothello sigillata to ca. 400,000 individuals in Pachycephala schlegelii. For the Moluccan species, current population sizes were higher in Seram, with an average of ca. 225,000 individuals per species, relative to Buru, with an average of ca. 175,000 individuals. Current population sizes ranged from ca. 140,000 individuals in Ceyx lepidus to ca. 350,000 in Ficedula buruensis in Seram, and from 50,000 individuals in Ceyx lepidus to 500,000 in Ficedula buruensis in Buru. Divergence times between populations in the Moluccan islands ranged between 0.1 and 1.4 My (mean = 0.583 My) and between 0.1 and 1.8 My (mean = 1.2 My) for the New Guinean populations.
Source: Ecology - nature.com
