Research Article
Print
Research Article
Temporal Altitudinal Biogeographic Shifts (tabs): R package for reconstructing biogeographic shifts in terrestrial and marine systems over time
expand article infoJohannes De Groeve‡, Kenneth F. Rijsdijk‡, Eline S. Rentier§, Suzette G. A. Flantua§, Sietze J. Norder|
‡ University of Amsterdam, Amsterdam, Netherlands
§ University of Bergen, Bergen, Norway
| Utrecht University, Utrecht, Netherlands

Corresponding author: Johannes De Groeve ( j.degroeve@uva.nl )

Academic editor: Robert Whittaker
© 2025 Johannes De Groeve, Kenneth F. Rijsdijk, Eline S. Rentier, Suzette G. A. Flantua, Sietze J. Norder.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation: De Groeve J, Rijsdijk KF, Rentier ES, Flantua SGA, Norder SJ (2025) Temporal Altitudinal Biogeographic Shifts (tabs): R package for reconstructing biogeographic shifts in terrestrial and marine systems over time. Frontiers of Biogeography 18: e151677. https://doi.org/10.21425/fob.18.151677
Open Access

Abstract

Past climate fluctuations have profoundly influenced the global distribution of biomes, ecosystems and species, altering their altitude, spatial configuration, area, and connectivity. Notable examples include island archipelagos and alpine biomes, where shifts in sea levels and forest lines respectively reshaped their spatial structures. To understand how such changes affected species distributions and biodiversity patterns, we require spatially explicit reconstructions over continuous time series. However, a comprehensive and reproducible methodology that captures their spatio-temporal dynamism is lacking. Here, we introduce the R package Temporal Altitudinal Biogeographic Shifts (tabs), a tool designed for reconstructing spatial configurations over time, focusing on biogeographic systems bounded by an altitudinal range. We demonstrate the use of tabs by modelling spatial configurations of island archipelagos and alpine biomes in response to fluctuations in sea levels and forest lines. Unique to tabs, it can also account for geophysical effects on regional sea levels, and for geotectonic topographic changes. Beyond past reconstructions, tabs can project spatial configurations shaped by future climatic conditions. This versatile package is easily adaptable to various altitude-bounded biogeographic systems influenced by long-term climatic variations, such as coral reefs and shelf seas. Studying the shifts in biogeographic systems through continuous spatial reconstructions, rather than snapshots in time such as the Last Glacial Maximum, captures the nuances of continuously changing environments and provides a more complete understanding of the biogeography of our planet.

Highlights

  • tabs R package enables detailed spatio-temporal reconstructions of biogeographic systems along altitudinal gradients.

  • Designed with user-friendliness in mind, it offers a streamlined setup, enabling researchers to start reconstructions quickly without extensive technical expertise.

  • tabs is applicable to diverse biogeographic systems, including islands, coastal zones, and alpine biomes.

  • It provides flexibility to tailor input data, parameters, and settings according to specific local conditions and research objectives.

  • tabs generates outputs using a standardised data structure and open data formats, enabling effective analysis and presentation of results.

Keywords

Climate-driven biome shifts, ecosystem modelling, paleogeographic reconstructions, Quaternary climate change, sea level fluctuations, spatio-temporal analysis

Introduction

Past climate fluctuations have been a significant driver of evolutionary and ecological processes, profoundly influencing the spatial distribution of biogeographic systems, biomes, ecosystems and their species communities (Rangel et al. 2018; Hewitt 2000; Xu et al. 2023). Climatic variations have repeatedly altered coastlines through sea-level fluctuations, leading to changes in the area and connectivity of terrestrial and marine systems (Fernández-Palacios et al. 2016; Kealy et al. 2017, 2018; Salles et al. 2021). These fluctuations have not only transformed land-sea boundaries (Ali and Aitchison 2014; Weigelt et al. 2016; Norder et al. 2019), but also reshaped the distribution of biomes and ecosystems across marine, coastal and mountainous regions (Sandel et al. 2011; Camoin and Webster 2015; Svenning et al. 2015; Flantua et al. 2019, 2020). In particular, climate-driven biogeographic shifts in mountainous areas, where species respond to warming by moving upslope and to cooling by moving downslope, highlight the dynamic relationship between climate and biome and ecosystem distribution (Rull and Nogué 2007; Dirnböck et al. 2011; Flantua et al. 2014; Flantua and Hooghiemstra 2018). Understanding these shifts is crucial for understanding ecological responses to climate fluctuations in the past, present, and future (Williams et al. 2007; Fernández-Palacios et al. 2016; Brown et al. 2020).

By examining how biogeographic systems shifted their distributions over time, we can reveal dynamics of expansion and contraction, isolation and connectivity, as well as their influence on ecological patterns and processes (Flantua et al. 2019, 2020; Norder et al. 2019). For instance, biogeographic shifts driven by Quaternary sea-level fluctuations and geological dynamics have shaped species distributions (Whittaker et al. 2008; Ali and Aitchison 2014; Rijsdijk et al. 2014; Fernández-Palacios et al. 2016), migrations (Comes and Kadereit 1998), and gene flows (Papadopoulou and Knowles 2015, 2017). Therefore, quantifying changes in the extent, configuration, and connectivity of biogeographic systems over time, is essential to enable a deeper understanding of how past environments influenced present-day biodiversity, and to provide better predictions of the effect of isolation and connectivity in response to future climate change (McGuire et al. 2016; Brown et al. 2020).

Despite a broad interest in spatio-temporal reconstructions and projections of biogeographic systems, there is currently no standardised, adaptable workflow capable of reconstructing spatio-temporal shifts across diverse biogeographic contexts, including both marine and terrestrial systems. Workflows are often developed to generate reconstructions for a specific biogeographic system, spatial or temporal scale, or research question. This limits their utility to other researchers who may wish to incorporate not only climate-driven shifts but also geological dynamics, or those who are focused on projected instead of reconstructed biogeographic shifts. Finally, current workflows often fail to provide outputs in open, compatible formats that facilitate further analysis, visualisation and collaboration. Thus, applying a proposed workflow also tends to require extensive data preparation and technical expertise, further limiting their usability. To address these gaps and advance our understanding of biogeographic dynamics across both space and time, the development of a versatile, user-friendly workflow is essential.

Here, we introduce a globally applicable, reproducible workflow for reconstructing and projecting biogeographic shifts in terrestrial and marine systems in response to climate, topography, and geological changes. The Temporal Altitudinal Biogeographic Shifts (tabs) R package models spatial configurations over time, focusing on biogeographic systems bounded by an altitudinal range, such as island archipelagos, alpine biomes, coral reefs and shelf seas. Additionally, to derive regional sea-level estimates, tabs allows to account for crustal deformations and gravitational forces generated by ice sheets, often neglected in earlier studies (e.g. Rijsdijk et al. 2014; Norder et al. 2018; Tan et al. 2023). It also allows for topographic corrections based on local geological conditions, including uplift and subsidence. tabs is highly versatile, enabling easy modification of input data and parameters to suit local conditions. It produces rapid, accurate reconstructions with minimal manual input, facilitating the modelling and quantification of spatial configuration changes across time. tabs generates outputs that are analytically robust and well-suited for further analysis and publication. We illustrate tabs by reconstructing the spatial configurations in two biogeographic systems, island archipelagos and alpine biomes, in response to fluctuations in sea level and forest lines. In the Methods, we first describe the main workflow of tabs, including its capabilities, functions and arguments. Its application is further demonstrated through the two case studies presented in the Results.

Methods: tabs workflow

Functionality

The functionality of tabs builds on and automates the foundational work laid out by Norder et al. (2018), who developed an R based workflow to produce a global island and archipelago configuration database. The study combined a global eustatic sea-level curve (e.g., cutler, lambeck; see section ‘curve parameter‘) with a topographic-bathymetric model to reconstruct island area change and configurations, including labelling of individual present and past island polygons. While the workflow from Norder et al. (2018) required manual re-runs to incorporate new datasets, tabs introduces a dynamic, customisable framework. It allows iterative reconstructions that can be fine-tuned based on varying climate and geological conditions, representing a major leap forward in modelling biogeographic dynamics. tabs further builds on early work of Simaiakis et al. (2017) who presented a reconstruction of island change and configurations for the Ionian and Aegean sea extents. Their model tracks island dynamics over the past 26,000 years, but different from Norder et al. (2018), this study accounts for regional variations in sea-level changes due to crustal deformations and gravitational forces generated by ice sheets. A global raster of sea level changes is available in tabs and allows to account for regional variations in sea levels induced by these geophysical processes. tabs also provides the tools to account for geotectonic topographic changes. For instance, the Ionian and Aegean seas exhibit regional variations in tectonic crustal uplift and subsidence ranging from 1 mm/yr to -1 mm/yr, which may thus result in an absolute topographic change of 26 meters in 26,000 years (Ganas and Parsons 2009). The integration of such complex reconstruction efforts in tabs demonstrates its capacity to model highly complex biogeographic systems with high temporal and spatial resolution. tabs further integrates previous work reconstructing past spatio-temporal dynamics of a high altitude alpine biome (Flantua and Hooghiemstra 2018; Flantua et al. 2019). In these studies, the altitudinal shifts and resulting connectivity dynamics of the high Andean páramo were reconstructed in response to temperature changes during the Quaternary period, mirroring sea-level driven island dynamics. tabs also supports applications for marine biogeographic systems, including the reconstruction of potential past coral reef and shelf sea distributions since the Last Glacial Maximum (De Groeve et al. 2022a). Using a comprehensive set of global datasets and parameter settings, tabs enables the creation of spatial-explicit models representing biogeographic systems worldwide. The next paragraphs describe the main datasets and function parameters used by the core function of tabs:reconstruct().

Global datasets

tabs allows for generating reconstructions using region-specific custom data sources, but various built-in datasets can be used to create reconstructions anywhere on the globe. To access these datasets, first-time users are encouraged to run the setup() function, which creates a directory structure in a specified or default local path, and automatically downloads these datasets (requiring 15 GB of disk space). Datasets include: i) a tiled version of the global bathymetry and topography (GEBCO Bathymetric Compilation 2024), ii) a spatio-temporal sea-level curve, iii) labelling datasets, including the Global Shoreline Vector (GSV, Sayre et al. 2019), bioclimatic and physical characterisation of the worlds’ islands (Weigelt et al. 2013a, 2013b), tectonic plates and orogens (Bird 2003), an open source labelling points dataset for many types of features (GeoNames 2023) and the Global Mountain Biodiversity Assessment (GMBA) Mountain Inventory v2 (Snethlage et al. 2022). These datasets are made available via Figshare to facilitate and manage any future changes (De Groeve et al. 2022b, 2023, 2024). For instance, the GEBCO model is updated on an annual basis and we aim to upload and announce the integration of most recent versions with upcoming tabs version releases. For more details on dataset sources and examples, see the vignettes on the tabs website (De Groeve et al. 2025c).

reconstruct() function

The reconstruct() function is the core function of tabs and operates with several integrated helper or utility functions (Fig. 1). The purpose of this function is to reconstruct and label past or future changes in terrestrial or marine biogeographic systems for a specific region in response to changes in climate and geology. These reconstructions are generated by intersecting a bathymetric or topographic model with a (spatio-temporal) curve and optionally correcting it for regional topographic changes. The resulting reconstructed spatial configurations (hereafter “shapes”) of the biogeographic system of interest are optionally named using a labelling dataset. The reconstruct() function has six main input parameters (region, reclabs, topo, curve, iso, correction) and several additional parameters (buffer, noise(rm), aggregate, fact, units, fillholes, metrics) that can further fine-tune reconstructions in accordance with the users’ needs.

region parameter:

To generate a tabs reconstruction, it is important to first define and select the region of interest. This can be done using one of five approaches (Table 1): (i) interactive drawing of an extent, (ii) extent coordinates, (iii) path to a local spatial dataset, (iv) a sf or spatVector R-object, and (v) region-name definitions including islands, archipelagos, countries, tectonic plates and mountains. In total 32,255 region name-definitions are available in tabs and can be explored using the built-in regions dataset. For more details on region selection, see the vignettes on the tabs website (De Groeve et al. 2025c).

Figure 1. 

The infographic provides a summary of the tabs R package and the steps incapsulated by the function reconstruct(). (A) A region of interest on the globe can be selected using a wide range of methods, including by extent, filename, R objects, or geographic names such as islands, countries, archipelagos, mountains, and tectonic plates. (B) For the selected region, a harmonised input dataset compilation is prepared, which includes topo (I1; topography), curve (I2), correction (I3), labs (I4; labels) and iso (I5; altitudinal offset). Each dataset can be customised, and both curve and correction support various formats, including numeric vector(s) or raster(s). The latter two datasets are intersected with the region of interest and resampled to the resolution of the topography. If needed, the correction dataset is also aligned to the temporal resolution of the curve. The labs dataset, specified through the reclabs argument, can be one of the default global datasets (e.g., islands, mountains, geonames), a custom column name, or labelling can be ignored altogether. The iso variable defines the altitudinal offset from the curve for a biogeographic system (e.g., coral reefs between 0–30 m below sea level), and in this example, it is set to 0. (C) Reconstructed shapes are calculated as topographic pixels higher than the sum of the curve and correction. A line graph illustrates the interaction between these three variables (topo, curve, correction) for a single pixel (cross-marked), showing reconstructed shapes for three timestamps (t0, t1, t2). At t0, the pixel remains below the combined curve-correction threshold, but exceeds it at t1 and t2, resulting in reconstructed shapes. (D/E) The generated shapes are then labelled by identifying the highest point in the reconstructed shapes and intersecting them with the labelling dataset. Emerging shapes (e.g., at timestamp t2) are named using a standard convention (e.g., ‘S’ for shape, time of disappearance, and an identifier). Each shape is also assigned spatial attributes, such as the period, curve and iso, the area (in m² and raster cells), [x, y, z] coordinates, a list of intersecting present-day labels from the labs dataset, and a record of all reconstructed shapes intersecting with the shape at timestamp tx. (E) The reconstruct() function results in three output datasets: reconstructed shapes in vector and raster formats (recvect, recrast), and reconstructed area estimates for each shape (recarea). (F) The final dataset collection includes all the necessary datasets to reproduce the workflow, including all in-and outputs. The dataset collection can be exported, imported, and explored visually using various open data formats, including a (zipped) directory and R object.

Table 1.
Download as
CSV
XLSX

Descriptions and values for the region argument.

Description Value
Default; Interactive drawing of an extent or polygon within R session region=NULL
Specifying the extent as vector or ext-object region=c(xmin,xmax,ymin,ymax)
A path to a locally stored spatial dataset (e.g., shapefile, geopackage) region=’filename’
A spatial (sf, SpatVector) object containing points or shapes region=’spatvector’
A name of an island, archipelago, country, tectonic plate or mountain region=’Greece’
region=’Sporades’

reclabs parameter:

tabs facilitates labelling of the biogeographic shapes, i.e., naming of each reconstructed shape. To assign such labels, the user can specify the reclabs parameter, which includes three approaches (Table 2). Specifically, (i) labelling using the islands, mountains or geonames global datasets (see 2.2. global datasets), (ii) custom labelling and (iii) automated labelling. With the parameter aggregate (boolean; TRUE/FALSE), the spatial scale at which to perform labelling can be defined, specifically for the whole region or for individual reconstructed shapes. For more information on the automatic assigning of labels, see the vignettes on the tabs website (De Groeve et al. 2025c).

Table 2.
Download as
CSV
XLSX

Descriptions and values for the reclabs argument.

Description Value
Labelling using the islands, mountains or geonames global datasets reclabs=‘isls’
reclabs=‘mnts’
reclabs=‘peak’
Custom labelling by specifying the column name of a user-defined region-object reclabs=‘name’
Automatically generated labels reclabs=FALSE
Default; No labelling dataset defined; Labelling depends on the input of the region-object reclabs=NULL

topo parameter:

By default, the GEBCO raster (2024) is used as the topographic (topo) object and consists of both bathymetry and topography at a global scale at a resolution of 30 arc seconds. Alternatively, the user can specify a regional digital elevation model as a path (e.g., topo=’path_to_file.tif’). tabs automatically extracts the highest point identified in the topographic object which is subsequently used for labelling reconstructed shapes. The spatial resolution of the topo object will be assigned to the tabs outputs (as produced by the reconstruct() function), unless the resolution factor (fact) is modified. Outputs are assigned the projection system of the input topo object.

curve parameter:

The curve can be defined as any dataset which expresses the relative shift in altitude over time for the biogeographic system of interest. Curves can vary both in space and in time within a region, thus it can be a vector of one or more altitude values, or one or more rasters. A typical example of a curve is a sea-level curve, expressing the change in sea level compared to the present-day reference (i.e., 0 m below sea level). Another example is the upper forest line (UFL), expressing the altitude above present-day sea level at which continuous forest persists. When a curve is represented as a raster, local variations can be incorporated. For instance, spatial differences in sea level due to gravitational pulling by large ice masses, or regional shifts in the upper forest line driven by slope, aspect and temperature gradients across mountain latitudes. A curve can also be defined by a single value, for instance, the median sea level across the Pleistocene or a future scenario of sea-level rise. Hence, tabs allows a flexible and customisable definition of a curve. tabs has several built-in curves that can be called by their name (Table 3). These curves represent different global sea-level curves (past and future), and an example of a UFL curve for the Northern Andes (funza; Torres et al. 2005; Flantua et al. 2019). Note that the shift in altitude should use the same unit as the topo object, which is usually expressed in meters. For more details, see the vignettes on the tabs website (De Groeve et al. 2025c).

Table 3.
Download as
CSV
XLSX

Curves incorporated in tabs.

Name Description Time span Temporal resolution Spatial resolution Source
Lambeck Global eustatic sea-level curve 35 ky BP 1000 yr Lambeck et al. 2014
curve=’lambeck’
Cutler Global eustatic sea-level curve 120 ky BP 1000 yr Cutler et al. 2003
curve=’cutler’
Bintanja Global eustatic sea-level curve 3 Myr BP 1000 yr Bintanja and van de Wal 2008
curve=’bintanja’
De Groeve Global spatio-temporal sea-level curve 26 ky BP 500 yr 0.2° De Groeve et al. 2022a
curve=’st_curve’
IPCC Global average sea level for four IPCC scenarios (SSP1, SSP2, SSP3, SSP5). Present to 2100 20 yr Intergovernmental Panel on Climate Change (IPCC) 2023
curve=’IPCC’
Funza Upper forest line altitudes derived from fossil pollen record Funza (Colombia) 1 Myr BP – 22 ky BP 1000 ky Torres et al. 2005, Flantua et al. 2019
curve=’funza’
User-defined A vector with the altitude(s) or a data frame containing altitudes at different time steps. Any Any Any
NULL If no argument is provided, construction will be produced for the present

iso parameter:

The iso (hypse) argument defines the offset from the curve and covers the altitudinal range of the biogeographic system to reconstruct. By default, this argument is set to 0 meter to reconstruct configurations of coastlines and forest lines, but any altitudinal range that describes the altitudinal distribution of a biogeographic region can be used. For example, the availability of coral and shelf seas varies synchronously with sea-level change, thus to reconstruct their spatial configurations we need to specify the minimum and maximum altitude range at which they may occur, respectively 0 m – 30 m below sea level (BSL) and 0 m – 140 m BSL.

correction parameter:

The reconstruction of biogeographic systems and their spatial configurations over time is influenced not only by specific climatic conditions but also by geotectonic and geophysical processes that can result in topographic changes during the study period. For instance, whether an island is above or below sea level will not only depend on the redistribution of water from ice sheets, but also local effects that change the topography. Examples are uplift and subsidence in tectonic active areas, deposition and erosion of sediment in fluvial deltas and land loss or land gain due to catastrophes such as volcanic eruptions (e.g. Rijsdijk et al. 2020). To account for these types of local topographic changes, the correction parameter can incorporate (1) temporally constant rates of change (e.g., linear uplift or subsidence in mm/yr), or (2) time-varying changes that differ across periods. Constant rates can be provided as a single value or raster, while time-varying changes are defined as a vector of values or a list of rasters. Although a database with global correction rates or values does not yet exist, user-defined corrections can be provided for specific regions. As an example, tabs has one built-in regional correction dataset for uplift and subsidence rates in the Sporades. For more details, see the vignettes on the tabs website (De Groeve et al. 2025c).

Other parameters

In addition to these six main input parameters of reconstruct(), there are other settings that can be tailored to specific purposes. A buffer can be defined to increase the extent of the region definition. With fillholes, remaining holes in reconstructed shapes are removed. Further, it is possible to identify and remove reconstructed shapes below a specific size (with the noise and noiserm parameters respectively) and to generate outputs at lower spatial resolution (fact). Both noise(rm) and fact will reduce computation time, since the former will decrease the number of spatial features, and the latter reduces the number of output cells in rasters. The metrics parameter specifies which metrics are by default computed, currently only including the area change over time. Moreover, several steps of the reconstruct() function can be run separately using specific helper-functions, including get_region(), get_curve() and get_correction().

tabs outputs

The two main outputs of the reconstruct() function are 1) the reconstructed biogeographic shapes at each time step in raster (recrast) and vector (recvect) format, and 2) a table with the calculated surface area for each reconstructed shape (recarea) per time step and a unique ID assigned to each shape. These IDs enable tracking which present-day shapes were merged into larger shapes across time steps. In addition, all preprocessed input datasets are also provided for the defined region (topo, labs, curve, correction). The reconstruct() output can be exported in three file formats: (1) directory, where rasters are exported as GEOTIFF and vectors in GeoPackage, (2) zip, a zipped version of directory-type export and (3) rds/qs2, as the default R export data format. Outputs can be exported by specifying the argument filename, including the path and optional extension in the reconstruct() function or in the export() function. Finally, the import() function allows loading a saved reconstruction (e.g. from a directory or rds/qs2) as an R object of class ‘tabs’.

In addition to these in- and outputs that are stored locally, an R object of class ‘tabs’ will be generated that can be used to explore the outputs interactively with the explore() function. This function will generate interactive maps to visually explore biogeographic shapes over time, or to compare reconstructed shapes to their present-day distribution.

Results: examples of case studies

To demonstrate the application of tabs, we present two case studies that reflect different biogeographic contexts: a coastal system in the Sporades (Greece) and a high-elevation system in the Venezuelan Andes. These examples illustrate how tabs can be used to reconstruct spatial configurations under different altitudinal scenarios. We provide sample code demonstrating the use of the reconstruct() function to illustrate the core of the functionality of tabs. Although sample code for the reconstruct() function is provided, generating the final figures requires additional steps and relies on code and data that are openly available at De Groeve et al. (2025a).

Case study 1: Coastlines - Sporades (Greece)

As a first case study, we reconstructed the coastlines of the Sporades archipelago in Greece since the LGM (26 ky BP) to the present-day using the default input datasets (topo, labs, curve). In addition, we defined a correction raster to account for local changes in uplift and subsidence (sample dataset in tabs). Such a complex reconstruction can be run simply by defining four parameters: (1) region, as the extent of the correction, (2) reclabs, to specify the labelling dataset, (3) correction, to account for topographic changes over the last 26 ky BP and (4) curve, set to the spatio-temporal sea-level curve by De Groeve et al. (2022; Table 2). Note that inputs topo and labs are also available as built-in datasets for the extent of the Sporades, but if the default datasets have been downloaded using setup() it is not necessary to specify the arguments.

# load built-in correction raster

correction <- sporades()$correction

# reconstruct

sporades <- reconstruct(region = ext(correction),

reclabs = ’isls’,

correction = correction,

curve = ’st_curve’)

Using the reconstruct outputs of tabs, we can create, for example, coastline change maps or graphs representing area change and island connectivity since the LGM due to sea-level fluctuations (Fig. 2). Two island groups remained disconnected from each other since the LGM (Fig. 2A), here labelled after the highest island of each group, namely Skopelos (Fig. 2B) and Skantzoura (Fig. 2C). The map shows that most islands within the Sporades archipelago were connected between 25–18 ky BP. More specifically, 12 of the 18 present-day islands were interconnected 24 ky ago, forming the Skopelos island group, while the other six present-day islands formed the Skantzoura island group 14 ky ago. In addition, four paleo-islands (non-existing islands in present) emerged within the Sporades, of which two merged with the Skopelos island group (S-BP0008000-9, S-BP0019000-1), one with the Skantzoura island group (S-BP0013500-3) and one emerged as a separate island 19.5 ky ago (S-BP0019500-3).

Figure 2. 

Reconstruction using the built-in dataset of the Sporades archipelago (Greece). (A) Coastline change map. (B/C) Area change of the island groups Skopelos (B) and Skantzoura (C). Lines indicate the island areas (y-axis, km2) at different moments in time (x-axis, yr BP). Separating lines indicate different islands as connections are lost due to rising sea levels. (D) Sea-level fluctuations at different moments in time (x-axis, yr BP) with and without topographic corrections. The colour gradient in all figures indicates time, where green represents the LGM and white the present.

Case study 2: Mountains – Andes, Venezuela

In our second case study, we reconstruct the area change of the ‘páramo’, i.e., high altitude alpine biome in Venezuela (Northern Andes), during the last 1 million years following Flantua et al. (2019). The authors used the reconstructed altitudes of the UFL, based on the long fossil pollen record Funza (Torres et al. 2005), to delimit the lower boundary of the páramo through time. For the implementation of tabs, we used reconstruct() without labelling and correction raster, while using the default topo. Instead of using the entire curve, we calculated decile statistics, representing the proportion of time that the UFL is higher than a given threshold using proportion steps of 0.1. Here, each proportion of 0.1 indicates that the altitude is above a given UFL-value for 10% of the time (approx. 100 ky). For 10 decile-curve values summing to 100%, we ran the reconstruct() function that subsequently provides the data outputs necessary to generate the decile map, the curve, and area change in relation to the forest line deciles (Fig. 3A, B, C). In addition, we also ran the reconstruct() function from 130–29 ky BP to map the area change over this period (Fig. 3D) (Note: Funza does not cover the last 29 ky BP).

# define region using extent coordinates

region <- terra::ext(c(-72.5,-69.2,7.5,10.2))

# curve

curve <- funza[2:length(funza)]

# calculate forest line deciles from curve

deciles <- quantile(curve, probs = seq(0, 1, by = 0.1))

# reconstruct forest line deciles

andes <- reconstruct(region = region,

curve = deciles,

reclabs = FALSE,

aggregate = TRUE)

# reconstruct 0-130 ky BP

andes130 <- reconstruct(region = region,

curve = funza[1:(130-27)],

reclabs = FALSE,

aggregate = TRUE)

We can see that the present-day UFL altitude (here standardised to UFL = 3250 m a.s.l.) is within the upper decile (UFL = 3105–3404 m a.s.l.), resulting in one of the smallest areal extensions (i.e., 2203 km2) of the páramo in comparison to the last 1 million years. More precisely, only for c. 12 ky, or 1.2% of the time, the UFL was higher and the corresponding area smaller than the present-day. Half of the time (500 ky; middle decile), the UFL was below 2750 m a.s.l. resulting in the páramo to be 1.84 times larger (> 4044 km2) than the present. For 10% of the time (100 ky; lower decile), the UFL was as low as 2213 m a.s.l., resulting in an areal extension c. 3.4 times larger (>7510 km2) than the present-day (2203 km2 for a UFL of 3250 m a.s.l.).

Figure 3. 

Reconstruction of upper forest line (UFL) deciles, representing the proportion of time during the last 1 million years that the UFL was above an altitudinal threshold. (A) The map represents the proportion of time (1000–29 ky BP) that a pixel was above the UFL threshold. (B) Graph shows the altitudinal variation of the Funza UFL from 1000–29 ky ago with the coloured deciles on the background showing the UFL altitudinal range of each decile. (C) Area change in relation to UFL deciles, where the dots show the number of individual “páramo sky islands” at each decile. (D) Area change curve for the last 130–29 ky BP.

Discussion and future directions

Here we present tabs as a workhorse for modelling and quantifying spatio-temporal dynamics of biogeographic systems in response to climate and geology. It is built in such a way that it can be used for any biogeographic system that can be roughly defined by an altitudinal range. Specifically, we modelled shorelines of the Sporades archipelago (Greece) and the upper forest lines delimiting the lower boundary of the páramo (high altitude alpine biome of the Northern Andes) as case studies to showcase tabs’ functionality. These case studies illustrate that tabs requires very limited lines of code and a limited number of parameter settings to reconstruct biogeographic shapes over time. While it is easy to generate such reconstructions using default settings, we showed that tabs is also highly flexible and can incorporate diverse curve representations including sea level rasters in the first case study and (aggregated) upper forest line vectors in the second. Thanks to tabs’ interoperable and standardised output structure, results can be turned into insightful visualisations to investigate patterns of changing biogeographic systems over time. For example, the first case study illustrates the ease of performing automated labelling procedures for both present-day and submerged islands, and how these outputs can be transformed into maps and graphs to analyse past dynamics of island connections. The second case study illustrates how statistics (i.e., deciles) formalising boundaries of a biogeographic system can be turned into insightful maps and graphs, representing surface area dynamics. In short, we illustrate tabs’ potential to facilitate and enrich scientific analyses with strong visualisations supporting a broad range of biogeographic studies.

Some existing software tools, especially used in the field of island biogeography, are complementary to tabs. The PleistoDist R package (Tan et al. 2023) focuses on visualising and quantifying the effects of eustatic sea-level changes on islands throughout the Pleistocene, providing tools for generating maps and calculating geomorphological metrics. To reconstruct past island areas, PleistoDist relies on a global eustatic and uniform sea-level curve. As a result, it does not account for the geophysical complexity and other types of topographic changes due to geological and geomorphological processes (Tan et al. 2023). In contrast, tabs has a suite of built-in curves, including a spatio-temporal curve that accounts for local sea-level variations (De Groeve et al. 2022a). Furthermore, new and improved spatio-temporal sea-level models, or other curves (e.g., UFL curves), can be easily integrated in tabs. The R package DAISIE (Valente et al. 2015) models island biota change using phylogenetic trees to estimate colonisation, speciation, and extinction rates. In the future, we aim to link or integrate tabs with complementary tools like PleistoDist and DAISIE, enabling the use of tabs outputs to improve predictions of colonisation, speciation, and extinction rates, as well as calculations of geomorphological metrics.

The performance of tabs is dependent on the complexity of the reconstruction. For instance, generating reconstructions using GEBCO and an eustatic curve, which has been the practice in literature until now (Norder et al. 2018; Tan et al. 2023), will be significantly faster than a reconstruction using high-resolution, local topographic models and spatio-temporal curves. Accounting for local variations by defining a correction raster will further decrease the speed of calculations but will result in more accurate reconstructions. Hence, there is a tradeoff between performance and accuracy of the reconstruction and its derived metrics such as area change. The resolution of the input topographic dataset will largely define the quality of the labelling and recognition of shapes. The higher the resolution of the topo dataset, the more likely that labelling is done correctly, and the more likely that shapes are identified as separate units. For instance, islands separated by narrow sea streets, like Sicily (850 m) from mainland Italy and Evia (40 m) from mainland Greece, will not be distinguished as separate shapes, ending up as part of the continent if no high-resolution topographic model is made available. Also, the identification of small shapes, such as islets, will benefit from using accurate topographic data.

Improvements and future developments tabs

Global vs regional reconstructions : Currently, tabs is customised for regional reconstructions across the globe and not for continuous global reconstructions, such as the global shoreline and shelf sea rasters (De Groeve et al. 2022a). Such global reconstructions are also aimed to be integrated in a future release of tabs. However, tabs can be used for global studies by subsequent definition of extents, for instance a list of archipelagos. Using standard iteration functions (e.g., sapply(), lapply()), biogeographic systems can then be reconstructed for each listed area of interest.

Naming : Automated labelling of biogeographic shapes is currently only integrated for shoreline shapes (continents, islands, islets) and mountain ranges. The default labelling option is based on GSV (Sayre et al. 2019) for shoreline shapes or GMBA (Snethlage et al. 2022) for mountain shapes. Hence, labelling based on other geotags, such as coral bed names, requires the definition of custom labelled points or shapes. When using custom points, it is essential to define locations that have a high probability to fall within the present-day biogeographic shapes (such as highest points or centroids), otherwise, points might not intersect with the reconstructed shapes. Another approach could be to export first a present-day reconstruction, to explore which reconstructed shapes are recognised, edit the shape names, and use this input as the labelling dataset.

Furthermore, labelling of reconstructed shapes is primarily based on the location of the highest point observed within the underlying topography. This is a valid naming convention for biogeographic systems which only have a lower range limit, such as shorelines as well as upper forest lines, but not for biogeographic systems defined by a range, such as corals or shelf seas. This is because the highest point is not constant and can shift horizontally with every sea-level change. A more advanced procedure is required to enable labelling of shapes during horizontal biogeographic shifts. Reconstructions that do not require labelling or only require area change reconstructions at the extent level are applicable for any biogeographic range definition.

Corrections : While topographic corrections can be integrated in the model, either through definition of a (list of) rasters or values, it is expected from the user to prepare this input. This is particularly time consuming if corrections are defined as rasters requiring potential data collection, integration and manipulation. Given that there are no integrated databases including rates of topographic change (e.g., uplift, subsidence, erosion, sedimentation), a user will need to first screen the literature to collect known rates for locations within the study area. These rates can be digitised as linear features (in case of linear breakpoints, typically occurring along a tectonic fold and thrust belt) or as point features (in case rates are known for a specific location). Further, it might be necessary to define support points to enhance and facilitate the interpolation to obtain a smoothed raster accounting also for breakpoints. In case rates are not linear and vary over time, the process will need to be repeated for every period. In the future, we aim to facilitate this procedure by creating (1) a point and linear feature database with topographic rates of change and (2) an interpolation function, accounting for breakpoint features.

Isolation metrics and visualisations : We aim to further integrate more functionalities to tabs, including other pre-calculated metrics than area change, such as distance between shapes and timing of isolation, and functions to create publication-ready visualisations including maps, connectivity metrics and curve changes for selected regions or reconstructed shapes.

We encourage researchers interested in using tabs to contact our team for potential collaborations, allowing us to better address user-specific needs and further optimise the tool for diverse applications. For reporting bugs and recommendations, we advise users to create an issue in the tabs gitlab repository (De Groeve et al. 2025d).

Conclusion

In conclusion, tabs is just as intuitive as pressing the tab key—streamlining the complex process of reconstructing spatio-temporal biogeographic configurations. By effortlessly mapping surface area changes, connectivity, and isolation over time, tabs generates structured geospatial datasets ready-to-use for further analysis and visualisation. Its customisable parameters can be adapted to any local environment, making it a valuable tool for researchers aiming to model and quantify biome and ecosystem changes over temporal scales, anywhere on the globe.

Acknowledgements

JDG was supported by institutional resources from the Computational Support Team of the Institute for Biodiversity and Ecosystem Dynamics (IBED), University of Amsterdam (Netherlands). ESR and SGAF acknowledge financial support from Trond Mohn Stiftelse (TMS) and the University of Bergen for the startup grant ‘TMS2022STG03’ to SGAF. We would like to thank scientific illustrator Miranta Kouvari for preparing the graphical abstract supported by the grant ‘TMS2022STG03’ to SGAF and by IBED institutional budgets.

Author contributions

JDG: Writing – Original draft, Writing – Review and Editing, Visualisation, Software, Data curation, Methodology, Project Administration, Conceptualisation. SJN: Writing – Original draft, Writing – Review and Editing, Resources, Methodology, Project Administration, Conceptualisation. KFR: Writing – Review and Editing, Project Administration, Methodology, Resources, Conceptualisation. ESR: Writing – Review and Editing, Project Administration, Validation, Conceptualisation. SGAF: Writing – Review and Editing, Project Administration, Resources, Funding Acquisition, Conceptualisation.

Data accessibility statement

The tabs R package is available on CRAN. The development version (v0.1.1), as used in this publication, is accessible via GitLab [https://gitlab.com/uva_ibed_piac/tabs/-/releases/0.1.1] and archived on Figshare (De Groeve et al. 2025b). All default global datasets used by tabs are available on Figshare (De Groeve et al. 2022b, 2023, 2024). R code for replicating the case studies can be found on Figshare (De Groeve et al. 2025a).

References

  • Ali JR, Aitchison JC (2014) Exploring the combined role of eustasy and oceanic island thermal subsidence in shaping biodiversity on the Galápagos. Journal of Biogeography 41: 1227–1241. https://doi.org/10.1111/jbi.12313
  • Brown SC, Wigley TML, Otto-Bliesner BL, Rahbek C, Fordham DA (2020) Persistent Quaternary climate refugia are hospices for biodiversity in the Anthropocene. Nature Climate Change 10: 244–248. https://doi.org/10.1038/s41558-019-0682-7
  • Camoin GF, Webster JM (2015) Coral reef response to Quaternary sea-level and environmental changes: State of the science. Sedimentology 62: 401–428. https://doi.org/10.1111/sed.12184
  • Cutler KB, Edwards RL, Taylor FW, Cheng H, Adkins J, Gallup CD, Cutler PM, Burr GS, Bloom AL (2003) Rapid sea-level fall and deep-ocean temperature change since the last interglacial period. Earth and Planetary Science Letters 206: 253–271. https://doi.org/10.1016/S0012-821X(02)01107-X
  • De Groeve J, Rijsdijk KF, Rentier ES, Flantua SGA, Norder SJ (2025a) DATA FROM: Temporal Altitudinal Biogeographic Shifts (tabs): R package for reconstructing biogeographic shifts in terrestrial and marine systems over time. figshare. https://doi.org/10.6084/m9.figshare.28891673.v3
  • De Groeve J, Kusumoto B, Koene E, Kissling WD, Seijmonsbergen AC, Hoeksema BW, Yasuhara M, Norder SJ, Cahyarini SY, van der Geer A, Meijer HJM, Kubota Y, Rijsdijk KF (2022a) Global raster dataset on historical coastline positions and shelf sea extents since the Last Glacial Maximum. Global Ecology and Biogeography 31: 2162–2171. https://doi.org/10.1111/geb.13573
  • De Groeve J, Kusumoto B, Koene E, Kissling WD, Seijmonsbergen AC, Hoeksema BW, Yasuhara M, Norder SJ, Cahyarini SY, Geer A van der, Meijer HJM, Kubota Y, Rijsdijk KF (2022b) Spatio-Temporal Relative Sea Level Curve (RSL). https://doi.org/10.21942/uva.20029991.v1
  • Fernández-Palacios JM, Rijsdijk KF, Norder SJ, Otto R, de Nascimento L, Fernández-Lugo S, Tjørve E, Whittaker RJ (2016) Towards a glacial-sensitive model of island biogeography. Global Ecology and Biogeography 25: 817–830. https://doi.org/10.1111/geb.12320
  • Flantua SGA, Hooghiemstra H (2018) Historical connectivity and mountain biodiversity. In: Hoorn C et al. (Eds) Mountains, Climate and Biodiversity. Wiley-Blackwell, 171–186.
  • Flantua SGA, O’Dea A, Onstein RE, Giraldo C, Hooghiemstra H (2019) The flickering connectivity system of the north Andean páramo. Journal of Biogeography 46: 1808–1825. https://doi.org/10.1111/jbi.13607
  • Flantua SGA, Hooghiemstra H, Van Boxel JH, Cabrera M, González-Carranza Z, González-Arango C (2014) Connectivity dynamics since the Last Glacial Maximum in the northern Andes; a pollen-driven framework to assess potential migration. Paleo­botany and Biogeography: A Festschrift for Alan Graham in His 80th Year. Missouri Botanical Garden, St. Louis: 98–123.
  • Flantua SGA, Payne D, Borregaard MK, Beierkuhnlein C, Steinbauer MJ, Dullinger S, Essl F, Irl SDH, Kienle D, Kreft H, Lenzner B, Norder SJ, Rijsdijk KF, Rumpf SB, Weigelt P, Field R (2020) Snapshot isolation and isolation history challenge the analogy between mountains and islands used to understand endemism. Global Eco­logy and Biogeography 29: 1651–1673. https://doi.org/10.1111/geb.13155
  • Ganas A, Parsons T (2009) Three-dimensional model of Hellenic Arc deformation and origin of the Cretan uplift. Journal of Geophysical Research: Solid Earth 114. https://doi.org/10.1029/2008JB005599
  • Intergovernmental Panel on Climate Change (IPCC) [Ed.] (2023) Atlas. In: Climate Change 2021 – The Physical Science Basis: Working Group I Contribution to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge University Press, Cambridge, 1927–2058. https://doi.org/10.1017/9781009157896.021
  • Kealy S, Louys J, O’Connor S (2017) Reconstructing Palaeogeography and Inter-island Visibility in the Wallacean Archipelago During the Likely Period of Sahul Colonization, 65–45 000 Years Ago. Archaeological Prospection 24: 259–272. https://doi.org/10.1002/arp.1570
  • Lambeck K, Rouby H, Purcell A, Sun Y, Sambridge M (2014) Sea level and global ice volumes from the Last Glacial Maximum to the Holocene. Proceedings of the National Academy of Sciences 111: 15296–15303. https://doi.org/10.1073/pnas.1411762111
  • McGuire JL, Lawler JJ, McRae BH, Nuñez TA, Theobald DM (2016) Achieving climate connectivity in a fragmented landscape. Proceedings of the National Academy of Sciences 113: 7195–7200. https://doi.org/10.1073/pnas.1602817113
  • Norder SJ, Baumgartner JB, Borges PAV, Hengl T, Kissling WD, van Loon EE, Rijsdijk KF (2018) A global spatially explicit database of changes in island palaeo-area and archipelago configuration during the late Quaternary. Global Ecology and Biogeography 27: 500–505. https://doi.org/10.1111/geb.12715
  • Norder SJ, Proios K, Whittaker RJ, Alonso MR, Borges PAV, Borregaard MK, Cowie RH, Florens FBV, de Frias Martins AM, Ibáñez M, Kissling WD, de Nascimento L, Otto R, Parent CE, Rigal F, Warren BH, Fernández-Palacios JM, van Loon EE, Triantis KA, Rijsdijk KF (2019) Beyond the Last Glacial Maximum: Island endemism is best explained by long-lasting archipelago configurations. Global Ecology and Biogeography 28: 184–197. https://doi.org/10.1111/geb.12835
  • Papadopoulou A, Knowles LL (2015) Genomic tests of the species-pump hypothesis: Recent island connectivity cycles drive population divergence but not speciation in Caribbean crickets across the Virgin Islands. Evolution 69: 1501–1517. https://doi.org/10.1111/evo.12667
  • Papadopoulou A, Knowles LL (2017) Linking micro- and macroevolutionary perspectives to evaluate the role of Quaternary sea-level oscillations in island diversification. Evolution 71: 2901–2917. https://doi.org/10.1111/evo.13384
  • Rangel TF, Edwards NR, Holden PB, Diniz-Filho JAF, Gosling WD, Coelho MTP, Cassemiro FAS, Rahbek C, Colwell RK (2018) Modeling the ecology and evolution of biodiversity: Biogeographical cradles, museums, and graves. Science 361: eaar5452. https://doi.org/10.1126/science.aar5452
  • Rijsdijk KF, Hengl T, Norder SJ, Otto R, Emerson BC, Ávila SP, López H, van Loon EE, Tjørve E, Fernández-Palacios JM (2014) Quantifying surface-area changes of volcanic islands driven by Pleistocene sea-level cycles: biogeographical implications for the Macaronesian archipelagos. Journal of Biogeography 41: 1242–1254. https://doi.org/10.1111/jbi.12336
  • Rijsdijk KF, Buijs S, Quartau R, Aguilée R, Norder SJ, Ávila SP, de Medeiros SMT, Nunes JCC, Elias RB, Melo CS, Stocchi P, Koene EFM, Seijmonsbergen ACH, de Boer WMT, Borges PAV (2020) Recent geospatial dynamics of Terceira (Azores, Portugal) and the theoretical implications for the biogeography of active volcanic islands. Frontiers of Biogeography 12: e45003. https://doi.org/10.21425/F5FBG45003
  • Salles T, Mallard C, Husson L, Zahirovic S, Sarr A-C, Sepulchre P (2021) Quaternary landscape dynamics boosted species dispersal across Southeast Asia. Communications Earth & Environment 2: 1–12. https://doi.org/10.1038/s43247-021-00311-7
  • Sandel B, Arge L, Dalsgaard B, Davies RG, Gaston KJ, Sutherland WJ, Svenning J-C (2011) The Influence of Late Quaternary Climate-Change Velocity on Species Endemism. Science 334: 660–664. https://doi.org/10.1126/science.1210173
  • Sayre R, Noble S, Hamann S, Smith R, Wright D, Breyer S, Butler K, Van Graafeiland K, Frye C, Karagulle D, Hopkins D, Stephens D, Kelly K, Basher Z, Burton D, Cress J, Atkins K, Van Sistine DP, Friesen B, Allee R, Allen T, Aniello P, Asaad I, Costello MJ, Goodin K, Harris P, Kavanaugh M, Lillis H, Manca E, Muller-Karger F, Nyberg B, Parsons R, Saarinen J, Steiner J, Reed A (2019) A new 30 meter resolution global shoreline vector and associated global islands database for the development of standardized ecological coastal units. Journal of Operational Oceanography 12: S47–S56. https://doi.org/10.1080/1755876X.2018.1529714
  • Simaiakis SM, Rijsdijk KF, Koene EFM, Norder SJ, Van Boxel JH, Stocchi P, Hammoud C, Kougioumoutzis K, Georgopoulou E, Van Loon E, Tjørve KMC, Tjørve E (2017) Geographic changes in the Aegean Sea since the Last Glacial Maximum: Postulating biogeographic effects of sea-level rise on islands. Palaeogeography, Palaeoclimatology, Palaeoecology 471: 108–119. https://doi.org/10.1016/j.palaeo.2017.02.002
  • Snethlage MA, Geschke J, Ranipeta A, Jetz W, Yoccoz NG, Körner C, Spehn EM, Fischer M, Urbach D (2022) A hierarchical inventory of the world’s mountains for global comparative mountain science. Scientific Data 9: 149. https://doi.org/10.1038/s41597-022-01256-y
  • Svenning J-C, Eiserhardt WL, Normand S, Ordonez A, Sandel B (2015) The Influence of Paleoclimate on Present-Day Patterns in Biodiversity and Ecosystems. Annual Review of Ecology, Evolution, and Systematics 46: 551–572. https://doi.org/10.1146/annurev-ecolsys-112414-054314
  • Tan DJX, Gyllenhaal EF, Andersen MJ (2023) PleistoDist: A toolbox for visualising and quantifying the effects of Pleistocene sea-level change on island archipelagos. Methods in Ecology and Evolution 14: 496–504. https://doi.org/10.1111/2041-210X.14024
  • Torres V, Vandenberghe J, Hooghiemstra H (2005) An environmental reconstruction of the sediment infill of the Bogotá basin (Colombia) during the last 3 million years from abiotic and biotic proxies. Palaeogeography, Palaeoclimatology, Palaeoecology 226: 127–148. https://doi.org/10.1016/j.palaeo.2005.05.005
  • Valente LM, Phillimore AB, Etienne RS (2015) Equilibrium and non-equilibrium dynamics simultaneously operate in the Galápagos islands. Ecology Letters 18: 844–852. https://doi.org/10.1111/ele.12461
  • Weigelt P, Jetz W, Kreft H (2013a) Bioclimatic and physical characterization of the world’s islands. Proceedings of the National Academy of Sciences 110: 15307–15312. https://doi.org/10.1073/pnas.1306309110
  • Williams JW, Jackson ST, Kutzbach JE (2007) Projected distributions of novel and disappearing climates by 2100 AD. Proceedings of the National Academy of Sciences 104: 5738–5742. https://doi.org/10.1073/pnas.0606292104
  • Xu W-B, Guo W-Y, Serra-Diaz JM, Schrodt F, Eiserhardt WL, Enquist BJ, Maitner BS, Merow C, Violle C, Anand M, Belluau M, Bruun HH, Byun C, Catford JA, Cerabolini BEL, Chacón-Madrigal E, Ciccarelli D, Cornelissen JHC, Dang-Le AT, de Frutos A, Dias AS, Giroldo AB, Gutiérrez AG, Hattingh W, He T, Hietz P, Hough-Snee N, Jansen S, Kattge J, Komac B, Kraft NJB, Kramer K, Lavorel S, Lusk CH, Martin AR, Ma K-P, Mencuccini M, Michaletz ST, Minden V, Mori AS, Niinemets Ü, Onoda Y, Onstein RE, Peñuelas J, Pillar VD, Pisek J, Pound MJ, Robroek BJM, Schamp B, Slot M, Sun M, Sosinski ÊE, Soudzilovskaia NA, Thiffault N, van Bodegom PM, van der Plas F, Zheng J, Svenning J-C, Ordonez A (2023) Global beta-diversity of angiosperm trees is shaped by Quaternary climate change. Science Advances 9: eadd8553. https://doi.org/10.1126/sciadv.add8553
login to comment