Incorporating the historical effects of Pleistocene sea-level change into models of island biogeography is thus essential to better understanding the patterns of community assembly and gene flow observed in island archipelagos (Lohman et al., 2011). However, while many studies have qualitatively and quantitatively described the impact of transient land bridges in island systems (Cros et al., 2020;Reilly et al., 2019), few software packages exist that allow for such analyses to be performed in a generalised and standardised manner.

To address the need for generalised methods that quantify the effects of sea-level change on islands, we present PleistoDist (https:// github.com/g33k5 p34k/Pleis toDistR), an R package that allows users to visualise sea-level change over time and calculate multiple metrics of island shape, inter-island distance, net inter-island migration and inter-island visibility, all normalised over time. This tool generalises workflows akin to those developed by Norder et al. (2018) and Rijsdijk et al. (2014) by simplifying the effect of sea-level change using a eustatic sea-level reconstruction, which assumes that sealevel change is globally uniform, and applying it to high-resolution bathymetric datasets to reconstruct the geographical extents of islands at different time points. This software is designed to be applicable across a wide range of geographical contexts, and the maps and island metrics generated by this workflow can be used to generate null expectations and hypotheses for ecological and evolutionary analyses. We demonstrate this with a nonexhaustive series of worked examples spanning a variety of taxa and localities, showing how PleistoDist can help address various biogeographic questions and be used in conjunction with other software packages.

PleistoDist requires two main input files: a bathymetry raster to generate maps of island extents at different historical sea levels and a shapefile of source points from which to calculate island metrics and/or pairwise distances. By default, users are advised to use bathymetry data from the General Bathymetric Chart of the Oceans (GEBCO: https://www.gebco.net), a global database of Earth's terrestrial and undersea terrain at 15-arc-second resolution, although PleistoDist can accept any type of ASCII-formatted bathymetry data as input. Users must also specify a map projection appropriate to their area of interest, as well as a time cut-off to set the temporal scope of the analysis.

PleistoDist models historical sea-level change by decomposing eustatic Pleistocene sea levels (by default from Bintanja and van de Wal ( 2008)) into a number of discrete intervals for a user-defined period of time (3000 kya by default, depending on the range of the sea-level reconstruction). These intervals can be calculated in two different ways: by binning over time with the getintervals_time function or by binning over sea level with the getintervals_sealvl function (Figure 1). Binning over time involves dividing the timespan of interest into several equal time intervals and calculating the mean sea level for each interval (Figure 1a). In contrast, binning over sea level involves calculating the total range of sea-level change for the timespan of interest and dividing this range into several equal bins (Figure 1b). For the latter method, the mean sea level of each bin is the average of the minimum and maximum sea level for that bin, without respect to time. Once calculated, the intervals are written to an 'interval file' (Figure 1) in comma-separated value (CSV) format in the output folder.

Users can generate maps of island extents based on mean sea levels specified in the interval file (Figure 1) using the makemaps function. This function reprojects the input bathymetry raster into the user-specified map projection using a bilinear resampling method and generates raster and shapefile outputs of island extents in three formats: a shapefile of island polygons, an ASCII flat raster with no topography and an ASCII topographic raster that preserves the original elevations of each land pixel. Since the makemaps function reads intervals and mean sea levels directly from the interval file, users can customise the output maps generated by this module by manually editing the interval file. Users can also specify an offset value that accounts for constant rates of tectonic uplift or subsidence across the entire modelled region (Ali & Aitchison, 2014;Kealy et al., 2017).

Based on the maps generated, PleistoDist can calculate a variety of intra-and interisland metrics for each time or sea-level interval specified in the interval file, as well as calculate the weighted mean of these values over time. For metrics of island shape, PleistoDist includes functions that calculate island perimeter, area and surface area for each interval. As for interisland distances, PleistoDist can calculate three kinds of island-to-island distances (centroidto-centroid, least shore-to-shore and mean shore-to-shore distances; Figure 2), as well as two kinds of point-to-point distances (Euclidean and least-cost distances; Figure 3) for each interval.

Islands that merge at particular intervals are considered fused for these metrics.

In addition to metrics of island shape and interisland distance, PleistoDist includes two higher level functions that calculate the expected equilibrium net migration between island pairs based on a model described by MacArthur and Wilson (1967), as well as the visibility of an island relative to an observer on another based on an estimate of the maximum sightline (i.e. horizon distance) of the observer (Figure 4). Details on how the calculations are performed can be found in the Supporting Information. Papadopoulou and Knowles (2015) studied the effects of Pleistocene island connectivity on the genetic differentiation patterns of Amphiacusta sanctaecrucis populations in the Virgin Islands, and found that population divergence times broadly correlate with a period of fluctuating sea levels and interisland land-bridge connections (~75-115 kya). In addition, all Virgin Island populations appeared to exhibit a pattern of isolation-by-distance, except those from St. Croix.

Because these analyses were based on present-day Euclidean distances, we used PleistoDist to replicate the landscape genetic analyses using geographical distance matrices that account for sea-level change over time, for a time cut-off of 20 kya (corresponding to the last glacial maximum [LGM]) and 20 time intervals. Similar to Papadopoulou and Knowles (2015), Mantel tests show no significant correlation between genetic and Euclidean distance (R 2 = 0.228, p = 0.0521; Table 1) when all populations are included in the analysis (Table 1). None of the other geographic distance metrics showed a significant Mantel correlation with genetic distance either. However, given the susceptibility of Mantel tests to spatial autocorrelation (Legendre et al., 2015), it is probable that genetic variation is nonlinearly distributed across geographic space. Using distance-based redundancy analyses (dbRDA), which account for spatial autocorrelation, we found that unlike Papadopoulou and Knowles (2015), we observed a highly significant correlation between genetic distance and timecorrected least shore-to-shore and centroid-to-centroid distances (p = 0.00661 and 0.00704, respectively; Table 1), while Euclidean distance showed a marginally significant correlation with genetic distance (p = 0.0348; Table 1). Further model selection indicated that the time-corrected centroid-to-centroid distance best fit the variation observed in the genetic data, suggesting that genetic divergence patterns were likely structured by broad-scale isolationby-distance driven by interisland overwater dispersal between panmictic populations. Using PleistoDist, we reanalysed angiosperm species-area relationships on the true islands of the Aegean Sea using both contemporary and time-corrected estimates of island area for 40 depth intervals over the last 20 kya. We found that species richness is significantly positively correlated with island area for all angiosperm chorotypes (Table 2), consistent with the findings of Hammoud et al. (2021). However, unlike Hammoud et al. (2021), our extended analyses suggest that present-day 2D island area is the best predictor of angiosperm species richness instead of time-corrected 2D area (Table 2). Our re-evaluation suggests that species turnover rates on true islands may be fast, resulting in a rapid equilibration of species richness to contemporary island extents.

One of the predictions of the theory of island biogeography (MacArthur & Wilson, 1967) is that the relative rate of migration between island pairs can be predicted by the size, relative orientation and distance of source islands. Gyllenhaal et al. (2020) tested this F I G U R E 2 PleistoDist calculates three interisland distances: centroid-to-centroid distance, least shore-to-shore distance and mean shoreto-shore distance, illustrated here with Viti Levu (left) and Gau islands, Fiji. Note how interisland distances are asymmetric for the mean shore-to-shore distance.

TAN et al.

expectation by estimating the rates of interisland migration in the Fiji bush warbler Horornis ruficapilla between the four large islands of Fiji, concluding that rates of interisland migration are largely consistent with neutral expectations. This neutral expectation can now be easily and accurately calculated using PleistoDist while accounting for the effect of sea-level change over time.

Re-running this analysis using time cut-off of 115 kya, corresponding with the start of the last glacial period, for 40 sealevel depth intervals, we found that empirical ratios of migrants for Kadavu → Viti Levu Viti Levu → Kadavu broadly correspond with neutral expectations, suggesting that net migration rates are unlikely to have changed much over the last 115,000 years (Figure 5). In contrast, empirical F I G U R E 4 PleistoDist estimates the visibility of a destination island relative to an observer on an origin island by calculating the horizon distance and performing a viewshed analysis to estimate the visible nonoccluded area of the destination island, illustrated here with Ovalau and Gau islands, Fiji. migration ratios for Taveuni → Vanua Levu Vanua Levu → Taveuni seem consistent with neutral expectations based on present-day sea levels but not expectations averaged over the last 115,000 years (Figure 5). The empirical net migration ratios for Viti Levu → Vanua Levu Vanua Levu → Viti Levu are challenging to interpret given the broad confidence intervals of the empirical data, but the mean empirical migration ratio from the model removing a putatively admixed individual is similar to the value inferred from the centroid-to-centroid distance model averaged over 115,000 years (Figure 5).

To demonstrate the utility of PleistoDist for spatially explicit simulations, we ran forward-in-time, spatially explicit Wright-Fisher simulations in SliM3 (Haller & Messer, 2019) across the Samoan archipelago using maps generated by PleistoDist, with a metapopulation size of 20,000, for 50,000 generations. We ran simulations both with and without sea-level change for eight mean dispersal distances, for 100 replicates per sea-level regime/dispersal distance combination (full model parameters in Supporting Information). As expected, F ST was negatively correlated with dispersal distance (Figure 6), and π was positively correlated. The impact of accounting for change in sea level was most noticeable for islands that were much closer at glacial maxima (i.e. Savai'i and Upolu islands, Samoa), but not when dispersal distance was high.

To demonstrate the utility of PleistoDist for modelling simple archipelagic colonisations, we simulated the colonisation of the Solomon Archipelago using the SLiM GUI (Haller & Messer, 2019), with initial colonists arriving on Buka/Bougainville islands (Mayr & Diamond, 2001). To simulate an expanding colonising population, we used a non-Wright-Fisher simulation, starting with 50 individuals on the easternmost large island in the archipelago. We used a mean long-distance dispersal distance of ~46.7 km (100 pixels), an exponential long-distance dispersal kernel, and offspring long-distance dispersal probability of 1%. For all simulations (two shown at https://github.com/g33k5 p34k/Pleis toDis tR/tree/main/simul ation_ scripts), Makira Island was colonised last due to its relative isolation and location at the far eastern end of the archipelago, demonstrating the utility of the script in generating null expectations against which complex empirical patterns can be compared.

To demonstrate the utility of accounting for migration rate variation between nonconnected populations, we used msprime v1.0.0 (Baumdicker et al., 2022) to model dispersal between Viti Levu and Kadavu islands, Fiji. We used PleistoDist to estimate per interval migration rates using least shore-to-shore distance (Figure 2) and an exponential dispersal kernel with means of 0 (no dispersal), 10, 20, 40 and 80 km. We calculated population size by multiplying island area for each interval by a density of 10 individuals/km 2 . We either held migration rates and population sizes constant (based on current sea level) or changed them based on sea level in intervals of 5000 years (2500 generations). The initial population split was set to 200,000 generations. At the end of the simulation, we calculated F ST and per-island π for each sea-level regime. As with the SLiM simulations, the impact of sea-level change was more notable when dispersal was low or not present (Figure 7). These simulations also demonstrate how sea-level change and interisland dispersal can benefit π on less isolated small island populations (i.e. Kadavu).

Despite its broad utility, PleistoDist nonetheless relies on a simplification of Pleistocene sea-level change, and users should be aware of its limitations. One major assumption is that bathymetry remains constant, or changes at constant rates, throughout the time-scale of the analysis, and PleistoDist therefore does not account for the effect of localised tectonic activity or transportbased processes such as sedimentation and erosion. For example, PleistoDist analyses of Taveuni island in Fiji with a cut-off time close to and exceeding 700 kya are unlikely to be accurate since the island likely emerged around that time (Cronin & Neall, 2001).