Quaternary landscape dynamics boosted species dispersal across Southeast Asia

Surface processes model and forcing mechanisms

Landscape evolution over the last one million years interval is performed with the open-source modelling code Badlands³⁴. It simulates the evolution of topography induced by sediment erosion, transport, and deposition (Fig. 1a). Amongst the different capabilities available in Badlands, we applied the fluvial incision and hillslope processes, which are described by geomorphic equations and explicitly solved using a finite volume discretisation. In this study, soil properties are assumed to be spatially and temporally uniform over the region, and we do not differentiate between regolith and bedrock. It is worth noting that the role of flexural responses induced by erosion and deposition is also not accounted for. Under these assumptions, the continuity of mass is governed by vertical land motion (U, uplift or subsidence in m/yr), long-term diffusive processes and detachment-limited fluvial runoff-based stream power law:

where z is the surface elevation (m), t is the time (yr), κ is the diffusion coefficient for soil creep³⁴ with different values for terrestrial and marine environments, ϵ is a dimensional coefficient of erodibility of the channel bed, m and n are dimensionless empirically constants, that are set to 0.5 and 1, respectively, and PA is a proxy for water discharge that numerically integrates the total area (A) and precipitation (P) from upstream regions³⁴.

Both κ and ϵ depend on lithology, precipitation, and channel hydraulics and are scale dependent³⁴. All our landscape evolution simulations are running over a triangular irregular network of ~18. e⁶ km² with a resolution of ~5 km, and outputs are saved every 1000 yr.

The detachment-limited fluvial runoff-based stream power law is computed with a ({{{{{{{mathcal{O}}}}}}}}(n))-efficient ordering approach⁵⁴ based on a single-flow-direction approximation where water is routed down the path of the steepest descent. The flow routing algorithm and associated sediment transport from source to sink depend on surface morphology, and sediment deposition occurs under three circumstances: (1) existence of depressions or endorheic basins, (2) if local slope is less than the aggregational slope in land areas and (3) when sediments enter the marine realm³⁴. Submerged sediments are then transported by diffusion processes defined with a constant marine diffusion coefficient³⁴.

All landscape simulations are constrained with different forcing mechanisms, and five scenarios were tested (Supplementary Table 2).

First, we impose precipitation estimates from the PaleoClim database^38,39,40. These estimates are products from paleoclimate simulations (coupled atmosphere-ocean general circulation model) downscaled at approximately the same resolution as our landscape model (~5 km at the equator). Annual averages of precipitation rates are then used to provide rainfall trends in our simulations based on the ten specific snapshots available (from the mid-Pliocene warm period to late Holocene and present day). Between two consecutive snapshots, we assume that precipitation remains constant for the considered time interval. For exposed regions that are considered flooded in the PaleoClim database, we define offshore precipitations using a nearest neighbour algorithm where closest precipitation estimates are averaged from PaleoClim inland regions. To evaluate the role of precipitation variability on landscape dynamics, we also run a uniform rainfall scenario (2 m/yr obtained by averaging the annual precipitation rates from the PaleoClim database).

Secondly, the models are forced with sea level fluctuations known to play a major role in the flooding history of the Sunda Shelf^11,13,53. Two sea level curves are tested (Supplementary Fig. 1d). To account for the inherent uncertainty in reconstructed sea level variations, we chose a first curve³⁷ obtained from a sea level stack constructed from five to seven individual reconstructions that agrees with isostatically adjusted coral-based sea level estimates at both 125 and 400 ka. The second one is taken from the global sea level curve reconstruction³⁶ based on benthic oxygen isotope data and has been recently used to reconstruct the subsidence history of Sundaland^17,18.

The last forcing considered in our study is the tectonic regime. First, we chose to explore a non-tectonic model based on the default assumption of stability for the shelf¹⁷. Secondly, we assumed a uniform subsidence rate of −0.25 mm/yr recently derived from a combination of geomorphological observations, coral reef growth numerical simulations and shallow seismic stratigraphy interpretations¹⁷. Then, to represent the regional variations in the tectonic regime, we have compiled and digitised a number of calibration points (Supplementary Fig. 1b and Table 1) that were used to produce a subsidence and uplift map by geo-referencing calibration points and available tectonic polygons, and by Gaussian-smoothing and normalising the uplift and subsidence rates between the calibrated range to avoid sharp transitions in regions without observations. The resulting map does not account for fine spatial scale tectonic features such as fault systems^43,55 or orogenic and sedimentary related isostatic responses. It rather represents a regional vertical tectonic trend with an overall uplift of Wallacea and NW Borneo regions and long-wavelength subsidence of Sunda Shelf and Singapore Strait¹⁷.

Landscape evolution model calibration

The landscape models start during the Calabrian in the Pleistocene Epoch, one million years before the present. At each time interval, the landscape evolves following Eq. (1) and the surface adjusts under the action of rivers and soil creep (Fig. 1a). In addition to surface changes, we extract morphometric characteristics such as drainage basins extents, river profiles lengths (Fig. 3 and Supplementary Fig. 2), distance between main rivers outlet (Supplementary Fig. 3) and tracks the cumulative erosion and deposition over time (Fig. 1b and Supplementary Fig. 1d).

For model calibration, we perform a series of steps consisting in adjusting the initial elevation and the erosion–deposition parameters (i.e., κ and ϵ in Eq. (1)) to match with different observations.

The initial paleo-surface is obtained by applying the uplift and subsidence rates backwards to calculate the total change in topography for the 1 Myr interval. Then, we test the simulated paleo-river drainages against results from a combination of phylogenic studies^9,13 and paleo-river channels and valleys found from seismic and well surveys^41,42,44. Iteratively, we modify our paleo-elevation to ensure those main river basins (e.g., Johore, Siam, Mekong, East Sunda) encapsulate the paleo-drainage maps reconstructed using lowland freshwater taxa described in¹³ (Supplementary Fig. 1a and Table 4) and that the major rivers follow paleo-rivers systems derived from both 2D and 3D seismic interpretations (Fig. 1b).

For surface processes parametrisation, we tested different ranges of diffusion and erodibility coefficients and compared the final sediment accumulation across the Sunda Shelf (Fig. 1b) using estimated deposit thicknesses^41,42,43,44. The Sunda Shelf is predominantly experiencing deposition over the past 500 kyr and increases in deposition are positively correlated with periods of sea level rise (i.e., Pearson’s coefficients for correlation with sea level above 80%, Supplementary Fig. 1d). After exploring a range of values, we set κ values to 1. e⁻² and 8. e⁻² m²/yr for terrestrial and marine environments and ϵ between 2.5 and 7.5 e⁻⁸yr⁻¹ for the different scenarios to fit with chosen surveys dataset (Supplementary Table 2 and 3).

Upon uniform subsidence case (−0.25mm/yr), flooding is limited, and the shelf only undergoes two full marine transgressions (>60% of the shelf flooded) around 125 ka and during the last 10–20 kyr (Supplementary Fig. 1c). Upon spatially variable tectonics (non-uniform subsidence), partial flooding events are more pervasive, with higher magnitudes and greater temporal durations. Due to the shallow and flat physiography of Sundaland, we also note that even small increases in sea level amplitudes (<10 m, as bore between our two sea level curves^36,37) could trigger up to a 30% increase in shelf area inundation (for instance the red and blue curves in Supplementary Fig. 1c during MIS7, ~200 ka). In addition, we further tested our combination of forcing mechanisms, initial paleo-surface and model parameters by counting the number of flooding events simulated within the Malay Basin on the Sunda Shelf (box A in Fig. 1b). For both sea level curves and variable tectonic map, we recorded at least five marine incursion episodes, similar to the number of events found for the same area in⁴¹ based on interpreted facies characteristics site survey borehole data and 3D seismic sections analysis.

Landscape elevational connectivity

The landscape elevational connectivity (LEC) is a measure of the energy required by a pool of adapted species to move outside its usual niche width, up and down their initial elevation range, to spread and colonise any other habitats³³. This metric can predict local species richness (α diversity) obtained from full metacommunity models when applied to real landscapes³³. In a recent study⁵², comparisons between generated Badlands model landscapes using both uniform and orographic precipitation conditions have shown similar results when the metric is applied over geological time scales.

Landscape elevational connectivity at node i (LEC_i) relies on the following set of equations³³:

where C_ji quantifies the closeness between sites j and i with respect to elevational connectivity and measures the cost for a given species adapted to cell j to move to cell i. This cost is a function of elevation and evaluates how often species adapted to the elevation of cell j have to travel outside their optimal species niche width (σ) to reach cell i⁵². p = [k₁, k₂, . . . , k_L] (with k₁ = j and k_L = i) are the cells in the path p from j to i.

Solving Eq. (2) relies on Dijkstra’s algorithm accounting for diagonal connectivity between cells and implemented using the scikit-image library⁵⁶. Even with a balanced parallel implementation, the LEC_i calculation is slow as a Dijkstra tree for all nodes must be created, and least-cost distances between each node and all others have to be calculated.

In this study, to decrease computation cost, we have modified the initial algorithm⁵² to limit the number of nodes that needs to be visited when calculating least-cost distances. We gave each cell an initial amount of energy that is consumed as the pool of species adapted to a specific cell is moving across the topography. The energy expenditure depends on (1) the differences in elevation and (2) the distances between neighbouring cells. Here, the initial amount of energy is set to 2000 and we weight the energy dissipation based on travelled distances by 0.4% assuming that species have the ability to move easily over long distances if they stay in their optimal niche range. This assumption is valid in our case as we are simulating the connectivity associated with a large area (~18. e⁶ km²) and over geological time scale.

Connectivity mapping from current density modelling

By considering all possible pathways, circuit theory offers a better modelling alternative to least-cost path approaches as it captures all the dynamics at play in travel decisions based on provided resistance maps. We chose Circuitscape³⁵ to model multiple pathways. The software uses random walk and electric current running through a circuit. Electric current runs across our cost surfaces between predefined source points. We position these points across Sundaland (approximately 250 km apart) chosen along the outer margin (≃100 m above sea-level) of the maximum fully submerged shelf coastline. Circuitscape functions as a graph, where each cell centre is a node connected to neighbouring nodes with links³⁵. The graph is interpreted as a circuit, and links have cost values combining the laws of electricity to estimate species flow. Effective resistance, voltage and current is then calculated over the landscape between the prescribed points with Ohm’s and Kirchhoff’s laws³⁵.

In Ohm’s law, voltage V applied to a resistor R gives the current I through I = V/R; as a result, a lower resistance (e.g., low geomorphic cost) in the landscape will correspond to higher flow of species. Kirchhoff’s law deals with effective resistance; when nodes are connected to several resistors the effective resistance will be the sum of the resistances: R = R₁ + R₂.

Connectivity maps statistical analyses

To evaluate the relative contribution of each of the three geomorphological features, we compute current density maps for different levels of Sunda Shelf exposure using different combinations of resistance surfaces. We then tested each feature independently (e.g., only landscape elevational connectivity, only distance to rivers (with rivers as barriers and as corridors) and only local slope), then pairwise costs (e.g., landscape elevational connectivity and distance to rivers, landscape elevational connectivity and local slope, local slope and distance to rivers) and finally all costs combined (Supplementary Fig. 5a).

Qualitative comparisons of current density maps rely on visual interpretation and are often altered by the choice of colour scale used to distinguish regions of high connectivity⁴⁷. To perform a better evaluation, we first express current flow values (c) as z-scores (({z}_{sc}=(c-overline{c})/sd)) by subtracting the mean current value (overline{c}) and dividing it by the standard deviation sd. We then used three different thresholds to estimate regions that have positive mean values (i.e., z_sc > 0) and values higher than one and two standard deviations (z_sc > 1 and 2, respectively, Supplementary Fig. 1b). The approach provides a quantitative assessment of flow maps sensitivity to the chosen resistance maps.

To gain additional insights into the distribution of connectivity regions across the shelf, we also employed a local spatial autocorrelation indicator, namely the Getis-Ord Gi⋆ index⁵⁷. This hotspot analysis method assesses spatial clustering of the obtained current density maps, and the resultant z-scores provide spatially and statistically significant high or low clustered areas. The approach consists in looking at each local current value relative to its neighbouring one. From this spatial analysis, we extract both statistically significant hot and cold spots for each combination of resistance surfaces (Supplementary Fig. 5c). To extract statistically significant and persistent biogeographic connectivity areas across the exposed Sunda Shelf, we then combine all hotspots together and define preferential migration pathways as regions having a positive Gi⋆ z-scores for all resistance surfaces combination.

We used the function zscore in the SciPy stats package to obtain the z-scores and the ESDA library for the Gi⋆ indicator computation.

Modelling assumptions and limitations

There are a number of important caveats for interpreting our modelling results.

First, we made several assumptions related to our transient landscape evolution simulations. A single-flow direction algorithm⁵⁴ was used to simulate temporal changes in river pathways. Recent work⁵⁸ has shown that this algorithm might lead to numerical diffusion, fast degradation of knickpoints and underestimation of river captures particularly in flat regions. One way to address this would be to use a multiple flow direction method⁵⁹ which allow for a better representation of flow distribution across the landscape. In this study, we also assumed a uniform and invariant soil erodibility coefficient for the entire domain and a detachment-limited erosion law. Even though the erodibility coefficient was calibrated independently for each simulation (Supplementary Table 3), soil cover and properties vary notably between Borneo, Sumatra, Java and the Malay Peninsula and soil conditions for the exposed sea floor would have changed significantly over successive flooding events¹². Badlands software³⁴ allows for multiple erodibility coefficients representing different soil compositions to be defined, and this functionality could be used to evaluate the impact of differential erosion on physiographic changes. Similarly, several transport-limited laws are also available and could be compared against our detachment-limited simulations.

A second set of simplifications lies in the climatic conditions (i.e., rainfall regimes) used to force our simulations. We relied on the PaleoClim database⁴⁰ which contains nine high-resolution paleoclimate dataset^38,39,40 corresponding to specific time periods (4.2–0.3 ka, 8.326–4.2 ka, 11.7–8.326 ka, 12.9–11.7 ka, 14.7–12.9 ka, 17.0–14.7 ka, ca. 130 ka, ca. 787 ka and 3.205 Ma). The climate simulations from which these time periods are extracted do not consider emerged Sunda Shelf for the oldest inter-glacial events which can result in incorrect climatic pattern⁶⁰. From 0.3–17 ka, precipitation fields in PaleoClim are obtained from the TRaCE21ka transient simulations of the last 21 kyr run with the CCSM3 model⁴⁰. Although Fordham et al.³⁹ show that precipitation errors range from 10–200% in their modern experiment, the paleoclim dataset provides a statistical downscaling method that includes a bias correction (namely the Change-Factor method, in which the anomaly between the modern simulation and observations is removed from the paleoclimate experiment) allowing the use of the model for paleoclimate studies⁴⁰. The very same technique is applied for 130 ka and 787 ka fields that were obtained with different GCMs (namely HadCM3 and CCSM2). Given the absence of a million-year long transient climate simulation, we oversimplified the climatic conditions by considering that precipitation distribution and intensity remain constant between two consecutive intervals, generating an artificial stepwise evolution of rainfall through time. To evaluate the sensitivity of physiographic responses on the Sunda Shelf to precipitation, we ran a model with uniform rainfall over 1 Myr (scenario 4). Despite changes in the timing and extent of basins reorganisation (Supplementary Fig. 2 and Fig. 3b), we found limited differences in terms of flooding history and erosion/deposition patterns when compared with scenario 5 (Supplementary Fig. 1c, d and Supplementary Table 2). Recent work⁶⁰ suggests clear regional responses induced by the emerged Sunda Shelf with seasonal enhancement of moisture convergence and continental precipitation induced by thermal properties of the land surface. This could significantly impact our simulation results. However, and at the time of writing, more continuous high-resolution paleoclimatic simulations considering the shelf as an emerged continental platform were still unavailable. Using high-resolution multi-model outputs would allow to target the uncertainty on climatic maps⁴ and will surely represent a significant advance for future studies. One approach would have used the orographic rainfall capability⁶¹ available in Badlands. The method is better suited to run generic simulations but falls short when applied to real cases as it relies on imposing paleo-environmental boundary conditions (e.g., temporal changes in wind direction and speed, moisture stability frequency or depth of moist layer) difficult to obtain for Earth-like model applied over geological time scales.

Finally, our species-agnostic approach assumes an equally weighted cost between the three considered geomorphic features and does not account for additional factors (temperature, vegetation cover, solar radiation to cite a few), which are all important when assessing landscape connectivity for wildlife. Most importantly, we model connectivity at very large scales (5 km resolution). Often, species are highly influenced by microclimates and small-scale topography⁴⁷. From our regional-scale simulations and hotspot analysis (Fig. 6), higher resolution models focusing on highly connected regions (across the Gulf of Thailand and Siam basin) could be applied to produce more detailed representations of species migration in the region. In addition, current flow field calculations from Circuitscape³⁵ rely on randomly selecting nodes around the region of interest. For connectivity analysis, we used 33 terrestrial points located around the perimeter of the buffered Sundaland area (white contour line in Fig. 1b). Using a selection of nodes in a buffered region allows to reduce the bias in current density estimates⁴⁶. However, bias might depend on the buffer size compared to the study area as well as the number of nodes selected^46,47. Because of memory limitations and the great number of computed grids used to cover the past 500 kyr, we made a trade-off between buffer size and the number of selected points for pairwise calculations. Additional experiments could possibly be tested to evaluate bias in the proposed connectivity maps potentially using a tilling approach to reduce cell number⁴⁵.

Source: Ecology - nature.com