|Home | About | Journals | Submit | Contact Us | Français|
The species composition of plankton, insect and annual plant communities may vary markedly from year to year. Such interannual variability is usually thought to be driven by year-to-year variation in weather conditions. Here we examine an alternative explanation. We studied the effects of regular seasonal forcing on a multi-species predator–prey model consisting of phytoplankton and zooplankton species. The model predicts that interannual variability in species composition can easily arise without interannual variability in external conditions. Seasonal forcing increased the probability of chaos in our model communities, but squeezed these irregular species dynamics within the seasonal cycle. As a result, the population dynamics had a peculiar character. Consistent with long-term time series of natural plankton communities, seasonal variation led to a distinct seasonal succession of species, yet the species composition varied from year to year in an irregular fashion. Our results suggest that interannual variability in species composition is an intrinsic property of multi-species communities in seasonal environments.
Aquatic and terrestrial communities are often characterized by a complex waxing and waning of species driven by the seasonal cycle. Plankton communities show some regularity in the form of an annually recurring spring bloom. Yet the height, timing and species composition of the spring bloom often vary strongly from year to year (figure 1; see also Talling 1993; Harris & Baxter 1996; Smayda 1998; Philippart et al. 2000). Similar year-to-year variability in species composition has been observed in multi-species communities of insects (Wolda 1988; Raimondo et al. 2004), soil fauna (Giller & Doube 1994; Berg et al. 1998) and annual plants (Guo et al. 2002). Often, this interannual variability in species composition is attributed to year-to-year variation in weather conditions (i.e. exceptionally cold winters, wet springs or hot summers). However, mathematical models (Hastings & Powell 1991; Vandermeer 1993; Huisman & Weissing 1999; Brose 2008) and laboratory experiments (Becks et al. 2005; Graham et al. 2007; Benincà et al. 2008) have shown that interactions between species may generate striking chaotic fluctuations in species abundances even without external forcing. One might therefore hypothesize that interannual variability in species composition may not require year-to-year variation in weather conditions or other external factors. Interannual variability could be an intrinsic property of multi-species communities in seasonal environments. To investigate this hypothesis, it is interesting to assess to what extent complex dynamics in multi-species communities are modified by the seasonal cycle.
Time series (bi-weekly averages) of marine phytoplankton species in the Marsdiep tidal inlet between the North Sea and Wadden Sea, The Netherlands, from January 1985 to December 1990. (a) Leptocylindrus minimus, dark blue line; Phaeocystis globosa, green ...
The effects of regular seasonal variation on population models of two or three interacting species have been studied extensively (Kot & Schaffer 1984; Doveri et al. 1993; Rinaldi et al. 1993; Steffen et al. 1997; Huppert et al. 2005). These studies have shown that periodically forced populations can display a rich repertoire of dynamical behaviours, including simple and complex periodic cycles, quasi-periodicity and chaos (Rinaldi et al. 1993; King & Schaffer 1999; Vandermeer et al. 2001). However, the parameter space in which chaotic behaviour occurs is usually small. Typically, the population dynamics show repeatable patterns. Slow-growing species may fluctuate on seasonal or multi-annual time scales, as exemplified by the famous cycles of voles and lemmings at northern latitudes (Stenseth 1999; Turchin 2003). Fast-growing species such as bacteria or plankton may display multiple ups and downs per year. The frequency of the population fluctuations can be remarkably persistent as a result of frequency locking (e.g. Scheffer et al. 1997; Vandermeer et al. 2001). Seasonal forcing tends to ‘lock’ the frequency of population oscillations, such that populations oscillate at the same frequency as the seasonal cycle or integer multiples of it.
While many theoretical studies have examined effects of seasonality on model systems of only a few species, seasonal forcing of multi-species communities has received surprisingly little theoretical attention (Ebenhöh 1992). Yet bacterial, plankton and insect communities may contain tens, hundreds and sometimes even thousands of species (Hutchinson 1961; Erwin 1982; Irigoien et al. 2004; Venter et al. 2004). Generally speaking, multi-species models display more complex dynamics than models with only two or three species (May 1973; Ellner & Turchin 1995; Huisman & Weissing 2001). From a conceptual perspective, multi-species food webs can be interpreted as systems with several interacting oscillations (e.g. several predator–prey cycles). Coupled oscillations are known to generate complex dynamics, including chaos (Hastings & Powell 1991; Vandermeer 1993, 2004; Huisman & Weissing 2001; Benincà et al. 2008). The prevalence of complex dynamics is of interest, because these non-equilibrium dynamics may help to sustain the biodiversity of natural communities (Armstrong & McGehee 1980; Huisman & Weissing 1999; Brose 2008) and also because complex dynamics can induce regime shifts in ecosystems, with important implications for their management (Scheffer et al. 2001; Hsieh et al. 2005; Ives et al. 2008).
Here, we study the effect of seasonal forcing on the dynamics of a multi-species predator–prey model, using phytoplankton and zooplankton as our model organisms. We use this model to assess to what extent a regular seasonal cycle will modify chaos in multi-species communities. Our results show that regular seasonal forcing can promote year-to-year variability in species composition. In addition, our results suggest that this interannual variability in species composition is affected by ecosystem productivity in a manner analogous to Rosenzweig's (1971) classical paradox of enrichment.
We study a minimal model that is sufficiently complex to investigate the impact of seasonal forcing on multi-species communities, yet sufficiently simple to produce generic insights. The model is based on a straightforward multi-species version of the classic Rosenzweig–MacArthur predator–prey model (Rosenzweig & MacArthur 1963; Vandermeer 1993; van Nes & Scheffer 2004), extended with seasonal forcing (Rinaldi et al. 1993; Scheffer et al. 1997). In our interpretation, the model represents a plankton community, although our findings can probably be generalized to other multi-species communities in which organisms have fast growth rates and short generation times compared with the length of the growing season (e.g. microbial food webs, soil fauna, tropical insects). Let Pi and Zk denote the population abundances of phytoplankton species i and zooplankton species k, respectively. Then the model reads
The phytoplankton species (equation (2.1)) grow logistically with maximum specific growth rates ri, carrying capacities Ki and competition coefficients αij to describe competition between species. The phytoplankton species are consumed by zooplankton species, as described by a multi-species functional response (of Holling type II) with a fixed half-saturation constant Hk and maximum grazing rate gk. Selective predation (Chesson 1978) is introduced through the selectivity coefficient Sik of zooplankton species k for phytoplankton species i and it can take values between 0 and 1, indicating the preference of the predator for its prey (van Nes & Scheffer 2004). The factor u accounts for small levels of immigration and is introduced to reduce the probability of heteroclinic cycles. Heteroclinic cycles are considered to be mostly biologically irrealistic, since species reach extremely low population abundances during these cycles without becoming extinct (May & Leonard 1975). The zooplankton species (equation (2.2)) grow on the consumed phytoplankton with an assimilation efficiency ek, suffer a mortality rate mk and immigrate at a small rate u similar to the phytoplankton.
Many biological parameters are sensitive to seasonal forcing. One might thus argue that seasonal forcing should be applied to all model parameters, perhaps with different parameters affected by seasonality in different ways depending on the species. However, this would yield a rather complex model, while we aim at a simple model that captures the essence of multi-species dynamics in a seasonal environment. Accordingly, we choose a simple way to incorporate seasonal forcing following earlier contributions (Doveri et al. 1993; Scheffer et al. 1997). In particular, seasonal variation in temperature and light conditions has a major impact on the growth rates and mortality rates of plankton species (Raven & Geider 1988; Litchman & Klausmeier 2001) and on the seasonal development of total plankton biomass (Sommer et al. 1986; Longhurst 2006). We therefore assume that seasonal fluctuations in the species' growth rates, mortality rates and carrying capacity (ri, Ki, gk, mk) can be described by a sinusoidal forcing function σ(t) (equation (2.3)), which can be interpreted as the environmental forcing imposed by seasonal variation in temperature or light. Factor a determines the amplitude of the seasonal forcing (Rinaldi et al. 1993) and takes values between 0 and 1. The cosine function is chosen to produce maximum rates in summer and minimum rates in winter (t = 0 is 1 January), and the period is set to 365 days (Scheffer et al. 1997).
We parameterized the model for 10 competing phytoplankton species (i = 1,…,10) grazed by 6 zooplankton species (k = 1,…,6). The parameter values assigned to the different species were selected from the ranges indicated in table 1, which are representative for plankton communities (Scheffer et al. 1997; Reynolds 2006). Phytoplankton intraspecific competition was set to unity (αii = 1 for all i), while the interspecific competition coefficients (αij) were drawn randomly from the interval (0.5, 1.5). Differences in grazing rate were introduced through the selectivity coefficients Sik, which were drawn randomly from the interval (0, 1) to create a food web of generalists where predators use prey species with average selectivity Savg,k = ∑iSik/10 = 0.5. We assumed that the carrying capacities of all phytoplankton species are equal (i.e. Ki = K for all species i), following the rationale that K is an environmental parameter reflecting the local nutrient and light conditions.
Parameter ranges used in the model simulations. The exact parameter values of each individual simulation presented in the figures are given in the electronic supplementary material, appendix S3.
We analysed the model without seasonal forcing (a = 0) and with seasonal forcing (0 < a < 1); the time-averaged parameter values in the model simulations with seasonal forcing were equal to the fixed parameter values used in the model simulations without seasonal forcing. We investigated the model communities at different levels of productivity (K = 2, 5, 10, 20, 50 mg l−1), to compare the species dynamics in a range from oligotrophic to eutrophic conditions.
We assessed how frequently the model communities displayed chaos and we calculated the corresponding values of the Lyapunov exponent by assembling 100 randomly generated model communities for every model scenario that we investigated. For this purpose, the parameter values of the 6 predator and 10 prey species in each model community were drawn randomly from uniform distributions covering the ranges indicated in table 1, and the initial biomasses of the species were drawn randomly from the interval (0, 10 mg l−1). The model communities were first simulated for 1000 years to ensure that the population dynamics had reached an attractor. Thereafter, we continued the model simulation for another 40 years, calculated the Lyapunov exponent and determined the nature of the attractors as stable, simple periodic (period-one limit cycles), complex periodic (multiple-period cycles), quasi-periodic and chaotic. The Lyapunov exponent quantifies the rate of exponential divergence (or convergence) of nearby trajectories (Strogatz 1994; Sprott 2003). A positive Lyapunov exponent indicates chaos, and its magnitude is a measure of the system's sensitivity to initial conditions. Our calculation of the Lyapunov exponent is explained in the electronic supplementary material, appendix S1. We used visual inspection and Poincaré maps as additional methods to verify the computed nature of the attractors or to check for undetermined cases. All simulations were carried out in MATLAB using our software package Grind (freely available at http://www.dow.wau.nl/aew/grind).
Without seasonal forcing, the model predicts various kinds of asymptotic regime, including stable equilibria (figure 2a), simple limit cycles (figure 2b), complex periodic cycles (figure 2c) and chaos (figure 2d). At first sight, seasonal forcing seems to have little influence on the dynamical repertoire of the model. With seasonal forcing, the model also displays simple limit cycles (figure 2e), complex periodic cycles (figure 2f) and chaos (figure 2g). However, a closer look reveals differences between the model behaviour with and without seasonal forcing. With seasonal forcing, the periodic solutions are ‘locked’ within the seasonal cycle: the same pattern repeats each year (figure 2e) or after some years (figure 2f). In addition, the model can also produce quasi-periodic cycles, where solutions are entrained within the seasonal cycle yet never repeat themselves as they slightly shift phase every year. Chaotic communities seem to experience similar seasonal patterns. However, the fluctuations of phytoplankton and zooplankton species in chaotic communities remain irregular even when entrained in a regular seasonal environment (figure 2g).
Community dynamics predicted by the model. The two top panels indicate the nature of the environmental forcing. Without seasonal forcing, the model produces (a) stationary equilibria, (b) simple cycles, (c) complex periodic cycles or (d) chaotic dynamics. ...
These dynamics can be illustrated by Poincaré maps sampling the model communities once per year for many consecutive years. Model communities with a periodicity of 1 year return to exactly the same species composition year after year, which appears as a single point on the Poincaré map. Communities with a periodicity of N years produce N points on the Poincaré map, quasi-periodicity produces a closed curve (figure 3a), while chaos produces a complex fractal structure (figure 3b).
Poincaré maps with annual snapshots of the model community collected over many years. More specifically, the maps plot the biomasses of two plankton species sampled from the model community at the 1st of January of each year for 100 000 years. ...
Many of the model communities exposed to seasonal forcing displayed chaos with remarkable synchronization patterns at the species level (figure 4). The species fluctuations are irregular, yet these irregular fluctuations are squeezed within the seasonal cycle. As a consequence, species enter the winter season in different proportions, and this affects the species composition of the next spring bloom. For instance, figure 4a shows a typical phytoplankton spring species. It reaches peak abundance in March, although its peak abundance varies from year to year, and some years it does not peak in spring at all. Figure 4b shows another phytoplankton species from the same plankton community. This species could be called a typical summer species. It is present every summer. However, some years it peaks twice, with a first peak in May–June and a second smaller peak in September. In other years, it peaks in September only. The zooplankton species show similar seasonal patterning. For instance, some zooplankton species are mainly present in winter (figure 4c), while others dominate during the summer period (figure 4d). The example in figure 4d is particularly interesting. In some years, this zooplankton species shows little variability from March to September, while in other years, it fluctuates wildly during the same period. Accordingly, the species composition in our model communities shows distinct patterns of seasonal organization, but with strong year-to-year variability.
Year-to-year variability in the population dynamics of (a) phytoplankton species 9, (b) phytoplankton species 6, (c) zooplankton species 5 and (d) zooplankton species 3. All four species are from the same model community, simulated for a total period ...
Which species traits and environmental conditions are responsible for the widespread chaotic dynamics in our model communities? A complete answer to this question is beyond the scope of this paper. However, some insight can be obtained by modifying the model assumptions systematically. This shows that more than 50 per cent of the model simulations produced chaos when using our default parameter settings (table 2, first row). The occurrence of chaos was not very sensitive to the relative magnitude of intraspecific versus interspecific phytoplankton competition (table 2). In contrast, modifying zooplankton predation had a striking effect on the occurrence of chaos. When zooplankton was removed from the model, very few simulations showed chaotic dynamics and they did so only under seasonal forcing (table 2). Similarly, inefficient zooplankton grazing and specialist zooplankton reduced the occurrence of chaos. This shows that predator–prey oscillations, and the nature of predation, played a key role in the generation of complex dynamics in our model communities.
Occurrence of chaos in our simulated communities under different model assumptions. The first row (in italics) shows the percentage of chaotic communities predicted by the reference model used in our study, both without seasonal forcing (a = 0) and with ...
Productivity also had a clear effect on the occurrence of chaotic dynamics. At low productivity (K = 2 mg l−1), stationary dynamics prevailed in constant environments, simple periodic dynamics prevailed in seasonal environments and chaos occurred only in a few model communities with strong seasonal forcing (table 2; see also appendix S1 in the electronic supplementary material). Chaos was widespread at intermediate productivity (K = 5 and K = 10 mg l−1). At high productivity (K = 20 and K = 50 mg l−1), the occurrence of chaos declined and the population dynamics often shifted to periodic cycles in both constant and seasonal environments.
In all cases summarized in table 2, seasonal forcing increased the occurrence of chaos. To investigate this aspect in further detail, we estimated whether the amplitude of seasonal forcing affected the occurrence of chaos in our model communities (figure 5). We focused on the intermediate productivities (K = 5 and K = 10 mg l−1). At K = 5 mg l−1, the amplitude of seasonal forcing increased the occurrence of chaos (figure 5a; linear regression: R2 = 0.45, N = 11, p = 0.024). At K = 10 mg l−1, mild forcing (0.1 < a < 0.4) caused a slight increase in the probability of chaos, but when the amplitude of seasonal forcing was further increased (a > 0.6), the probability of chaos declined (figure 5b; quadratic regression: R2 = 0.66, N = 11, p = 0.013). We further explored the predictability of these communities by calculating their Lyapunov exponents. A positive Lyapunov exponent indicates chaos. The inverse value of the Lyapunov exponent is often used as a simple metric of the predictability of chaotic systems (Strogatz 1994). In those simulations that displayed chaotic dynamics, the magnitude of the Lyapunov exponent was not affected by the amplitude of seasonal forcing (figure 5c,d). This indicates that the predictability of the chaotic plankton communities was neither enhanced nor reduced by a stronger seasonality. However, the median values of the Lyapunov exponents were significantly higher at K = 10 mg l−1 than at K = 5 mg l−1 (figure 5c,d; t-test: t = −3.77, d.f. = 20, p < 0.002), which indicates that the predictability of the model communities was affected by productivity.
(a,b) Relative frequency at which randomly generated model communities display chaos, plotted as a function of the amplitude of seasonal forcing. Results are shown for model communities grown at two productivities: (a) K = 5 mg l−1 and (b) K = ...
Our model results show that interannual variability in species composition can be an intrinsic property of multi-species communities in a seasonal environment. It does not require year-to-year variability in weather conditions. In many simulations, the timing and abundances of different plankton species varied strongly, both within years and among years. Some species peaked only once per year, while others peaked two or three times; some species were present every year, while other species popped up only occasionally (figure 4). An often-invoked and seemingly straightforward intuitive explanation for this interannual variability in species composition is that winter ‘resets’ population densities, whereas stochastic variation in weather conditions during spring and summer determines interannual differences in community composition. However, this idea of winter resetting is obviously an oversimplification. Each autumn, species enter the winter season in different proportions. Thus, the species composition from the previous autumn affects the species composition of the next spring bloom. Our model results show that this mechanism of seasonally entrained chaos can easily create interannual variability in species composition without invoking year-to-year differences in external environmental conditions.
Interannual variability in species composition does not imply that the species composition varies at random. The model predicts temporal organization of the species, even if their population dynamics is chaotic, in the form of a seasonal succession. Some species occur mainly in spring, while others dominate in summer (figure 4). Which mechanisms are responsible for this seasonal pattern? For instance, what makes a species a typical spring species? The model assumes similar thermal physiologies for all species (i.e. they are all forced by the same function σ(t)). Therefore, interspecific differences in thermal physiology or other species-specific seasonal cues cannot explain the seasonal succession predicted by the model. Instead, the spring species in figure 4a peaks in March, because its main predators (e.g. zooplankton species 5) have just declined while its key competitors (phytoplankton species 6) are still low in abundance. Apparently, seasonal forcing locks the species interactions, such that the only window of opportunity for this species is restricted to the spring. This illustrates that seasonal succession can be an emergent property of the underlying community dynamics, in which species are sorted according to their positions in the complex network of species interactions.
The classical predator–prey model of Rosenzweig and MacArthur (1963), which provided the point of departure in our study, is clearly a major simplification of reality. For instance, the model does not specify the underlying mechanisms of phytoplankton competition for nutrients and light (Tilman 1977; Huisman et al. 1999), ignores induced defences and other forms of phenotypic plasticity (Vos et al. 2004; Stomp et al. 2008), does not detail the population structure and life history of plankton species (De Roos et al. 1992; Nelson et al. 2005), neglects the potentially stabilizing effect of planktivorous fish on zooplankton dynamics (Scheffer et al. 1997; Gliwicz & Wrzosek 2008), does not take into account species-specific adaptations to the seasonal cycle (such as resting stages; Marcus & Boero 1998), and ignores many other factors that may play a role in the population dynamics of natural plankton communities. However, the model does describe the essence of multi-species predator–prey interactions. As such, it provides the core of more complex plankton models widely used in aquatic ecology and oceanography (e.g. Fasham et al. 1990; Franks 2002; Lima et al. 2002). It is therefore interesting to assess to what extent the key qualitative predictions of the model, most notably the interannual variability in species composition, are consistent with empirical data.
The predicted patterns of seasonal organization with interannual variability at the species level are well in line with observations from real plankton communities. A closer look at the time series of the Dutch coastal zone in figure 1 reveals typical spring species like the diatom Asterionella kariana (figure 1b). If present, this species blooms in March, although its peak abundance varies from year to year, and some years it does not peak in spring at all. The prymnesiophyte Phaeocystis globosa, a nuisance species that can leave large layers of foam on the beach, reaches its maximum in late spring or early summer, and in some years has a smaller second peak in late summer (figure 1a). The diatoms Rhizosolenia imbricata and Guinardia delicatula bloom mainly in the period between June and August and can thus be called summer species. Whereas R. imbricata typically blooms only once per year, G. delicatula can display several peaks per year. All species in this time series show striking year-to-year variability in timing and/or peak abundance.
Similar examples of interannual variability have been documented in many studies. Maberly et al. (1994) report considerable year-to-year variation in the timing and peak abundance of the diatom Asterionella formosa in a 45-year time series in Lake Windermere, UK. Smayda (1998) recognized different patterns of species variability in a 37-year time series of the plankton in Narragansett Bay, USA. For instance, the diatom Asterionellopsis glacialis displayed episodic irregular blooms, while the diatom Thalassiosira nordenskioeldii peaked at 5-year intervals. Interannual variation in species composition can have dramatic consequences. The spring bloom in the Kattegat between Denmark and Sweden is usually dominated by diatom species. In the late spring of 1988, however, the prymnesiophyte Chrysochromulina polylepis produced a major bloom with severe toxic effects on higher organisms, including fish, molluscs, ascidians and cnidarians (Nielsen et al. 1990; Lekve et al. 2006). Since the 1988 event, large-scale blooms of Chrysochromulina have not returned in the area. Numerous other studies have described similar patterns of interannual variability in plankton community composition (Reynolds & Bellinger 1992; Talling 1993; Harris & Baxter 1996; Arhonditsis et al. 2004; Huisman et al. 2006; Valdés et al. 2007; Smetacek & Cloern 2008). In many of these case studies, the underlying causes for the observed interannual variability were not apparent.
Not all model simulations produced interannual variability. For the same environmental setting (i.e. same values of a and K), some simulations generated chaos whereas other simulations with different species combinations generated simple periodic solutions (figure 2). This seems in line with natural plankton communities, where some time series display less interannual variability than other time series (Smetacek & Cloern 2008). For example, phytoplankton in Lake Kinneret showed little interannual variability over a period of 20 years, with dinoflagellate species dominant in late winter–spring and small chlorophytes in summer–autumn (Berman et al. 1992). In a brackish lagoon along the Baltic Sea, a copepod and polychaete species showed very stable seasonal succession over a period of 22 years, while rotifers displayed high interannual variability (Feike & Heerkloss 2008). Thus, it would be interesting to investigate to what extent interannual variability of plankton communities depends on community composition. For instance, would plankton communities dominated by buoyant cyanobacteria be more regular than communities dominated by diatom species?
One might think that seasonality would enhance the predictability of multi-species communities. Indeed, our results show that seasonal forcing enables temporal organization of the species. Species are entrained within the seasonal cycle and become dominant during the specific periods in the year that match their highest growth potential (figure 4). In this sense, seasonality generates recurrent patterns in species composition. Yet the seasonal cycle also interferes with intrinsic species interactions, which can have both stabilizing and destabilizing effects (e.g. Rinaldi et al. 1993; King & Schaffer 1999). In particular, mild seasonal forcing increases the likelihood of chaotic dynamics (figure 5a,b), whereas strong seasonal forcing may lead to synchronization of the species dynamics (figure 5b). Strikingly, for chaotic communities, the magnitude of the Lyapunov exponent is not affected by the strength of seasonal forcing (figure 5c,d). Thus, the predictability of the species trajectories is independent of the strength of seasonal forcing, which suggests that seasonality per se does not necessarily affect the time horizon for accurate prediction of changes in plankton community structure.
Our model predicts that ecosystem productivity affects the nature of the species fluctuations. This result can be explained by Rosenzweig's (1971) ‘paradox of enrichment’ (see appendix S2 in the electronic supplementary material for a complete discussion). In short, this classic work showed that stable predator–prey systems start to display oscillations when productivity is increased beyond a certain threshold value. In our multi-species communities, with 10 phytoplankton and 6 zooplankton species, there are 60 predator–prey pairs. With increasing productivity, many of these predator–prey pairs will start to oscillate, each with its own characteristic frequency. It is well known that such systems of coupled nonlinear oscillators have a strong tendency to generate chaos (Vandermeer 1993, 2004; Huisman & Weissing 2001; Benincà et al. 2008), which explains the widespread occurrence of chaos in model communities at intermediate productivity (table 2). Interestingly, we found that many of these chaotic predator–prey fluctuations coalesced to periodic cycles at high productivity (i.e. in hypertrophic environments). These model predictions suggest that changes in productivity (e.g. through human-induced eutrophication) are likely to alter patterns of interannual variability in multi-species communities.
We thank G. Sugihara, K. D. Jöhnk, T. Fukami and the anonymous reviewers for their constructive comments on earlier versions of the manuscript, and C. Sprott for his valuable advice. C. G. Cadée and J. M. van Iperen are acknowledged for the plankton counts of the Marsdiep. V. D. was supported by a Bodossakis Foundation Fellowship. E. B. and J. H. were supported by the Earth and Life Sciences Foundation (ALW), which is subsidized by the Netherlands Organization for Scientific Research (NWO).