Do latitudinal gradients exist in New Zealand stream invertebrate metacommunities?

Introduction

The latitudinal diversity gradient is among the most well-known patterns in ecology (Hillebrand, 2004; Jocque et al., 2010). While general patterns of increasing richness from the poles to the equator are common, there are many exceptions (Gaston & Blackburn, 2000; Hillebrand, 2004; Heino, 2011). The potential mechanisms behind this gradient are broad, whether non-biological (e.g. mid-domain effect hypothesis; Colwell & Lees, 2000), ecological (e.g. species-energy hypothesis; Currie, 1991), or evolutionary/historical (e.g. evolutionary rate and effective evolutionary time hypotheses; Mittelbach et al., 2007), but incorporating variation among local communities can provide additional insight (Qian & Ricklefs, 2007; Qian, Badgley & Fox, 2009; Leprieur et al., 2011; Astorga et al., 2014). As a key mechanism behind the latitudinal diversity gradient, climate increases in harshness with increasing latitude (Stevens, 1989). However, many other factors influence local climate including island size and the level of isolation. Isolated oceanic islands, for instance, have lower seasonality and predictability than continental locations at similar latitudes (Tonkin et al., 2017a; Fig. 1). Jocque et al. (2010) argue that a shift in climatic stability with latitude drives a dispersal–ecological specialisation trade-off at the metacommunity level, producing gradients in dispersal ability, ecological specialisation, range size, speciation, and species richness. In particular, increased temporal variability in environmental conditions promotes increased dispersal ability of organisms (Dynesius & Jansson, 2000; Jocque et al., 2010).

Community differences attributable to latitude are therefore likely to be driven by underlying metacommunity processes. Four metacommunity archetypes have been synthesised to summarise the relative roles of local (niche) and regional (dispersal) processes in community assembly (Leibold et al., 2004; Holyoak et al., 2005; Leibold & Chase, 2018): neutral, patch dynamics, species sorting, and mass effects. What remains to be tested, however, is the influence that latitude has on the roles of different metacommunity processes (Jocque et al., 2010). In a testable hypothesis, Jocque et al. (2010) predicted a stronger role of dispersal limitation in the tropics accompanied by a shift to more species sorting and mass effects with increasing latitude.

Situated at mid-latitudes, New Zealand comprises a series of islands spanning a large latitudinal gradient. With a climate reflecting its oceanic position, rainfall (Fig. 1) and river flow regimes are typically unpredictable (Winterbourn, Rounick & Cowie, 1981; Winterbourn, 1995). Although most streams tend to be perennial, the high variability in rainfall (Heine, 1985) produces considerable variation in flows, with frequent, but typically short-duration, spates and floods (Duncan, 1987). Coupled with their flashy flow regimes comes a lack of seasonality in some food resources because of a predominantly evergreen flora (Winterbourn, Rounick & Cowie, 1981; Thompson & Townsend, 2000). These factors, combined with its highly dynamic geological history, making the country particularly sensitive to sea-level fluctuations during the Quaternary, ultimately lead to a largely generalist, opportunistic, and seasonally asynchronous stream fauna adapted to coping with these harsh conditions and climatic unpredictability (Winterbourn, Rounick & Cowie, 1981; Winterbourn, 1995; Thompson & Townsend, 2000). Most notably, New Zealand streams feature a predominance of endemic genera, invertebrates with poorly synchronised and flexible life histories, and a predominance of non-specialist ‘collector-browser’ species (Winterbourn, 1995). Consequently, New Zealand stream communities provide an interesting test case for investigating latitudinal controls on community structure.

To test a series of hypotheses related to metacommunity structuring across a broad latitudinal gradient, we explored gradients of stream invertebrate metacommunity structure (spatial structuring and environmental filtering) spanning the length of New Zealand’s two largest islands (∼1,300 km). As a secondary exploration, we examined the best-fit idealised ‘metacommunity types’ assigned through the Elements of Metacommunity Structure framework (EMS; Leibold & Mikkelson, 2002). To further disentangle the role of dispersal, we deconstructed assemblages into strong and weak dispersers. Doing so can be fruitful for exploring processes behind latitudinal diversity gradients (Kneitel, 2016). Taking a multi-faceted approach across latitudinal gradients allows for identifying complementary patterns in factors shaping metacommunities, compared to local community structure, advancing our understanding of how communities assemble in such dynamic landscapes.

We tested the following primary hypotheses based on the predictions of Jocque et al. (2010): Metacommunities are primarily structured by environmental variables (in line with the species sorting archetype; H_1a) and spatial variables increase in importance from north to south (representing increasing dispersal and in line with the mass effects archetype; H_1b). The alternative to this hypothesis (H₁A) is that, given the highly dynamic nature of New Zealand streams (Winterbourn, Rounick & Cowie, 1981), they are so unpredictable (at different stages of post-flood succession) that metacommunity structuring is context dependent from region to region. Because environmental heterogeneity and the spatial extent of metacommunities are important regulators of the relative strength of species sorting compared to dispersal limitation and surplus (both of which should increase the spatial signature in the metacommunity) (Heino et al., 2015b), we also explored the influence of these factors on observed patterns. Using the deconstructed dispersal groups, we tested the secondary hypothesis, based on the predictions of Jocque et al. (2010) (H₂), that strong dispersers increase from north to south. The EMS analysis was used as an additional exploratory analysis, thus we did not form any specific hypotheses.

Methods

Study sites

We used data previously collected (Astorga et al., 2014) from 120 streams in eight regions (15 sites in each region), spanning a latitudinal gradient of 12° (Fig. 2). Values of regional γ and β diversity, and mean α diversity are reported in Astorga et al. (2014). These eight datasets span across the five biogeographic regions of the New Zealand mainland (Di Virgilio et al., 2014): two in northern North Island, two in southern North Island, one in central New Zealand, two in mid-South Island, and one in southern South Island. Site selection followed a series of criteria, outlined in the following sentences, to minimise differences between regions. Streams were sampled primarily in protected areas (National or State Forest Parks) and were restricted to those with maximum of 14% exotic forestry and 30% pasture in the upstream catchment. All sites had a minimum intact riparian buffer of 50 m (Freshwater Ecosystems of New Zealand (FENZ)) (Leathwick et al., 2010) and were selected in proportion to FENZ classes in regions. Sites were restricted to <7 m wide headwater streams (order 1–3), with similar aspect and with a rocky substrate, and sampling was confined to the riffle zone. Almost all of New Zealand’s landmass belongs to the temperate oceanic (Cfb) climate zone (Peel, Finlayson & McMahon, 2007). Although the large majority of our sampling sites were geographically situated in this zone, some of the South Island sites (e.g. in Fiordland and Arthur’s Pass) likely fell on the border of temperate oceanic, subpolar oceanic (Cfc) and tundra (ET) zones. Nevertheless, all sites were situated in areas without shortage of rainfall, which has been described as ‘plentiful’ but temporally variable in New Zealand (Heine, 1985), although there are more arid regions in eastern zones such as Hawke’s Bay and central Otago. Therefore, all sites had permanent flow and the large majority of streams were runoff fed.

Results

The Fiordland and Northland metacommunities had the greatest spatial extents (Fig. 2E), but there was little difference in environmental heterogeneity between the regions (Fig. 2F). The gradient in environmental conditions was weak across the eight regions, with a low percentage of variance explained (37%) by the first two principal components (Fig. 2B), and no variables contributing more than 15% to either of the first two components. Invertebrate communities differed significantly between the eight regions, with a clear latitudinal trend in assemblage structure (PERMANOVA: F_7,112 = 7.30, R² = 0.313, P = 0.001; Fig. 2C). Regional richness tended to be highest at the north of each island and decline towards the southern zones (Fig. 2D), as demonstrated in Astorga et al. (2014). The regional pool of most regions were well sampled. However, Kahurangi did not reach a clear asymptote and had the steepest species accumulation curve. Moreover, the North Island regions’ curves tended to be less steep compared to those in the South Island. However, Chao’s estimated values did not differ in a systematic manner, with differences between sampled and projected richness not being consistently higher in the South Island.

Metacommunity structuring and the role of dispersal

There was no gradient with latitude in the relative importance of environmental or spatial control for all species combined and for individual dispersal groups (Fig. 3) suggesting H₁ can be rejected. The influence of spatial extent and its interaction with dispersal ability did not resolve this lack of pattern in the relative role of spatial or environmental components in the variation partitioning models (Figs. 2–4). Finally, contrary to H₂, the ratio of strong to weak dispersers decreased from north to south (Fig. 4).

When considering all species together, only three of the eight regions were significantly structured by both environmental and spatial components together, and thus could be considered for variation partitioning (Fig. 3). In the deconstructed dispersal group datasets, only one of the eight regions had combined significant environmental and spatial components. Environmental control was more commonly important than spatial in structuring both strong and weak disperser metacommunities. Northland exhibited no spatial or environmental structure for any of the datasets.

Considering all models (including those assigned 0% explained), environmental variables explained more of the variation when the whole community was considered (mean adj. R² = 0.134; 13.4% variance explained) compared to breaking into high (7.1%) and low (4.8%) dispersal ability groups. This result was particularly evident for certain regions, such as Westland, which could be explained well when considering the full community (strongest model), but not for the dispersal groups. However, strong dispersers had on average higher adj. R² values (adj. R² = 0.191; 19.1% explained) when only considering the significant models, than all combined (18.0%) or weak dispersers (9.6%). Spatial variables explained less of the variation in community structure than environmental, when non-significant models were included (adj. R²—all: 0.047; high: 0.049; low: 0.054) but not when only considering significant models (adj. R²—all: 0.126; high: 0.200; low: 0.143).

Forward-selected environmental variables were highly variable in the RDA models, with no particular variable consistently important across the eight metacommunities (Table 2; Table S1 in Appendix S1).

Table 2:

Forward-selected environmental variables used in the variation partitioning analysis when a global RDA model was significant.

Subset	Region	F	P	Variables
All	U	2.57	0.001	Temp, pH
All	E	2.96	0.001	OHCov, Elev, SI, Depth
All	K	2.25	0.001	Cond, OHCov
All	A	2.64	0.026	Temp
All	W	4.55	0.001	Cond, pH, Slope
All	F	2.13	0.01	Order
Strong	E	3.83	0.001	OHCov, Elev, SI
Strong	K	2.64	0.005	Cond, Chla
Strong	A	3.20	0.037	Temp
Weak	U	3.32	0.001	Temp, pH
Weak	T	2.57	0.001	OHCov, Pfankuch_bottom, Chla, Depth
Weak	K	2.20	0.024	Cond
Weak	F	2.13	0.018	Order

DOI: 10.7717/peerj.4898/table-2

Note:

Only if a global model was significant, was forward selection performed. Forward-selected variables are given in the ‘Variables’ column. Subset, subset of species (all species, and strong and weak dispersers). Full results of both global and forward-selected models, including spatial variables can be found in Table S1.

Discussion

As a result of the relatively high latitude of New Zealand and based on the hypotheses of Jocque et al. (2010), we hypothesised (H₁) a dominant role of species sorting and dispersal surplus (reflecting the mass effects archetype) in structuring these assemblages (H_1a) and an increasing dispersal surplus from north to south (H_1b). However, despite a latitudinal gradient present in assemblages at the community level overall and within each island for regional γ diversity (as well as α and β diversity, Astorga et al., 2014), what emerged at the metacommunity level was more idiosyncratic. In particular, there was no latitudinal trend in either environmental vs. spatial control (rejecting H_1b) or the idealised metacommunity types tested through the EMS analysis at both the full community level and for dispersal groups. Lack of fit to the hypothesis of Jocque et al. (2010) likely reflects the dynamic, unpredictable nature of New Zealand streams (partially supporting H₁A).

New Zealand comprises a series of mid-latitude islands, with a typically unpredictable climate (Fig. 1) and flashy river flow regimes (Winterbourn, Rounick & Cowie, 1981) reflecting its oceanic position. At a single time-point, communities are therefore most likely at different stages of post-flood recolonisation (H₁A). Antecedent conditions are not only important for dynamic systems like these, but also for more continental climates. For instance, preceding-year climatic conditions have been demonstrated to be more important in shaping European stream invertebrate communities than long-term climatic trends (Jourdan et al., 2018). The dynamism of streams, particularly in oceanic climates, represents a fundamentally important factor controlling metacommunity dynamics, with assembly mechanisms varying temporally in dynamic streams (Datry, Bonada & Heino, 2016; Sarremejane et al., 2017). The relative roles of local and regional processes will depend on the amount of time that has passed for dispersal and colonisation to play out (Brendonck et al., 2014). With the central importance of natural cycles of flooding and drought in streams (Poff et al., 1997; McMullen et al., 2017; Tonkin et al., 2018b), it stands to reason that antecedent flow conditions play a key role in structuring metacommunities in streams (Campbell et al., 2015).

The lack of seasonality and predictability in New Zealand’s climate likely plays a strong role in the low predictability in metacommunity structuring. The hypothesis of Jocque et al. (2010) does not take into account differences in island size and isolation, fundamental aspects controlling biodiversity (MacArthur & Wilson, 1967). Island and mainland locations at similar latitudes do not comprise the same climatic patterns (Tonkin et al., 2017a), with continental locations having much greater predictability in their seasonality compared to islands. To demonstrate this point, we compared a 30-year sequence of monthly rainfall totals from the central North Island of New Zealand with Western Australia, a Mediterranean climate, using wavelet analysis (Fig. 1) (Torrence & Compo, 1998). Although this is just one of the locations examined in our study, which vary in their rainfall regimes, we use this simple comparison to demonstrate the extent of climatic unpredictability present in this region compared to a predictable climatic zone. Figure 1 demonstrates clearly the strongly seasonal and predictable pattern apparent in Western Australia, with a significant and repeatable cycle at the one-year time period over the full sequence. By contrast, central New Zealand’s climate exhibits no repeatability in the rainfall, with very few time points in the sequence indicating any power at the one-year period.

New Zealand streams have other features that may limit their fit to our primary hypotheses, some of which are shared by other island localities, including: rivers tend to be short, swift, and steep due to the narrow landmass and tectonically active nature; evergreen vegetation dominates the flora; and riparian vegetation is scarce for much of their length leading to a predominance of autochthonous rather than allochthonous control of river food webs (Winterbourn, Rounick & Cowie, 1981; Thompson & Townsend, 2000). As such, New Zealand streams are considered as being physically, rather than biologically, dominated systems (Winterbourn, Rounick & Cowie, 1981). These factors, in conjunction with its highly dynamic geological history, have led to the evolution of a stream invertebrate fauna with flexible and poorly synchronised life histories, and generalist feeding behaviour (Winterbourn, Rounick & Cowie, 1981; Thompson & Townsend, 2000; Scarsbrook, 2000). Although New Zealand stream invertebrate communities are not necessarily less species rich or different in terms of food web structure to overseas locations, there is a clear paucity of shredder species in particular, with generalist browsers predominating communities (Thompson & Townsend, 2000). Under these circumstances, it is not surprising that metacommunity dynamics can be difficult to predict, as we clearly demonstrate, without a strong temporal resolution in the data. Thus, in support of our alternative first hypothesis, despite the large latitudinal gradient examined, predictable metacommunity dynamics appear to be masked by short-term unpredictability in environmental conditions.

Results were highly idiosyncratic between different regions, with considerable variability in the relative roles of environmental and spatial structuring, important environmental variables, and the idealised metacommunity types, with no real match between the two approaches. This context dependence did not reflect an interaction between spatial extent and dispersal ability. Although much of this unpredictability may be related to the unpredictable characteristics of New Zealand streams, it is pertinent to recognise that this is a challenge facing many stream metacommunity studies globally, where patterns differ considerably between different catchments (Heino et al., 2012, 2015a; Tonkin et al., 2016a). Lawton (1999) pinpointed this problem in ecology over a decade ago suggesting that community ecology is rife with contingency, so much so that generality is unlikely. Lawton goes on to highlight that the problem is indeed most severe at the intermediate organisational level of communities, compared to more predictable lower (e.g. populations) or higher levels (e.g. macroecology). Metacommunities are indeed difficult systems to predict, with processes affecting different subsets of organisms and operating at specific times (Driscoll & Lindenmayer, 2009). One source of context dependence in metacommunity structuring is differences between different trait modalities, such as dispersal modes (Thompson & Townsend, 2006; Canedo-Arguelles et al., 2015; Tonkin et al., 2016b). Thus, if spatial extent and dispersal limitation were interacting to structure the metacommunity, deconstructing the full assemblage into dispersal groups (e.g. strong vs. weak dispersers) should have helped to explain discrepancies in our predictions, but this was not the case. Nevertheless, we must also entertain the possibility that greater spatial replication would have strengthened the observed patterns.

Finally, contrary to the expectation of Jocque et al. (2010) that dispersal ability increases moving away from the equator (H₂), we found a decrease in the ratio of strong to weak dispersers moving from north to south. Theoretically, temporal variability in environmental conditions promotes increased dispersal ability of organisms (Dynesius & Jansson, 2000; Jocque et al., 2010); an hypothesis strongly tied with Rapoport’s rule of increasing range size with increasing latitude (Stevens, 1989) and one that receives support from the population genetics literature via increased genetic divergence among populations nearer the equator (Eo, Wares & Carroll, 2008). However, it is important to note that while dispersal ability can play a strong role in determining species range sizes, its influence may be less common than previously thought (Lester et al., 2007). Although there is evidence that weak dispersers have stronger latitudinal diversity gradients than strong dispersers in Europe, the mechanisms behind this are related to the ability of organisms to recolonise northern sites following glaciation (Baselga et al., 2012); a different issue to that experienced in New Zealand. The conflicting result we observed may reflect several factors. (1) Lack of time for dispersal and colonisation to play out post-disturbance (Brendonck et al., 2014; Campbell et al., 2015). (2) The requirement of a longer latitudinal gradient for these mechanisms to play out. Over the length of New Zealand, the continuity of habitat availability in space and time, a key mechanism behind Jocque et al. (2010), likely differs very little. (3) Climatic idiosyncrasies not reflecting a north–south gradient and thus not selecting for a gradually increased dispersal ability at higher latitudes.

Conclusions

Jocque et al. (2010) highlighted the fundamental role of dispersal in driving the latitudinal diversity gradient, suggesting a climate-mediated dispersal–ecological specialisation trade-off as a key factor regulating this pattern. We tested several hypotheses based on those of Jocque et al. (2010) relating to how New Zealand stream invertebrate metacommunity structure changed along a broad latitudinal gradient, and examining the mediating role of dispersal. We rejected our primary hypotheses, finding that: (1) species sorting appears to be weak or inconsistent, and its influence did not change predictably with latitude; and (2) weaker dispersers increased with latitude. We associate this lack of fit to these hypotheses on the strong unpredictability of New Zealand’s dynamic stream ecosystems (supporting H₁A) and a fauna that has evolved to cope with these conditions. While local community structure turned over along this latitudinal gradient, metacommunity structure was highly context dependent and dispersal traits did not elucidate patterns.

These results, along with other recent findings (Heino et al., 2012, 2015a; Tonkin et al., 2016a), provide a cautionary tale for examining singular metacommunities. The considerable level of unexplained context dependency suggests that any inferences drawn from one-off snapshot sampling may be misleading. Given the importance of understanding metacommunity processes for the successful management of river ecosystems (Siqueira et al., 2012; Heino, 2013; Tonkin et al., 2014; Stoll et al., 2016; Swan & Brown, 2017), this level of unpredictability is a major cause for concern. While spatial replication of multiple metacommunities may elucidate some of this uncertainty, studies on temporal dynamics of metacommunity processes are clearly needed. We therefore urge researchers to consider the temporal dynamic, particularly in relation to seasonal cycles and their predictability.

Supplemental Information