Mapping landscape connectivity as a driver of species richness under tectonic and climatic forcing
Species distribution and richness ultimately result from complex interactions between biological, physical, and environmental factors. It has been recently shown for a static natural landscape that the elevational connectivity, which measures the proximity of a site to others with similar habitats, is a key physical driver of local species richness. Here we examine changes in elevational connectivity during mountain building using a landscape evolution model. We find that under uniform tectonic and variable climatic forcing, connectivity peaks at mid-elevations when the landscape reaches its geomorphic steady state and that the orographic effect on geomorphic evolution tends to favour lower connectivity on leeward-facing catchments. Statistical comparisons between connectivity distribution and results from a metacommunity model confirm that to the 1st order, landscape elevation connectivity explains species richness in simulated mountainous regions. Our results also predict that low-connectivity areas which favour isolation, a driver for in situ speciation, are distributed across the entire elevational range for simulated orogenic cycles. Adjustments of catchment morphology after the cessation of tectonic activity should reduce speciation by decreasing the number of isolated regions.
The idea that mountainous landscapes play a role in biological evolution has a long history that can be tracked back to Darwin and Wallace, when fauna boundaries were noted to correspond to physiographic discontinuities and gradients (Wallace, 1860). Over geological timescales (millions of years), surface processes including erosion and incision conspire to undo surface uplift driven by tectonic and geodynamic processes. These competing processes can convert landscapes of low elevation and relief, a homogeneous environment, and low resistance to migration into complex landscapes with sharp environmental gradients and fragmented habitats separated by migratory corridors (Steinbauer et al., 2016).
Several studies have shown that on geological timescales, changes in landscape morphology stimulate migratory behaviour, dispersal, and redistribution as species track their optimum habitats (Gillespie and Roderick, 2014; Craw et al., 2015; Muneepeerakul et al., 2019). Species running out of favourable habitats are forced to coexist with other ones and adapt, which leads to speciation, increasing endemism and biodiversity (Smith et al., 2014). As the mountainous landscape becomes more complex and diverse, environmental gradients increase and species have to move across shorter distances to track suitable habitats (Smith et al., 2014; Elsen and Tingley, 2015) and find refuges to survive both long-term mountain-building processes and short-term climatic changes. Hence, mountains host a disproportionately large fraction of terrestrial species (Badgley, 2010; Hoorn et al., 2010; Steinbauer et al., 2018), illustrating the tectonic influence on ecology and evolutionary biology.
Determining the underlying geophysical drivers of species richness requires accounting for factors (gradients of topography, habitat capacity, humidity, temperature, altitude, and solar exposure) that often covary with elevation (Lomolino, 2001; Kearney and Porter, 2009; Martensen et al., 2017; Liu et al., 2018). Recently, Bertuzzo et al. (2016) proposed a metric to assess biodiversity patterns. The metric, called landscape elevational connectivity (LEC), accounts for two fundamental landscape geomorphological properties: (1) within mountains ranges the maximum surface area peaks for mid-elevations rather than the lowest elevations, and (2) at mid-elevations species have the option to disperse up or down to well-connected suitable patches. As a consequence of these two properties, species richness (e.g. diversity) is predicted to be highest at mid-elevations, showing a hump-shaped pattern (Bertuzzo et al., 2016; Lomolino, 2001; Rahbek, 1995, 2005; McCain and Grytnes, 2010; Kessler et al., 2011) because diversity is promoted by increased area and connectivity (MacArthur and Wilson, 1967; Solé and Bascompte, 2006; Economo and Keitt, 2010; Carrara et al., 2012).
The results from Bertuzzo et al. (2016) were obtained on a static landscape. The aim of this study is to extend the analysis over the entire orogenic cycle from the mountain-building to relaxation phases using a synthetic example and to map over space and time the landscape connectivity. First we analyse the temporal and spatial LEC distribution based on the results of a landscape elevation model (Salles, 2016) under uniform tectonic and variable climatic conditions. Then, we run a zero-sum metacommunity model (Hubbell, 2001) on the simulated synthetic surfaces and, as in Bertuzzo et al. (2016), we find that the LEC quantity explains to the 1st order species richness in mountainous regions and can be used to infer changing patterns of biodiversity resulting from the effect of climate, tectonic, and geomorphic processes. The obtained results provide new insights on how the tempo imposed by landscape morphological changes over geological timescales affects the distribution of species richness in mountainous regions.
2.1 Coupled landscape evolution and orographic precipitation model
In this paper, the numerical simulations of mountain evolution are performed with the Badlands landscape evolution model (Salles, 2016; Salles et al., 2018a), and fluvial erosion rates and predicted sediment transport in rivers are solved using the stream-power law (SPL). The SPL relates the erosion rate I to the product of mean annual net precipitation rate (), drainage area (A), and local slope (S) and takes the form
where the erodibility coefficient κe is controlled by climate and lithology, and m and n are positive exponents (Chen et al., 2014) that mostly depend on the nature of the dominant erosional mechanism (Dietrich et al., 1995). Regardless of its simplicity, the SPL (Eq. 1) allows us to simulate the main geomorphic characteristics of mountainous regions, in which landscape evolution is dominated by a detachment-limited erosion regime (Whipple and Tucker, 1999).
where z is the local landscape elevation and κd is the diffusion coefficient. This formulation assumes a linear dependency of superficial sediment transport on the topographic gradient.
Accounting for the two processes defined above, the equation of mass conservation driving landscape temporal evolution is expressed as
where t is time and U is the rock uplift rate. It is worth noting that the landscape evolution scenarios presented in this study neglect spatial and temporal variations in tectonic evolution and rock strength, as well as the influence of sediment characteristics on river incision.
Precipitation patterns are known to impact geomorphic evolution at local (ridge–valley) to regional (full mountain range) scales (Anders et al., 2008; Bonnet, 2009). To estimate the controls of precipitation variability on topography evolution and on associated landscape connectivity, the set of landscape evolution equations presented above is coupled with a linear orographic precipitation model (Smith and Barstad, 2004). The approach already implemented in Badlands (Salles et al., 2018b) only accounts for the 1st-order physics of orographic precipitation, computes rainfall under idealised climatic conditions defined by a specific wind velocity and direction (Anders et al., 2008), and assumes steady, uniform, and saturated airflow (Fig. 1a).
In the study, analyses of catchment dynamics are performed using river longitudinal profiles as well as χ plots and maps. χ analysis is a method of extracting information from channel profiles that attempts to compare channels with different discharges (Perron and Royden, 2013). The longitudinal coordinate χ has dimensions of length and is linearly related to the elevation (Mudd et al., 2014). The quantity χ at any channel point depends on the river network geometry and estimates the dynamic state of river basins (Willett et al., 2014):
with x the distance along a given channel from the river base (the outlet) xb and the point considered, A0 a scaling area, m and n the coefficients of the SPL (Eq. 1), and A(x′) the upstream drainage area at position x′. Mapping χ along channel networks and comparing values of χ across drainage divides provides a measure of the equilibrium state of a river network under the influence of tectonic uplift and river erosion.
2.2 Landscape elevational connectivity
Landscape connectivity encapsulates the combined effects of (1) landscape morphological structure and (2) the species ability to move outside its usual niche width (Tischendorf and Fahrig, 2000). As such it is both species- and landscape-specific. Here, we generalise the concept and assume that any position in the simulated mountainous landscape corresponds to a particular habitat containing a pool of adapted species able to move in a limited elevation range up and down their initial and preferred habitat elevation.
We chose to map landscape connectivity and its distribution by measuring the landscape elevational connectivity (LEC) proposed by Bertuzzo et al. (2016). This metric does not account for a specific species but rather focuses on a pool of species assuming a niche width that can be a percentage of the mountain elevation range or fixed. Here, we use the bioLEC Python package to measure the LEC index (Salles and Rey, 2019). The LEC is calculated on the evolving landscape obtained from Badlands at discrete time intervals and measures how easily species living in other habitats can spread and colonise a given point.
Considering a 2-D lattice made of N squared cells, the LEC for cell i (LECi) is given by
where Cji quantifies the closeness between sites j and i with respect to elevational connectivity. Cji measures the cost for a given species adapted to cell j to spread and colonise 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 (Fig. 1b), assuming that species niche is represented by a Gaussian function of the elevation.
The estimation of Cji requires the computation of all the possible paths p from j to i and is defined as the maximum closeness value along these paths. Following Bertuzzo et al. (2016), Cji is expressed as
where (with k1=j and kL=i) represents the cells in the path p from j to i.
Equation (6) is solved for each cell j using Dijkstra's algorithm (Dijkstra, 1959) with diagonal connectivity between cells. The approach in bioLEC is based on the scikit-image Dijkstra algorithm (van der Walt et al., 2014). For each cell j, the algorithm builds a Dijkstra tree that branches the given cell with all the cells defining the simulated region. Edge weights are set equal to the square of the difference between the considered vertex elevation () and zj. The least-cost distance between j and i is then calculated as the minimum sum of edge weights obtained from the cells along the shortest path (Fig. 1b).
Calculation of LEC values over the entire simulated region (∼1 M points in this study) is slow. Here we use the parallel strategy available in the bioLEC package wherein Dijkstra trees for all paths are balanced and distributed over multiple processors using a message-passing interface (MPI). Using this approach, LEC computation on the synthetic landscapes obtained in this study is less than 4 min when distributed over 240 processors.
2.3 Experimental setting
Our experimental design consists of a 100 × 100 km surface at a resolution of 100 m (Fig. 1a). The initial topography is a flat area over which random noise (≤1 m) is applied. To estimate landscape connectivity through a complete orogenic cycle, simulations are run for 10 Myr and a uniform uplift of 1 mm yr−1 is applied during the first 5 Myr (e.g. mountain-building phase) on all surface nodes with the exception of boundary points that remain fixed over the runs (Fig. 1a). A second phase, corresponding to the cessation of tectonic activity over the remaining 5 Myr, simulates the lowering of interfluves and the erosional decay of the mountain system (e.g. mountain relaxation phase).
In this study and for simplicity, the SPL parameters (Eq. 1) are set constant and do not change between simulations. We assign a uniform erodibility coefficient κe of yr−1, and exponents m and n are set to 0.5 and 1, respectively; these exponents values are commonly used for eroding fluvial systems (Whipple and Tucker, 1999; Ferrier et al., 2013). For the hillslope processes (Eq. 2), a constant diffusion coefficient of 0.1 m yr−1 is chosen.
To estimate the influence of climate on geomorphic changes and in turn on landscape connectivity, we perform two simulations (Fig. 1a) with (1) a uniform precipitation rate of 1 m yr−1 over the simulated 10 Myr and (2) an orographic precipitation model considering a prevailing wind direction from the south with a wind speed of 0.5 m s−1 and other orographic parameters set within the range proposed by Smith and Barstad (2004). In both cases, simulated mountain ranges are characteristic of dendritic, erosional landscapes.
The computation of LEC over time requires the definition of species niche width (σ). In the study conducted here, we limit our analysis to the case in which all species have the same niche width at a specific time assuming a ratio constant. As the landscape changes through time, the elevation range is modified and so is species niche width. Therefore, we implicitly assume that the characteristic response of the landscape to tectonic and climatic forces happens on a temporal scale much longer than the effective evolution and adaptation of individual species (Mokany et al., 2012; Steinbauer et al., 2016).
3.1 Connectivity distribution under uniform conditions
In this first set of results, we consider a geographic domain free of environmental gradients under a constant precipitation rate. As the region is uniformly uplifted at a rate of 1 mm yr−1 during the first 5 Myr, landscape dissection by riverine and hillslope processes leads to valley deepening and the development of self-organised and stable large drainage basins (Fig. 2a).
Following the initial orogenic phase of plateau uplift and landscape incision, river longitudinal profiles exhibit a concave upward shape, a common signature of bedrock river systems (red lines in Fig. 3b). This typical shape is a direct outcome of the detachment-limited stream-power law model and is the result of the integrated effect of tectonic, climatic, and bedrock erodibility. After 4 Myr of evolution, the landscape has fully adjusted to the imposed constant climatic and tectonic forcing, as shown by the increasing linearity of river χ plots (blue lines in Fig. 3b) (Perron and Royden, 2013). After the cessation of uplift at 5 Myr, erosion induces the rapid lowering of elevations (Fig. 2a) and slope gradients (Fig. 3d) through valley widening that ultimately leads to a low-relief peneplain.
From an elevational niche perspective, all connectivity (e.g. migratory) paths on a flat landscape are equal. This is not the case on a complex landscape where connectivity can only occur along a network of corridors providing species with their elevational requirements (Smith et al., 2014; Elsen and Tingley, 2015). The mapping of the spatial distribution of LEC (Fig. 2a) and associated gradients of LEC (Fig. 3a, c) reveals a migration of maximum LEC regions – during the initial orogenic phase – from the highest elevations (e.g. uplifted plateau) to mid-elevations as the landscape stabilises.
This redistribution leads to two regions of maximum LEC, with the highest one attributed to the symmetrical nature of the simulated synthetic landscape. During the relaxation phase, the LEC peak migrates to lower elevations as the floor of the valleys widens (Fig. 2b). After dissection of the uplifted plateau, the model predicts that higher LECs are found at intermediate to high elevations (Fig. 3a, thick blue line). This result agrees with the hump-shaped elevational gradients of mean species richness widely observed in nature (Grytnes and Vetaas, 2002; Colwell et al., 2004; Brehm et al., 2007). We also see that for areas at similar elevation, connectivity values could differ significantly as shown by standard deviation curves (highlighted blue regions at 5.0 and 6.0 Myr in Fig. 3a, c). It suggests that using only elevation as a predictor of landscape connectivity is unlikely to be a good proxy.
We also observe that, irrespective of the scale of observation, LEC distribution follows a similar pattern for specific time intervals even if the considered domains traverse different elevational extents. For example, the result in Fig. 2c shows similar mid-elevation maximum LEC (i.e. peak of diversity) over valley, catchment, and regional scales.
Our results show that under uniform tectonic and climatic forcing, landscape dynamics exert a 1st-order control on the landscape connectivity distribution, which is mostly distributed at mid-elevation and peaks on the highest flanks when geomorphic steady state is reached (Fig. 3a, thick blue line). The peak position is related to the symmetrical nature of the chosen boundary conditions and species niche width (σ=0.1) as shown in Sect. 4.1.
3.2 Orographic effect on patterns of landscape connectivity
To explore the influence of heterogenous precipitation on connectivity, we design a new set of experiments using a coupled climate and landscape evolution model in which precipitation responds through space and time to the build-up of topography (i.e. orographic precipitation; Fig. 4a, b).
In this experiment, a prevailing northward-directed wind condition creates a specific distribution of rainfall patterns across the mountain range, with precipitation maxima occurring predominantly on windward-facing areas. Induced changes in hydrological regime drive drainage-divide migration (Fig. 5a) and the development of asymmetric topography (Roe et al., 2003; Anders et al., 2008). Topography on the leeward side (rain shadow area) is formed by tectonic uplift with limited erosion, whereas topography on the windward side is predominantly erosional. Over geological timescales, the differential distribution of rainfall influences the overall topographic evolution of the simulated mountain range. First, it affects river profiles by changing concavity (Roe et al., 2003) and reducing erosional rates on the drier regions. As a consequence, a decrease in discharge produces steeper slopes to compensate for the imposed uniform uplift rate (Fig. 4c). Orographic rain also implies continuous disequilibrium conditions (Willett et al., 2014; Whipple et al., 2016), affecting river networks through reorganisation and concomitantly changing the topologies and geometries of drainage basins (Fig. 4e).
Similar to the uniform precipitation model, we observe a temporal migration of the LEC maximum from the uplifted plateau towards the bottom of the valleys. However, we find some specific features driven by the geomorphic responses to orographic precipitation (as shown by the LEC distribution at steady state between Figs. 2a and 4d) with a clear correlation between sites with high- or low-connectivity areas and windward- or leeward-facing catchments. It shows how uplift and its effect on regional climate influence landscape connectivity with the development of geographic connections or barriers as topography changes.
The study of two specific catchments sharing a common drainage divide further highlights these effects. In Fig. 5, we note that both leeward and windward catchments exhibit similar hump-shaped patterns of LEC, indicative of mid-elevation peaks, with higher connectivity on the windward-facing area compared to the leeward-facing one. Temporal evolution analysis through the complete orogenic cycle shows optimum landscape connectivity during the mountain-building phase at middle to high elevations as illustrated by the broader ranges of LEC values at 2.5 Myr compared to the ones at 5.5 Myr (Fig. 5b).
The results described above quantify the temporal and spatial evolution of landscape connectivity patterns during a complete orogenic cycle considering both uniform and orographic precipitation scenarios. In this section, we build upon the work from Bertuzzo et al. (2015, 2016). First, we test the influence of different boundary and species niche width conditions on landscape connectivity. Then, we explore the relationship between landscape elevational connectivity and biodiversity using a similar metacommunity model as the one proposed in Bertuzzo et al. (2016) and discuss the possible roles of geomorphological changes in patterns of species richness. By defining isolated regions as places with low connectivity, we then analyse how they evolve during simulated orogenic cycles.
4.1 Impact of forcing conditions and species niche width on connectivity mapping
In our experimental setting, we have deliberately chosen a simple model of an evolving landscape similar to existing analogue and numerical models (Bonnet and Crave, 2003; Bonnet, 2009; Willett et al., 2014; Salles and Hardiman, 2016). Under such conditions, the resulting landscape is symmetric with the formation of a mountain range composed of four low-angled ridges ending at the four corners of the model. In this section, we slightly modify these initial tectonic boundary conditions and perform a new simulation with outlets along only one of the domain boundaries (Castelltort and Yamato, 2013). The other landscape parameters (erodibility value κe, SPL exponents m and n, and hillslope diffusion κd defined in Sect. 2.3) remain unchanged and the simulated landscape is again purely erosional (similar to the homogeneous model proposed in Roy et al., 2015). Evaluating how LEC will change for constructional landscapes (Craw et al., 2015, 2019) requires different tectonic conditions and is beyond the scope of this paper.
The general pattern of geomorphological evolution under uniform precipitation is presented in Fig. 6a for the first 5 Myr and corresponds to the mountain-building phase. As for the previous cases, the evolution of the experiment involves a growth phase followed by a steady-state phase (Bonnet and Crave, 2003; Castelltort and Yamato, 2013). During the initial plateau uplift (e.g. growth phase), some topographic incisions form along the western open border of the model. As uplift continues, these incisions grow and propagate inward until there is complete dissection of the plateau (third map on the right-hand side of Fig. 6a). The obtained fluvial landscapes display dendritic networks, and river longitudinal profiles present a concave upward shape characteristic of bedrock river systems (as for Fig. 3b).
To evaluate the impact of these new boundary conditions on connectivity calculation, we compute the LEC at steady state (Fig. 6b) over the entire domain and for three quadrants (Fig. 7). In addition, we test two values of 0.1 and 0.2 for the species niche width ratio () defined in Eq. (6). As shown in the results (Sect. 3) and for the Swiss Alps in Bertuzzo et al. (2016), the frequency distributions of the elevation are hump-shaped regardless of the spatial extent of the considered domain (red lines representing the entire domain in Fig. 6 or three smaller regions as in Fig. 7). From the results presented in Fig. 7, we find that connectivity distribution varies significantly with respect to variations of niche width σ. As the niche becomes wider (higher σ), the connectivity maps (Figs. 6c and 7c) reveal a clear spatial pattern, with valleys and mountain tops characterised by lower LEC (Bertuzzo et al., 2016). As σ increases, any point on the map is less constrained by elevational barriers, and thus LEC is mostly controlled by the abundance of sites with similar elevation in the simulated domain (i.e. the elevation frequency distribution as shown by the blue lines in Figs. 6c and 7b). Bertuzzo et al. (2016) show that for very large σ, LEC becomes insensitive to elevation and elevational gradients tend to flatten out.
4.2 LEC as a measure of biodiversity
Bertuzzo et al. (2016) have shown that when applied to real landscapes, the LEC metric predicts the α diversity simulated by the full metacommunity model well. Here we perform a similar analysis on our simulated landscape and estimate for both scenarios (uniform and orographic rain) the correlations between LEC and simulated local species richness (α diversity).
To do so, we use a zero-sum metacommunity model (Hubbell, 2001) in which local communities (N) are defined on an regular 2-D elevation mesh. The zero-sum assumption states that each local community is saturated at all times, which means that the total number of individuals n in a particular point is constant over time. In such a case and for any given time, the entire region is made of a population of N⋅n individuals, where N is equivalent to the number of points on the grid. The population number is made of different species which have different elevational niches. Each of these niches describes the competitive ability of particular species with elevation and is defined with a Gaussian function following Rybicki and Hanski (2013) (Bertuzzo et al., 2016):
with fi(z) the competitive ability of species i at elevation z. The elevation is the optimal elevation for which the competitive ability equals and σi is the niche width for species i as defined in Eq. (6). In addition to the zero-sum assumption, we suppose that all species have the same niche width σi=σ and maximum competitive ability as well as similar rates of dispersal, death, and fertility.
The ecological interactions between individuals follow the same approach as the one proposed in Bertuzzo et al. (2016). For any given time step and on each local community, a randomly chosen individual dies. Based on the zero-sum assumption, this dead individual needs to be, replaced and two strategies are possible. First, the offspring can come from the pool of individuals living in the vicinity of the local community (i.e. coming from the community itself or from one of the neighbouring ones). In this case, the chosen offspring is selected with a probability proportional to the values of fi(z) for all the individuals living in the vicinity of the local community (where z is the elevation of the considered local community). Second, the offspring is an individual belonging to a new species that does not already exist in the system. This second option is probabilistically selected at every time step and allows us to model both speciation and immigration from external communities (Hubbell, 2001; Chave et al., 2002). When a competitor from a new species enters the system its optimal elevation is drawn from a uniform distribution spanning twice the relief of the system to avoid edge effects (Bertuzzo et al., 2016).
This metacommunity model is applied on the elevation grids obtained from our landscape evolution model. In all simulations, the system is initially populated by one single species and is run until a statistically steady state is reached (∼104 generations, where a generation is N⋅n time steps). As in Bertuzzo et al. (2016), we consider periodic boundary conditions. The parameter fmax does not affect the system dynamics and is, without loss of generality, set to 1.
Figure 8 presents the relationships between calculated LEC and local species richness for the simulated landscapes at steady state. Overall, we have Pearson's correlations (R) of >70 %, which is considered to be relatively strong. These results are in agreement with the Bertuzzo et al. (2016) study and consolidate the idea that the LEC calculation is able to predict the α diversity well, especially in the case of uniform rainfall. Therefore, if LEC captures the regional variability of community diversity, one can deduce from its distribution the temporal and spatial evolution of species richness.
We provide in Fig. 9 a statistical summary of the LEC evolution for the two climatic scenarios. LEC predicts biodiversity peaks when the landscape reaches a geomorphological steady state, and maximum species richness is observed at mid-elevation (Fig. 9a – ) as observed in nature (Grytnes and Vetaas, 2002; Colwell et al., 2004; Brehm et al., 2007). Interestingly, when using a coupled climate–landscape model, windward-facing regions host higher biodiversities. During relaxation phase, when tectonic forcing ceases, the peak in species richness shifts from middle to low elevations (bottom graph in Fig. 9a). From the figure, we infer that mapped areas showing strong connectivity values (i.e. high LEC) within a given elevational band have a direct effect on the biodiversity they host (Rahbek, 1997; McCain, 2007). We also hypothesise that these areas represent regions of high local species richness as local communities can be gathered from a more diverse regional pool of species that are fit to live at a similar elevation (Romdal and Grytnes, 2007).
In summary, LEC calculation can be used to estimate in a simple way the temporal and spatial evolution of biodiversity distribution induced by changes in landscape structure over geological time. Based on LEC distribution at steady state, we found that species richness in mountainous landscapes reaches its peak at mid-elevation not only because this is where the maximum surface area available to species is located, but also because species can move up or down to accommodate climate warming or cooling, respectively.
There are many biologic and abiotic properties that drive species richness in mountainous landscapes (Hoorn et al., 2010, 2013; Giezendanner et al., 2019). For example, mountains provide species with a rich variety of environmental conditions over restricted surface areas (e.g. range of temperature, solar irradiation, wind exposure, moisture and rainfall, soils thickness, and composition to cite a few). In this study, we have limited our exploration to the role of morphological changes imposed by uniform tectonic, riverine, and climatic processes, and we found that these changes alone could explain to the 1st order the biodiversity found in mountainous regions.
4.3 Dynamics of geomorphically driven isolation
From the LEC values, one can derive geomorphically driven isolation by finding spatial regions of low LEC compared to their surrounding areas. Here, we apply an approach similar to a depression-filling algorithm (Barnes and Lehman, 2013) using the LEC values instead of the elevation. Isolated areas in that sense represent depressions or inwardly draining regions of the LEC map which have no outlet. To prevent coalescence of these isolated regions (i.e. small depressions within bigger ones), one can apply a threshold on the filling limit of the LEC map. Increasing such a threshold produces a smaller number of isolated regions but with larger areas.
In Fig. 10, we evaluate through time these areas of low LEC values (normalised LEC below 20 %) and analyse series of metrics: mean LEC, mean area, and mean elevation. Our results show that the competition between mountain building and erosion increases the number of isolated regions during the orogenic phase (Fig. 10a). Orographic precipitation fosters faster isolation than the uniform precipitation model (red line trends in Fig. 10a), especially on the steeper and drier leeward sides (rain shadow zones) of the drainage divides (Fig. 10b). For the uniform rain model, we observe a decrease in the average area of these isolated regions concomitant with catchment stabilisation and the emergence of the steady state (light grey line in Fig. 10a). For both precipitation scenarios, the cessation of tectonic activity at 5 Myr results in a rapid decrease in the number of isolated regions associated with an initial episode of augmentation of mean area that lasts no more than 250 000 years. We attribute this peculiar effect to the fast widening of valley high- and mid-elevation zones related to the adjustment of catchment morphology and induced by the advective slope retreat predominant across the valley heads (Pelletier, 2011). Hence, our model predicts an increase in isolation during the mountain-building phase and a decrease during the relaxation phase.
Based on the results from Sect. 4.2 and as high LEC regions represent – to the 1st order – places of higher biodiversity, we could take the argument further and propose that low LEC areas should potentially correspond to places where species get trapped into isolated refuges, and such ecological niches should favour speciation and endemism.
Not surprisingly, it has been found that the potential for isolation increases in more topographically diverse areas (Ali and Aitchison, 2014; Gillespie and Roderick, 2014). In recent years, several studies have shown that isolation induced by uplift and climatic processes plays a central part in shaping biodiversity (Wang et al., 2009; Craw et al., 2015), and the mapping of geomorphically driven low landscape connectivity could potentially highlight these areas of speciation hotspots. Therefore, the LEC calculation could provide an efficient way to identify and manage regions of higher conservation value.
Studies have predicted higher endemism at higher elevation as a consequence of increasing isolation with elevation (Steinbauer et al., 2016, 2018). We see a similar trend with high-elevation zones less connected to each other and statistically smaller than lower-elevation ones. However, our results also suggest that low-connectivity regions form fragmented patches that are distributed across the entire elevational range during an orogenic cycle (zmean in Fig. 10a). Several other parameters and mechanisms not accounted for in this study will promote spatial variations in speciation rates, including temperature and biotic interactions (Emerson and Kolm, 2005; Brown, 2014). Nevertheless, it shows that when considering the morphological complexity of mountainous landscapes over geological time, zones of low connectivity are not distributed continuously over elevation gradients and the associated isolated regions will likely be found at different elevations. Therefore, the prediction of increasing endemism with higher elevation might not always be the rule.
This paper presents a methodology to quantify landscape connectivity under tectonic and climatic forcing. Over geological timescales (millions of years), tectonics and surface processes conspire to transform flat regions into complex dissected mountainous landscapes. We used landscape elevational connectivity (LEC) to measure the connectivity between sites of similar elevations (Bertuzzo et al., 2016). We show that this abiotic parameter is able to account for up to 80 % of the α diversity predicted by a zero-sum metacommunity model (Hubbell, 2001; Rybicki and Hanski, 2013), and therefore it could potentially explain the 1st-order distribution of biodiversity found in mountainous regions (Ali and Aitchison, 2014; Steinbauer et al., 2018). In the future, we can easily improve the predictive capacity of this metric by integrating fitness parameters other than elevation.
From the LEC calculation, one can quantify the role of geomorphology in landscape connectivity and topography-driven isolation. As the periodicity of isolation and connection dictates evolutionary outcomes, understanding this dynamic might be used to test models of biological diversification and to understand species distribution and biodiversity patterns through time. Our results suggest that peaks in species richness for mountainous landscapes can be derived from the analysis of the dynamic organisation of geomorphic features as they evolved in response to tectonic, climatic, and erosion processes. These features might be inferred from physiographic metrics such as χ maps, in situ observations, river longitudinal profiles, or landscape evolution modelling. We also found that geomorphically driven isolation has the potential to increase rates of speciation over the entire elevational range. This is an important outcome that could potentially lead to the reassessment of the spatial viability of existing conservation sites and could help us improve future conservation planning, policy, and practice (Mokany et al., 2012).
The results and data presented and discussed in this paper were simulated using the Badlands model (https://badlands.readthedocs.io/en/latest/; Salles, 2016), and the connectivity maps were obtained from the bioLEC Python package (https://biolec.readthedocs.io/en/latest/; Salles and Rey, 2019). All the model outputs were visualised using the Python Matplotlib library for standard graphics, Seaborn data visualisation library for statistical graphics, and the Paraview software (v5.2.0; https://www.paraview.org, last access: 27 September 2019) from Kitware, Sandia National Labs, and CSimSoft for 3-D graphics and 2-D maps.
TS developed the code, designed the experiment, and contributed to output analysis and paper writing; PR contributed to output analysis and interpretation, as well as paper writing; EB developed the code and contributed to output analysis, interpretation, and paper writing.
The authors declare that they have no conflict of interest.
The authors acknowledge the Sydney Informatics Hub and the University of Sydney's high-performance computing cluster Artemis for providing the high-performance computing resources that contributed to the research results reported within this paper. We thank Phaedra Upton, the anonymous reviewer, and the journal editor for their comments that greatly improved the paper.
This paper was edited by Robert Hilton and reviewed by Phaedra Upton and one anonymous referee.
Ali, J. R. and Aitchison, J. C.: Exploring the combined role of eustasy and oceanic island thermal subsidence in shaping biodiversity on the Galápagos, J. Biogeogr., 41, 1227–1241, https://doi.org/10.1111/jbi.12313, 2014. a, b
Barnes, R. and Lehman, M.: Priority-Flood: An Optimal Depression-Filling and Watershed-Labeling Algorithm for Digital Elevation Models, Comput. Geosci., 62, 117–127, https://doi.org/10.1016/j.cageo.2013.04.024, 2013. a
Bertuzzo, E., Carrara, F., Mari, L., Altermatt, F., Rodriguez-Iturbe, I., and Rinaldo, A.: Geomorphic controls on elevational gradients of species richness, P. Natl. Acad. Sci. USA, 113, 1737–1742, https://doi.org/10.1073/pnas.1518922113, 2016. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s
Bonnet, S. and Crave, A.: Landscape response to climate change: Insights from experimental modeling and implications for tectonic versus climatic uplift of topography, Geology, 31, 123–126, 2003. a, b
Brehm, G., Colwell, R. K., and Kluge, J.: The role of environment and mid‐domain effect on moth species richness along a tropical elevational gradient, Global Ecol. Biogeogr., 16, 205–219, https://doi.org/10.1111/j.1466-8238.2006.00281.x, 2007. a, b
Carrara, F., Altermatt, F., Rodriguez-Iturbe, I., and Rinaldo, A.: Dendritic connectivity controls biodiversity patterns in experimental metacommunities, P. Natl. Acad. Sci. USA, 109, 5761–5766, https://doi.org/10.1073/pnas.1119651109, 2012. a
Castelltort, S. and Yamato, P.: The influence of surface slope on the shape of river basins: Comparison between nature and numerical landscape simulations, Geomorphology, 192, 71–79, https://doi.org/10.1016/j.geomorph.2013.03.022, 2013. a, b
Chave, J., Muller‐Landau, H. C., and Levin, S. A.: Comparing Classical Community Models: Theoretical Consequences for Patterns of Diversity, Am. Nat., 159, 1–23, 2002. a
Chen, A., Darbon, J., and Morel, J.-M.: Landscape evolution models: a review of their fundamental equations, Geomorphology, 219, 68–86, 2014. a
Craw, D., Upton, P., Burridge, C. P., Wallis, G. P., and Waters, J. M.: Rapid biological speciation driven by tectonic evolution in New Zealand, Nat. Geosci., 9, 140–144, https://doi.org/10.1038/ngeo2618, 2015. a, b, c
Craw, D., King, T. M., McCulloch, G. A., Upton, P., and Waters, J. M.: Biological evidence constraining river drainage evolution across a subduction-transcurrent plate boundary transition, New Zealand, Geomorphology, 336, 119–132, https://doi.org/10.1016/j.geomorph.2019.03.032, 2019. a
Dietrich, W. E., Reiss, R., Hsu, M. L., and Montgomery, D. R.: A process-based model for colluvial soil depth and shallow landsliding using digital elevation data, Hydrol. Process., 9, 383–400, 1995. a
Giezendanner, J., Bertuzzo, E., Pasetto, D., Guisan, A., and Rinaldo, A.: A minimalist model of extinction and range dynamics of virtual mountain species driven by warming temperatures, PLoS ONE, 14, e0213775, https://doi.org/10.1371/journal.pone.0213775, 2019. a
Grytnes, J. A. and Vetaas, O. R.: Species Richness and Altitude: A Comparison between Null Models and Interpolated Plant Species Richness along the Himalayan Altitudinal Gradient, Nepal, Am. Nat., 159, 294–304, https://doi.org/10.1086/338542, 2002. a, b
Hoorn, C., Wesselingh, F. P., ter Steege, H., Bermudez, M. A., Mora, A., Sevink, J., Sanmartín, I., Sanchez-Meseguer, A., Anderson, C. L., Figueiredo, J. P., Jaramillo, C., Riff, D., Negri, F. R., Hooghiemstra, H., Lundberg, J., Stadler, T., Särkinen, T., and Antonelli, A.: Amazonia Through Time: Andean Uplift, Climate Change, Landscape Evolution, and Biodiversity, Science, 330, 927–931, https://doi.org/10.1126/science.1194585, 2010. a, b
Kearney, M. and Porter, W.: Mechanistic niche modelling: combining physiological and spatial data to predict species' ranges, Ecol. Lett., 12, 334–350, https://doi.org/10.1111/j.1461-0248.2008.01277.x, 2009. a
Kessler, M., Kluge, J., Hemp, A., and Ohlemüller, R.: A global comparative analysis of elevational species richness patterns of ferns, Global Ecol. Biogeogr., 20, 868–880, https://doi.org/10.1111/j.1466-8238.2011.00653.x, 2011. a
Liu, C., Dudley, K. L., Xu, Z.-H., and Economo, E. P.: Mountain metacommunities: climate and spatial connectivity shape ant diversity in a complex landscape, Ecography, 41, 101–112, https://doi.org/10.1111/ecog.03067, 2018. a
MacArthur, R. H. and Wilson, E. O.: The Theory of Island Biogeography, Princeton University Press, Princeton, New Jersey, 1967. a
Martensen, A., Saura, S., and Fortin, M.-J.: Spatio-temporal connectivity: assessing the amount of reachable habitat in dynamic landscapes, Methods Ecol. Evol., 8, 1253–1264, https://doi.org/10.1111/2041-210X.12799, 2017. a
McCain, C. M.: Area and mammalian elevation diversity, Ecology, 88, 76–86, https://doi.org/10.1890/0012-9658(2007)88[76:AAMED]2.0.CO;2, 2007. a
Mokany, K., Harwood, T. D., Williams, K. J., and Ferrier, S.: Dynamic macroecology and the future for biodiversity, Glob. Change Biol., 18, 3149–3159, https://doi.org/10.1111/j.1365-2486.2012.02760.x, 2012. a, b
Mudd, S. M., Attal, M., Milodowski, D. T., Grieve, S. W. D., and Valters, D. A.: A statistical framework to quantify spatial variation in channel gradients using the integral method of channel profile analysis, J. Geophys. Res.-Earth, 119, 138–152, https://doi.org/10.1002/2013JF002981, 2014. a
Muneepeerakul, R., Bertuzzo, E., Rinaldo, A., and Rodriguez-Iturbe, I.: Evolving biodiversity patterns in changing river networks, J. Theor. Biol., 462, 418–424, https://doi.org/10.1016/j.jtbi.2018.11.021, 2019. a
Rahbek, C.: The elevational gradient of species richness: a uniform pattern?, Ecography, 18, 200–205, https://doi.org/10.1111/j.1600-0587.1995.tb00341.x, 1995. a
Roy, S. G., Koons, P. O., Upton, P., and Tucker, G. E.: The influence of crustal strength fields on the patterns and rates of fluvial incision, J. Geophys. Res.-Earth, 120, 275–299, https://doi.org/10.1002/2014JF003281, 2015. a
Salles, T.: Badlands: A parallel basin and landscape dynamics model, SoftwareX, 5, 195–202, https://doi.org/10.1016/j.softx.2016.08.005, 2016 (model available at: https://badlands.readthedocs.io/en/latest/, last access: 27 September 2019). a, b, c
Salles, T. and Hardiman, L.: Badlands: An open-source, flexible and parallel framework to study landscape dynamics, Comp. Geosci., 91, 77–89, 2016. a
Salles, T. and Rey, P.: bioLEC: A Python package to measure Landscape Elevational Connectivity, Journal of Open Source Software, 4, 1498, https://doi.org/10.21105/joss.01530, 2019 (code available at: https://biolec.readthedocs.io/en/latest/, last access: 27 September 2019). a, b
Salles, T., Ding, X., and Brocard, G.: pyBadlands: A framework to simulate sediment transport, landscape dynamics and basin stratigraphic evolution through space and time, PLOS ONE, 13, 1–24, https://doi.org/10.1371/journal.pone.0195557, 2018a. a, b
Salles, T., Ding, X., Webster, J. M., Vila-Concejo, A., Brocard, G., and Pall, J.: A unified framework for modelling sediment fate from source to sink and its interactions with reef systems over geological times, Sci. Rep., 8, 5252, https://doi.org/10.1038/s41598-018-23519-8, 2018b. a
Smith, B. T., McCormack, J. E., Cuervo, A. M., Hickerson, M. J., Aleixo, A., Cadena, C. D., Pérez-Emán, J., Burney, C. W., Xie, X., Harvey, M. G., Faircloth, B. C., Glenn, T. C., Derryberry, E. P., Prejean, J., Fields, S., and Brumfield, R. T.: The drivers of tropical speciation, Nature, 515, 406–409, https://doi.org/10.1038/nature13687, 2014. a, b, c
Solé, R. V. and Bascompte, J.: Self-Organization in Complex Ecosystems (MPB-42), Princeton University Press, Princeton, New Jersey, 2006. a
Steinbauer, M. J., Field, R., Grytnes, J., Trigas, P., Ah‐Peng, C., Attorre, F., Birks, H. J. B., Borges, P. A. V., Cardoso, P., Chou, C., Sanctis, M. D., de Sequeira, M. M., Duarte, M. C., Elias, R. B., Fernández‐Palacios, J. M., Gabriel, R., Gereau, R. E., Gillespie, R. G., Greimler, J., Harter, D. E. V., Huang, T., Irl, S. D. H., Jeanmonod, D., Jentsch, A., Jump, A. S., Kueffer, C., Nogué, S., Otto, R., Price, J., Romeiras, M. M., Strasberg, D., Stuessy, T., Svenning, J., Vetaas, O. R., Beierkuhnlein, C., and Gillespie, T.: Topography‐driven isolation, speciation and a global increase of endemism with elevation, Global Ecol. Biogeogr., 25, 1097–1107, https://doi.org/10.1111/geb.12469, 2016. a, b, c
Steinbauer, M. J., Grytnes, J.-A., Jurasinski, G., Kulonen, A., Lenoir, J., Pauli, H., Rixen, C., Winkler, M., Bardy-Durchhalter, M., Barni, E., Bjorkman, A. D., Breiner, F. T., Burg, S., Czortek, P., Dawes, M. A., Delimat, A., Dullinger, S., Erschbamer, B., Felde, V. A., Fernández-Arberas, O., Fossheim, K. F., Gómez-García, D., Georges, D., Grindrud, E. T., Haider, S., Haugum, S. V., Henriksen, H., Herreros, M. J., Jaroszewicz, B., Jaroszynska, F., Kanka, R., Kapfer, J., Klanderud, K., Kühn, I., Lamprecht, A., Matteodo, M., di Cella, U. M., Normand, S., Odland, A., Olsen, S. L., Palacio, S., Petey, M., Piscová, V., Sedlakova, B., Steinbauer, K., Stöckli, V., Svenning, J.-C., Teppa, G., Theurillat, J.-P., Vittoz, P., Woodin, S. J., Zimmermann, N. E., and Wipf, S.: Accelerated increase in plant species richness on mountain summits is linked to warming, Nature, 556, 231–234, https://doi.org/10.1038/s41586-018-0005-6, 2018. a, b, c
Tucker, G. and Hancock, G. R.: Modelling landscape evolution., Earth Surf. Proc. Land., 35, 28–50, 2010. a
van der Walt, S., Schonberger, J., Nunez-Iglesias, J., Boulogne, F., Warner, J., Yager, N., Gouillart, E., and Yu, T.: Scikit Image Contributors – scikit-image: image processing in Python, PeerJ, 2, e453, https://doi.org/10.7717/peerj.453, 2014. a
Wallace, A. R.: On the Zoological Geography of the Malay Archipelago, Journal of the Proceedings of the Linnean Society of London. Zoology, 4, 172–184, https://doi.org/10.1111/j.1096-3642.1860.tb00090.x, 1860. a
Wang, I. J., Savage, W. K., and Shaffer, H. B.: Landscape genetics and least‐cost path analysis reveal unexpected dispersal routes in the California tiger salamander (Ambystoma californiense), Mol. Ecol., 18, 1365–1374, https://doi.org/10.1111/j.1365-294X.2009.04122.x, 2009. a
Whipple, K. and Tucker, G.: Dynamics of the stream-power river incision model: implications for height limits of mountain ranges, landscape response timescales, and research needs, J. Geophys. Res., 104, 17661–17674, 1999. a, b
Whipple, K. X., Forte, A. M., DiBiase, R. A., Gasparini, N. M., and Ouimet, W. B.: Timescales of landscape response to divide migration and drainage capture: Implications for the role of divide mobility in landscape evolution, J. Geophys. Res.-Earth, 122, 248–273, https://doi.org/10.1002/2016JF003973, 2016. a