Introduction

Tropical islands of the Indo-Australian Archipelago (IAA), the Indo-Pacific, and the Caribbean and Atlantic have had complex geological pasts that have affected their species distribution patterns. While several Pacific islands are true oceanic islands arising de novo from the seafloor, many islands of the Caribbean and the IAA were once connected with each other or to the continental landmasses during periods of global sea-level fall in the Quaternary (Lohman et al., 2011; Voris, 2000). These land bridges and the additional habitat they offered during periods of sea-level fall facilitated on one hand colonisations among islands and between islands and the mainland (Andersen et al., 2015; Cros et al., 2020; Irestedt et al., 2013; Moyle et al., 2009; Ng et al., 2017; Pujolar et al., 2022) and on the other, range expansions within islands in response to the fluctuations in suitable habitat space. The habitat type and quality of these intermittently exposed land bridges varied, ranging from open savannah biomes to lowland forests (Cannon et al., 2009; Lohman et al., 2011). In the case of stratified elevational gradients, with tropical upland montane forests expanding onto lower slopes during periods of global cooling, upland montane species could disperse across otherwise persistent dual barriers of land and sea (Cannon et al., 2009).

Tropical islands’ geologic pasts make them an important natural system to study species responses to changing habitat at both ecological and evolutionary timescales. Birds are an ideal clade for this because they are taxonomically well-characterised, are ecologically well-studied, and charismatic. Avian lineages have arisen on tropical islands during the Pleistocene (Andersen et al., 2015; Garg et al., 2018; Irestedt et al., 2013; Moyle et al., 2009), concomitant with complex, clade-wise colonisation and recolonisation trajectories (Filardi and Moyle, 2005;Jønsson et al., 2008, 2011b, 2014) including colonisations of the mainland (Jønsson et al., 2011a) and repeated, independent colonisation events (Cibois et al., 2011). In birds, these responses to changing habitat are governed by species traits in several continental taxa and non single-island endemics (Brüniche-Olsen et al., 2021, 2019), but responses in single-island endemics are poorly understood.

With such complex biogeographic pasts, tropical, single-island endemics represent either refugial populations, or in-situ diversifications as explained by taxon cycles (Ricklefs, 1970; Ricklefs and Cox, 1978, 1972). Being confined to single islands, these species are exceptionally vulnerable to ongoing anthropogenic climate change which is unprecedented and unlike past Pleistocene and early Holocene climatic changes (Crowley, 1990). This is exacerbated because we are currently experiencing a period of relatively high sea levels with reduced habitat availability for species resulting in population bottlenecks across taxa (Hewitt, 2000; Willis et al., 2004) including birds (Nadachowska-Brzyska et al., 2015; Smith et al., 2021). Information on tropical, single-island endemics’ demographic responses to past climate change can inform conservation efforts, owing to the genomic signatures that predispose a species to extinction (Mays et al., 2018; Spielman et al., 2004).

Pairwise Sequential Markovian Coalescent (PSMC) is a powerful method to reconstruct the effective population size (Ne) for a species using a single diploid genome (Li and Durbin, 2011) and has been used to infer demographic history across taxa (Chattopadhyay et al., 2019; Kim et al., 2016; Kozma et al., 2016; Murray et al., 2017; Nadachowska-Brzyska et al., 2015). Using this along with paleo-ecological niche modelling allows us to directly correlate the demographic history of a species with its distributional range at different time points (Chattopadhyay et al., 2019). We perform these analyses on a global panel of tropical single-island endemics to understand the effects of past climatic changes on species with assumed low genetic diversity owing to their status as refugial species or newly arisen lineages.

Results

Bird species panel

Our final species panel comprises 30 tropical single-island endemic bird species. Out of these, PSMC analyses were possible for 22 species and Ecological Niche Modelling (ENM) was possible for 28 species. Both PSMC and ENM analyses could be done successfully for 20 species (Table 1). These included seven Papuan species, five Philippine species, three Caribbean species, and three species from the Bismarck Archipelago. Of these 20 species, 13 species are passerines, with two from the white-eye family (Zosteropidae). Of the non-passerines, two species are parrots (family Psittacidae). All other families are represented by single species.

Table 1:

Paleo-habitat reconstruction

Among the 19 climatic variables used, we observed precipitation of the warmest quarter (BIO18) to contribute the most for each island endemic bird species except for the Caribbean species Amazona guildingii. 24 out of the 28 bird species experienced an increase in suitable habitat from the Last Interglacial (LIG) to the Last Glacial Maxima (LGM) (Supplementary table S4). 11 out of 13 Papuan species showed an increase in habitat availability from the LIG to the LGM, with Cicinnurus regius and Pseudorectus ferrugineus as exceptions. All Philippine species showed an increase in habitat from the LIG to the LGM as well (Figure 1, Supplementary information S2). Two out of the three Caribbean species experienced a decline in suitable habitat during this period. Of the remaining species, all showed an increase in available habitat. Upland Papuan species remained confined to upland montane regions at all the time points we reconstructed available habitat, with Rhagologus leucostigma as a notable exception (Supplementary information S2). Papuan lowland species remained largely confined to regions which were not newly-exposed landbridges at all times as well. From the LGM to the present day, habitat decreased for 19 out of the 28 bird species, and two species experienced no change (Supplementary table S3). For the present study, we primarily concentrate on habitat fluctuation between the LIG and the LGM as this is the period for which we have comparative evidence of fluctuations of both habitat as well as genetic diversity.

Figure 1:

Strong Quaternary fluctuations in genetic diversity

PSMC analyses could be successfully done for 22 species (Table 1). Results indicate large genetic diversity variations in almost all bird species, corresponding to the Quaternary climate shifts (Figure 2). The demographic history of these species reconstructed using Pairwise Sequentially Markovian Coalescent (PSMC) analyses extends back over a million years (Supplementary information S1). Reconstructed Ne values were generally concordant across the three PSMC settings used, except for five species where Ne values had large differences across PSMC settings for Ne values at the LIG and LGM. However, in these species as well, trends of Ne increase or decrease from the LIG to LGM were robust across all three PSMC models considered. Centropus unirufus, Dicaeum eximium, Irena cyanogastra, Sterrhoptilus dennistouni, Zosterops hypoxanthus were the six species which had large differences (Supplementary information S1). We found no significant difference between the different sets of PSMC settings used (Kruskal-Wallis test, χ²= 0.217, DF = 2, p = 0.897).

Figure 2:

The 14 passerine species we analysed showed overall lower Ne values in the LGM as compared to the LIG (Figure 1). We obtained similar results for non-passerines as well, with Rhynochetos jubatus as a notable exception using the –p “4 + 30 * 2 + 4 + 6 + 10” setting, displaying a large peak followed by a crash in Ne (Supplementary information S1). This is likely a PSMC artefact, because plots using other parameter settings did not show peak (Supplementary information S1). Along with overall lower values, we also see a smaller range of Ne values in the LGM as compared to the LIG (Figure 1).

Correlation between historical fluctuations in Ne and distribution

Habitat change was poorly associated with change in Ne for the 20 species for which both PSMC and ENM analyses were possible (Cramer’s V = 0.15). However, passerine species only showed a strong association (Cramer’s V = 0.96), while non-passerines showed a weak negative association (Cramer’s V = –0.15).

Bayesian multivariate regression models revealed an overall fluctuation in Ne was associated with species biology and habitat change during the Last Glacial Period (LGP) even after controlling for the geographical island group (Figure 3). Most confidence intervals and beta coefficients did not overlap with zero (Figure 3) suggesting significant relationships. We found a positive association between change in habitat area (LIG to LGM) and Ne (β = 9.45, 95% CI: [1.95, 21.1]). This means that species with increasing habitat area in the LGP also showed an increase in Ne. Change in Ne and body mass showed a significant negative association (β = −8.63, 95% CI: [16.27, –2.61]). That is, large-bodied species showed decreases in Ne during the LGP. Diet specialists like frugivores (β = –11.89, 95% CI: [–22.65, –2.99]) and invertivores (β = –11.7, 95% CI: [–22.66, –2.94]) also showed significant negative associations with population change. However, omnivores (β = –1.39, 95% CI: [–13.38, 11.5]) did not show any association with population change. Whether or not a species was a passerine was an important predictor of Ne only in combination with the change in habitat from LIG to LGM (β = 7.97, 95% CI: [1.82, 16.91], suggesting that passerines respond positively to habitat change, as suggested by Cramer’s V as well. The interaction between habitat change and body mass (β = 10.05, 95% CI: [–0.3, 24.41) suggests that there is no impact of habitat change in larger species. Finally, the random intercept for Country (sd (Intercept)) showed a marginal positive influence (β = 0.85, 95% CI: [0.04, 2.24]) (Figure 3, Supplementary table S6)

Figure 3:

Discussion

Island birds are vulnerable to the effects of climate change

Our results on the fluctuations in paleohabitat and genetic diversity (Ne) during the LGP highlight the vulnerability of island birds to climatic fluctuations. Climate fluctuations change habitat availability and, in turn, genetic diversity, making these species particularly sensitive to a changing environment. In response to landmass expansion, several species in our panel across islands showed an increase in their genetic diversity during the LGP (Figure 1, Supplementary information S1), although the overall relationship is weak (Cramer’s V = 0.15). This could be because large amounts of intermittently exposed landbridge habitat remain uninhabited by the species in our panel (Supplementary figure S2). For example, we know that large amounts of exposed landbridge habitat available to our panel’s Papuan species range from lowland evergreen to seasonal rainforest, with an expansion of upland habitats of evergreen to seasonal rainforest. Transitional hill forest formed an ecotone between these two (Cannon et al., 2009). However, most of our species inhabit dense forest (Supplementary table S7), which is unlikely to comprise newly exposed landbride habitat.

Our results also reveal that both passerine and non-passerine island endemics have entered the Holocene with low genetic diversity (Figure 1). A loss of genetic diversity predisposes species to extinction (Mays et al., 2018; Spielman et al., 2004) with a potentially causal relationship for birds (Evans and Sheldon, 2008). The current climate crisis has greatly imperilled biodiversity with the documented loss of many vertebrate species (Steadman and Martin, 2003; Tan et al., 2023; Willis et al., 2004). Fossil evidence has revealed the extinction of several species in the Caribbean since the LGM (Morgan and Woods, 1986; Orihuela et al., 2020; Steadman et al., 1984), and many of the species in our panel are Southeast Asian where rapid habitat loss is ongoing (Sodhi et al., 2004). Climate warming, sea level rise and changes in vegetation are all associated with this loss, and flightless birds and endemics are particularly prone to extinction (Fromm and Meiri, 2021).

Habitat loss, rise in sea level and warming temperatures can rapidly accelerate species extinction, particularly in the tropics (Şekercioğlu et al., 2012). Effective population size, and in extension genetic diversity are associated with survival of a species and its extinction risks (Frankham, 2005). Coalescent effective population size as used in PSMC is an important predictor of species vulnerability to extinction (Brüniche-Olsen et al., 2021; Wilder et al., 2023). Thus, an evolution-informed understanding of species vulnerability to climate change and their associations with paleohabitats can be an important tool to predict species vulnerability to the current climate crisis and future extinction risk (Brüniche-Olsen et al., 2021; Chattopadhyay et al., 2019; Gabrielli et al., 2024; Germain et al., 2023).

A strong, passerine association might be driven by their rapid diversification

Passerines are a hyperdiverse clade representing over 60% of extant avian diversity. The most recent, fossil-calibrated passerine phylogeny (Oliveros et al., 2019) shows that passerines began diversifying in the Middle (47 mya) to Late Eocene (38–39 mya) with increasing diversification rates suggested as we move towards the present. Crown passerines originated in the Australo-Pacific region (Oliveros et al., 2019), and this places the many Papuan and Australian species in our panel close to the passerine diversification centre which is known to have rapid diversification rates (McCullough et al., 2022). Because bird lineages are known to have arisen on tropical islands in the Pleistocene (Andersen et al., 2015; Irestedt et al., 2013; Moyle et al., 2009), it is likely that the Southeast Asian species in our panel represent newly arisen lineages over refugia, possibly belonging to currently rapidly-diversifying clades. This is supported by the fact that several of these species belong to small, often oligotypic families such as Oreoicidae (Aleadryas rufinucha), Ptilonorhynchidae (Amblyornis subalaris), Cnemophilidae (Cnemophilus loriae), Eulacestomatidae (Eulacestoma nigropectus), Ifritidae (Ifrita kowaldi), Machaerirhynchidae (Machaerirhynchus nigripectus), and Rhagologidae (Rhagologus leucostigma). Many of these species belong to poorly studied tropical genera for which species level divergence times are unavailable. The strong passerine association between change in Ne and available habitat area (Cramer’s V = 0.96; figure 3, β = 9.45) reveals that passerines are strongly plastic to the environment, most possibly as a consequence of their ability to rapidly diversify.

Non-passerines representing a paraphyletic outgroup do not show this trend (Cramer’s V = – 0.15). The seven non-passerine species for which both ENM and PSMC analyses were possible belong to six different families (Table 1), representing a diverse sample across the avian phylogeny and precluding non-passerine phylogenetic insight. Notably, we confirm (Cramer’s V = 0.15) that habitat area is positively correlated with genomic diversity in a clade non-specific manner for birds in general (Brüniche-Olsen et al., 2021, 2019).

Species traits influenced historic fluctuations in Ne during the LGP

The interplay of species traits and habitat availability determine how Ne values fluctuate with time in non-endemic birds (Brüniche-Olsen et al., 2021), and our results confirm this for tropical island endemics as well. Habitat change in the LGP was positively associated with Ne fluctuations (Figure 3, β = 9.45), that is, species which showed an increase in habitat in the LGP also showed a concurrent increase in Ne. While diet generalists like omnivores showed no significant association with Ne fluctuations, specialists like frugivores (Figure 3, β = –11.89) and invertivores (Figure 3, β = –11.7) showed significant negative association with fluctuation in Ne, implying that specialist species tended to show decreases in Ne in the LGP. Previous work has shown that large-bodied species have lower values of standing genetic diversity at a given time point (Brüniche-Olsen et al., 2021, 2019; Eo et al., 2011). We find large-bodied species to show a negative association with Ne fluctuations (Figure 3, β = −8.63) as well. This implies that Ne for large-bodied species have decreased in the LGP. Traits like body size and diet determine extinction risk (Ripple et al., 2017; Willis et al., 2004), and our results add to the understanding of how traits modulate Ne changes through time.

Assessing the relative contribution of species traits and habitat availability in determining Ne requires their quantitative measurements. Our methods generated values for available habitat at various time points (Supplementary table S3). However, Ne values can so far only be reconstructed using coalescent methods like PSMC. These values depend upon the parameter settings used. While we find these to be generally concordant across settings (Kruskal-Wallis, p = 0.897, Supplementary table S1), the estimates vary. We therefore chose to cautiously measure only the direction of change (increase or decrease) and not use numerical values.

Caveats of PSMC analyses

Coalescent methods for inferring demographic history like PSMC can trace a local population’s demographic history rather than the entire species history (Gattepaille et al., 2013; Heller et al., 2013; Ptak and Przeworski, 2002). Population structure is thus a confounding factor. However, this is unlikely to be a problem given that all our species are single-island endemics. Large peaks followed by apparent collapses in Ne have also recently been shown to be an artefact of PSMC analyses (Hilgers et al., 2024). We try three different parameter combinations, and look for consensus amongst them in order to try and mitigate this. Because reconstructed Ne values are inferred based on these parameter settings, we do not use them directly for analyses, instead looking at the direction of population size change (increase/decrease/no change). Lastly, we emphasise that Ne is a scaled value of the real population size (Wang et al., 2016), and is measured to be up to one-tenth of the real population size in wild populations (Frankham, 1995). We emphasise caution in interpreting Ne results as a proxy for real population size.

Methods

Species selection

We queried Avibase (Lepage et al., 2014) for all species endemic to tropical islands and selected all species for which assembled whole genome sequences were available on GenBank using short-read data from the Illumina platform. We only selected species that had not been analysed using similar methods already. We excluded genomes assembled from museum specimens. This resulted in a panel of 30 species (Table 1).

Pairwise Sequential Markovian Coalescent (PSMC) analyses

For each assembled diploid genome, we identified the contigs corresponding to the sex chromosomes using Satsuma ver. 2.0 (Grabherr et al., 2010) by aligning contigs to a Gallus gallus (chicken; BioSample: SAMN02981218) genome as a reference. All contigs mapping to the sex chromosomes were removed. Next, we obtained all available Illumina short reads for each species from the Sequence Read Archive (SRA) database of the NCBI, and checked them for errors using FastQC (Andrews, 2010). We used Trimmomatic (Bolger et al., 2014) for preprocessing and trimming the reads. The cleaned reads were aligned to the autosomes using BWA-MEM ver. 0.7.17 (Li, 2013). Next, we used SAMtools ver. 1.10 (Danecek et al., 2021) to merge all the resulting files to generate a single bam file for each species. We further sorted and removed duplicates using SAMtools and estimated the depth of coverage for each species. We excluded species with a depth of coverage < 18X. This is a sufficient value of coverage for PSMC analyses (Li and Durbin, 2011). Next, we implemented the SAMtools mpileup-bcftools pipeline to identify SNPs. The minimum and maximum depths for calling SNPs was set to one-third and double the mean depth respectively following Nadachowska-Brzyska et al., (Nadachowska-Brzyska et al., 2015).

We used three sets of parameters for PSMC analyses and looked for consensus trends amongst them because PSMC analyses can result in spurious peaks of effective population size (Hilgers et al., 2024). We used the following parameter sets for PSMC analyses: –t5 –b –r1 –p “4 + 30 * 2 + 4 + 6 + 10”, –t5 –b –r1 –p “2 + 2 + 30 * 2 + 4 + 6 + 10”, and –t5 –b –r1 –p “1 + 1 + 1 + 1 + 30 * 2 + 4 + 6 + 10”. We performed 30 iterations for parameter optimisation and ran 100 bootstrap replicates to judge the uncertainty in our estimates. Bootstrapping was done by randomly sampling chromosome segments with replacement and running PSMC on them. To plot the results from the PSMC analyses, we used previously estimated values of mutation rates and generation times for passerines (1.4 x 10e–9 years/site and 2 years respectively; (Ellegren et al., 2012; Nadachowska-Brzyska et al., 2016) and non-passerines (1.91 x 10e–9 years/site and 1 year respectively) (Nam et al., 2010) (Supplementary material S1, supplementary table S1).

Species and climate records

We performed reconstructions of paleohabitats through ecological niche models and reconstructed species distributions from four time periods: the Last Interglacial (LIG, approx. 110,000–130,000 years ago), the Last Glacial Maximum (LGM, approx. 20,000 years ago), the Mid-Holocene (MDH, approx. 6000 years ago) and Current (CUR, present day).

We used R (ver. 4.2.1; R Core Team, 2021) for all paleo-habitat modelling analyses. We accessed the location records of 20 bird species from the Global Biodiversity Information Forum (GBIF) in September, 2023 using the “rbgif” package (ver. 3.7.9) (Chamberlain and Boettiger, 2017) in R. We discarded duplicates and retained only human observed records using the “tidyverse” (v-2.0.0) package (Wickham and RStudio, 2023) followed by the ‘CoordinateCleaner’ package (v-3.0.1) (Zizka et al., 2023). We used the “spThin” package (Aiello-Lammens et al., 2014) to account for spatial autocorrelation and finally discarded spurious distribution records like records overlaying water bodies, buildings, roadways, and railways using QGIS (v-3.34+; https://www.qgis.org/).

We used 19 climate variables (labelled BIO1 – BIO19; Supplementary table S2) for the four different time periods. We downloaded CHELSA (https://chelsa-climate.org/; Karger et al., 2021) climate predictors at a resolution of 2.5 arc-minutes (∼5 km) for all regions except for Amazona guildingii from Saint Vincent where we used a resolution of 30 arc-seconds (∼1km). For the latter, we used WORLDCLIM (https://www.worldclim.org/; Fick and Hijmans, 2017) for data for the LIG and MDH while CHELSA was used for the LGM and CUR.

We followed Chattopadhyay et al., (2019) and used a global dataset approach to account for idiosyncratic biases due to smaller datasets to extract the location points of each endemic-island group against each climatic variable. We tested for multicollinearity using variance inflation factor (VIF) using the ‘usdm’ R package (Naimi, 2023) and considered variables for further analyses if their VIF value was ≤ 5.

Habitat Suitability Modelling and Area Calculation

We used R to reconstruct paleo-climatic suitable habitat for 28 endemic island birds through ecological niche models for the four time periods considered. For this, we accessed the Global Biodiversity Information Forum (GBIF) to extract our species’ occurrence records. We could not perform ENM analyses for two species due to a paucity in the number of GBIF occurrence points for them (Table 1).

Data partitioning and model evaluation

For each species, we used available occurrence records from the present day to generate pseudo-absence data points. We generated 500-10,000 pseudo-absence/background points for each species depending upon the area of the island (Supplementary table S3), except for Amazona guildingii endemic to Saint Vincent where we generated 40 points. Saint Vincent’s small land area could not accommodate more points than this. We used a subset of bioclimatic variables to generate a weighted average ensemble species distribution model (eSDM) for each island group. eSDMs were implemented in the R package “sdm” (v-1.2.37; Naimi and Araujo, 2016) applying the ‘MaxEnt’, ‘GLM’, and ‘BRT’ algorithms. eSDMs account for the limitations of different models by generating a weighted average of multiple models (Araújo and New, 2007; Dormann et al., 2018; Naimi and Araujo, 2016). We used k-fold cross-validation (CV) with replication for each method for training and test datasets across each endemic island group (Supplementary table S3). Model accuracy was measured using the Area Under the Curve (AUC) and True Skill Statistics (TSS) metrics. We first performed the above analyses for the present-day (CUR) distribution using the weighted-average of AUC and TSS. AUC and TSS values greater than 0.9 and 0.75 respectively have been shown to be indicators of superior model performance (Ahmad et al., 2019) and we chose a similar threshold (AUC ≥ 0.9 & TSS ≥ 0.8) to get the best ensemble model for most of the species except for 10 species for which data quality was poor (Supplementary table S3). Models for the other time points: the LIG, LGM, and MDH, were generated using the eSDM model generated based on CUR data.

Suitable area analysis

All spatial analyses were carried out using the R package ‘terra’ (ver. 1.7.13) (Hijmans et al., 2024). Equal area projections were used to calculate the absolute suitable area of each species across the four time-periods. To account for a spherical Earth and the limitations of landmass depiction of the original eSDM raster (World Geodetic System 1984) we reprojected it into a pseudocylindrical projection (Eckert IV). Further, we used the average quantile threshold pixels to measure the availability of absolute suitable areas for the LIG, LGM, MDH, and CUR periods (Supplementary table S4).

Statistical Analyses

We checked for associations between Ne and habitat change during the Last Glacial Period (from LIG to LGM) (LGP) using Cramer’s V implemented in the“polycor” package in R (Fox, 2022). Cramver’s V varies between −1 to 1, with values closer to 0 suggesting no-correlation. This was done for the entire dataset, and then for passerines and non-passerines separately.

We further explored the effects of species biology, habitat, and phylogenetic constraint on the fluctuation in Ne during the LGP using Bayesian Multilevel Models (MLMs), using the “brms” package in R (Bürkner, 2018). Our fixed effect predictors were habitat change from the LIG to LGM, diet (invertivore, frugivore, or omnivore), body mass, and clade identity (passerine and non-passerine). The latter allows us to check if belonging to the rapidly diversifying passerine clade results in a signal. We also included country (the island/archipelago a species is endemic to) as a random effect variable for the analysis to account for country-specific effects. For the response variable i.e., the change in Ne, a Bernoulli distribution with a logit link because it is a binary response variable.

We ran a series of models (n=12) to find the best-fit model (brms11, Supplementary table S6) based on leave-one-out cross-validation. We generated a total of 16,000 posterior samples by using the No-U-Turn Sampler (Hamiltonian Monte Carlo) algorithm with 4 chains with 5000 iterations each (with 1000 warm-up iterations; supplementary information S3). Further, we estimated the posterior parameters with 95% confidence intervals to find negative or positive associations between the Ne and its predictors. In addition, we assessed the posterior convergence of the sample through R-hat statistics (R^ = 1) where values close to 1 suggest an ideal convergence. We also reported the effective sample size for bulk and tail distributions (Supplementary table S5).

Data availability

The raw data used for this study is available on NCBI and details of the accession ID are provided in Table S1. All codes used in this study are provided on the Zenodo URL https://doi.org/10.5281/zenodo.14603966.

Acknowledgements

BC acknowledges the support of the Trivedi School of Biosciences. KMG acknowledges the support from the DBT-Ramalingaswami Fellowship (No. BT/HRD/35/02/2006). VI is supported by DBT grant (BT/PR42830/BRB/10/1995/2021) to KMG.

Additional files

Supplementary figures

Supplementary tables