Consistent sorting but contrasting transition zones in plant communities along bioclimatic gradients

Correlations between plant species occurrences and climate are used as evidence of ecological sorting and climate sensitivity. However, because observed patterns may be spatially and historically contingent, interpretations of compositional responses to spatial differences in climate should ideally consider past climatic fluctuations, edaphic factors, nonlinear turnover rates and transition zones (ecotones) in environmental and geographic space. We tested whether drivers of plant community compositional turnover and transition zones are consistent, rather than spatially contingent, across spatially isolated and biologically contrasting (independent) ecosystems. We used four broad-scale transects in Australia, which are spatially isolated and climatically and ecologically diverse: NATT (Northern Territory); BATS (New South Wales); SWATT (Western Australia) and TREND (South Australia). We used Generalised Dissimilarity Modelling to ascertain the relative contributions of space (geographic distance) and environmental variables to compositional turnover. We mapped transition zones using modelled nonlinear responses. Purely spatial contributions were estimated using variance partitioning. Models explained 26–52% of deviance in compositional dissimilarity. Mean annual rainfall and geographic distance were the most important explanatory variables, although topographic relief outranked rainfall for the topographically complex BATS transect. Purely spatial structure accounted for only 7–12% of explained deviance, pointing to the shared prevalence of species sorting along spatially structured rainfall gradients. Zones of rapid compositional turnover across bioclimatic gradients were evident across the contrasting ecosystems. These transitions zones tended to occur at the arid extremes of rainfall gradients and as plains transitioned into complex topography. Responses to other variables were inconsistent among transects. Transition zones were spatially contingent, influenced by local topography, and generally occurred within different subintervals (value ranges) along gradients. The somewhat idiosyncratic compositional responses suggest that assessment of broadscale compositional drivers including climate change should incorporate regional variation. However, the consistently strong species sorting along rainfall gradients demonstrates a common key response to moisture


Introduction
Correlations between plant species occurrences and climatic factors are frequently tied to ecological sorting and used to infer climate sensitivity. However, using a single system to model compositional turnover in plant communities has limitations because patterns observed locally along spatial environmental gradients can be confounded by historical and biogeographical factors, such as past disturbance, dispersal limitation, migration barriers, edaphic features, ecological drift and imprints of past climatic regimes (Dexter et al., 2012;Warren et al., 2014). These factors can lead to spatial structure, and boundaries in species composition, that can add to, or be confused for, those caused by current climate (Nekola and White, 1999;Jácome et al., 2007;Dexter et al., 2012).
The relative importance of spatial and environmental predictors of compositional turnover is related to dominant assembly processes (Qian et al., 2017). For example, if niche-based processes have dominated the sorting of species along gradients, we would expect to find that environmental variables are better predictors of composition than purely spatial variables. Additionally, plant communities differ in whether species composition responds principally to variation in temperature or rainfall , a critical factor for understanding spatial patterning and predicting future change. Consequently, compositional turnover patterns differ from region to region due to unique histories and spatial structuring of the environment (Peres-Neto et al., 2012).
In addition to differences in the relative importance of space, temperature and rainfall on composition, turnover rates are often nonstationary along environmental gradients (Oksanen and Tonteri, 1995;Ferrier et al., 2007;Blois et al., 2013). Ecotones delineating different ecosystems can result from either rapid turnover zones along environmental gradients or spatial and historical processes, such as dispersal limitation (Risser, 1995;Dexter et al., 2012;Caracciolo et al., 2014). Ecotones may represent thresholds that are likely to be particularly sensitive to climate change (Risser, 1995;Fitzpatrick et al., 2013;Caracciolo et al., 2014;Andersen et al. 2015aAndersen et al. , 2015bBaselga et al., 2015).
The identification of macroecological transitions zones (regional ecotones characterised by particularly rapid compositional turnover) along bioclimatic gradients is therefore of value for predicting where compositional change is likely to be greatest under a future climate (Risser, 1995;Fitzpatrick et al., 2013;Jones et al., 2016), and identifying foci for surveillance monitoring (Guerin et al., 2013). Increasingly realistic interpretations of such compositional patterns have been made possible by correlative approaches that recognise that turnover rates are rarely static along environmental gradients, at different scales or through space (Ferrier et al., 2007;Fitzpatrick et al., 2013;Guerin et al., 2013;Warren et al., 2014), and which can disentangle spatial from environmental effects (Goslee and Urban, 2007).
Given that the interpretation of climate-vegetation patterns is potentially confounded by local (e.g. coincidence with topography) and historical (e.g. dispersal from refugia) factors, it is pertinent to ask whether drivers of compositional turnover and associated transition zones are consistent among spatially independent ecosystems. Empirical studies show that compositional turnover in regions with different histories and spatial contingencies can have different drivers and rates (Peres-Neto et al., 2012). For example, Fitzpatrick et al. (2013) compared the drivers of compositional turnover in south-west Western Australia and northern Europe within the context of the Last Glacial Maximum (LGM), during which northern Europe experienced glaciation that forced plant into refugia (Hewitt, 2004), while Australia remained unglaciated but underwent climatic fluctuations (Byrne, 2008). Fitzpatrick et al. (2013) found that spatial structure and temperature were leading drivers of plant species composition in Europe, whereas rainfall was more important in Australia, suggesting that glaciation and subsequent recolonisation caused higher spatial structuring in Europe due to dispersal lag.
If compositional turnover patterns reflect species sorting along environmental gradients over dispersal limitation (Nekola, 1999), we would expect to find little evidence of purely spatial structure in compositional patterns. That is, we would expect the environment (including its covariance with space due to geographically arranged gradients) to be more important than space per se in the absence of large-scale disturbance. However, fire events are important determinants of spatial structure in vegetation type and composition in susceptible Australian ecosystems (Baker et al., 1991;Hill, 2017).
Here, we investigated how compositional turnover in diverse Australian plant communities is influenced by spatial and environmental factors, to test the extent of shared patterns. We used broadscale transects that are spatially isolated and situated within different Australian eco-regions and climate types, including semi-arid, tropical and Mediterranean. Repeated studies across bioclimatic gradients in statistically independent, or at least spatially separated, systems (as compared to local replication for statistical validity; Andersen et al., 2015b) have the ability to detect general patterns above the historically contingent idiosyncrasies of individual regions and ecosystems (Jácome et al., 2007;Fitzpatrick et al., 2013;Borer et al., 2014;Tuomisto et al., 2016;Caddy-Retalic et al., 2017).
We addressed the following research questions: 1) How does the relative importance of spatial structure, rainfall, temperature, edaphic and topographic factors vary along independent bioclimatic transects? 2) To what extent are zones of more rapid compositional turnover consistent among independent bioclimatic gradients?

Study systems
The study analysed four independent bioclimatic transects on mainland Australia (Fig. 1) from the plot network of the Terrestrial Ecosystem Research Network (TERN): BATS (Biodiversity and Adaptation Transect Sydney), NATT (Northern Australian Tropical Transect), SWATT (South West Australian Transitional Transect) and TREND (Transects for Environmental Monitoring and Decision-making) (Hutley et al., 2011;Guerin et al., 2014b;Caddy-Retalic et al., 2017;Gibson et al., 2017). The transects span a range of eco-climate types across the continent described by Hutchinson et al. (2005), as follows: seasonal tropics to desert (NATT); warm-wet to cool-wet temperate (BATS); and Mediterranean to semi-arid (SWATT and TREND). Although the floras of these regions overlap at genus and family levels, they are highly distinct at species level and differ in the relative proportion of plant functional types. For example, 79% of plant species in the SWATT region are endemic to Western Australia (Beard et al., 2000), and very few species on NATT occur on other transects. C4 grasses are more important in the north than in the south relative to C3 grasses (Johnson et al., 1999). These differences make it possible to test for generality, Fig. 1. Location of the four independent bioclimatic transects in Australia and their spatial definition for this study. The transects sample a wide range of climatic zones (cool-wet, desert, Mediterranean, seasonal tropics, semi-arid and warm-wet) but overlap in mean annual precipitation ranges, and their floras are largely independent at species level, making it possible to test for generality among ecosystem transitions along major temperature and rainfall gradients. independent of response types of particular species or the evolutionary or physiological characteristics of particular assemblages. This 'networked transect' approach also helps to counter the potentially confounding effects of topography and local colinear geographic gradients (Caddy-Retalic et al., 2017), partly because the transects vary substantially in topographic complexity (Appendix S2).
For this analysis, the NATT was defined by a ∼200 km wide quadrilateral with vertices at (134.5°E, 17.83°S), (133.5°E, 12.33°S), (131.5°E, 12.33°S) and (132.5°E, 17.83°S). BATS was defined by a ∼70 km wide quadrilateral with vertices at (150.1178°E, 32.83988°S), (151.4172°E, 33.49183°S), (151.1456°E, 34.08963°S) and (149.7930°E, 33.39708°S). However, the Cumberland Plain region within BATS was deliberately excluded from sampling in order to restrict the analysis to sandstone vegetation to minimise the confounding effects of the very different vegetation types known to occur on these landscapes. SWATT was defined by a ∼200 km wide quadrilateral with vertices at (116.

Composition datasets
For NATT and SWATT, analyses were based on herbarium records obtained via Australia's Virtual Herbarium (https://avh.chah.org.au/, accessed 17/08/2016). Records were filtered for spatial and taxonomic validity, and to records with coordinates accurate to within the spatial resolution used for the analysis. The datasets consisted of 79,306 records of 2,806 species for NATT; and 50,744 records of 4,068 species for SWATT. For BATS and TREND, species inventory data from vegetation surveys were available and were used in preference to herbarium records as they are expected to be less biased (Channon and Heard, 1997;Franklin et al., 2017;Guerin et al., 2018). For BATS there were 114,736 records of 1,308 species from 3,468 field survey sites, and for TREND there were 125,085 records of 2,145 species from 4,740 sites. Species occurrences were compiled into matrices representing map grid cells (Guerin et al., 2015) for each transect at a resolution of 0.2°for NATT and SWATT, and 0.1°for BATS and TREND, which had denser species data and more rapid turnover of climatic variables.

Model fitting and comparison
We used Generalised Dissimilarity Modelling (GDM) to assess compositional turnover with respect to environmental gradients and spatial structure. In this context, 'compositional turnover' is defined as the pairwise dissimilarities among sample units (in this case map grid cells) as determined by the Bray-Curtis index (Ferrier et al., 2007). GDM accounts for the curvilinear relationship of compositional dissimilarity with environmental distance and allows the rate of turnover to change at different points along gradients, by fitting a series of I-spline functions (Ferrier et al., 2007). In a previous analysis, GDMs were shown to reproduce similar patterns of compositional turnover between datasets based on herbarium records versus plot-based inventory data . However, it was also shown that plot-based inventory data were more robust to sampling bias and should be used in preference to herbarium data wherever possible . Nevertheless, the outputs of GDM need to be interpreted appropriately (Prober et al., 2012;Guerin et al., 2018;See Discussion).
A set of environmental predictor variables was compiled that included climate data from WorldClim (raster data at 30 arc-second or approximately 1 km resolution; Hijmans et al., 2005, Table 1) and edaphic/topographic variables from the Soil and Landscape Grid of Australia (raster data at 3 arc-second or approximately 90 m resolution; Grundy et al., 2015), with values extracted at the centroids of grid cells. We included candidate variables to represent major gradients typically expected to influence plant distribution. Pairwise and multiple colinearity among the predictor variables were assessed for each transect (Dormann et al., 2013). Variables were excluded to ensure that no pairs initially had a Pearson's r greater than 0.8, and, subsequently, that variance inflation factors (VIFs) were below 10. Where a selection among variables was necessary to reduce highly correlated pairs or high VIFs, those with weak importance in exploratory models were excluded, with a preference to retain at least one temperature and one rainfall variable.
We allowed five splines in the non-linear responses and included geographic distance as a covariate to account for spatial autocorrelation and to disentangle the influence of gradients from spatial processes (Goslee and Urban, 2007). Observations were weighted according to species richness to minimise the expected influence of under-sampling on observed richness within cells . A different subset of predictor variables was included in the final model for each transect following a process of backward elimination. Variables were sequentially removed when model coefficients summed to zero, indicating no effect on species turnover. Among variables that had at least one non-zero coefficient, we subsequently removed those found to be non-statistically significant using the method of Leitão et al. (2017). For final models, deviance explained was calculated and the variables were ranked by importance according to the height of their modelled functions on the y-axis, which represents the total amount of compositional turnover along that gradient and is equivalent to the sum of the I-spline basis function coefficients (Ferrier et al., 2007).
Modelled responses can be influenced by the inclusion of different co-variables, which limits comparability. We therefore repeated GDMs using a common subset of predictors with no further vetting procedure by including each predictor that was included in at least two final models. These 'standard' models for each transect were otherwise dealt with as for the above 'exploratory' models.
In addition, we used variance partitioning to break the explained deviance in compositional turnover into three components: purely spatial structure; purely environmental; and spatial-environmental covariation (Borcard et al., 1992). The relative importance of these components can be related to the importance of dispersal versus environmental filtering in explaining compositional turnover (Gilbert and Bennett, 2010). We used deviance explained in additional GDMs that included only spatial, only environmental, or both spatial and  Gilbert and Bennett (2010).

Identification of transition zones
For the purposes of this analysis, we defined transition zones in two ways and subsequently mapped them for each transect with respect to the most important rainfall and temperature variables, using the fitted 'exploratory' GDMs above. The definitions highlight variation within transects and also provide a basis for cross-transect comparisons. Note that defining transition zones by set cut-off points in compositional turnover is impractical for comparison since absolute values may differ between models and transects for reasons outlined in the Discussion.
1. Continuous spatial method: Using GDM-transformed predictor layers, in which cell values represent height on the y-axis of GDM functions, we calculated mean differences among grid cells in 11 × 11 cell moving neighbourhood windows. This process maps turnover as a combination of turnover rate per unit of the predictor and rate of spatial change in the predictor. 2. Pre-defined gradient breadth: For the second definition, we identified subintervals of fixed breadth (but flexible position) along the gradients over which the turnover rate per unit of the predictor was highest. This method disregards the rate of spatial change in the predictor and therefore also excludes any transition zones resulting purely from spatial discontinuities such as abrupt altitudinal changes, which are highlighted in method 1. We generated a 200long numeric sequence representing the full range of the predictor in the training data and calculated associated GDM-predicted values from fitted I-splines. Next, we calculated slope across sequential, overlapping subintervals that were one eighth of the total range, resulting in 176 unique subintervals. We then selected the subinterval with the highest slope.

Model fitting and comparison
For 'exploratory' GDMs, the deviance explained was 26% for NATT, 41% for BATS, 44% for SWATT and 49% for TREND (Table 2). Between three and nine environmental variables were retained ( Fig. 2; Table 2). Geographic distance (Geo) and mean annual rainfall (Bio12) were consistently the top ranked predictors, except for BATS, where topographic relief (Relief) out-ranked mean annual rainfall. The inclusion of variables in GDMs and their relative importance differed among transects as detailed in Table 2.
Geographic distance and mean annual rainfall were retained as predictors in each 'exploratory' model, while mean maximum temperature of the warmest month (Bio5), mean temperature (Bio1), topographic relief, soil effective cation exchange capacity (ECE) and soil available water capacity (AWC) occurred in two models each. These variables were therefore selected and combined as predictors for 'standard' models'. For 'standard' GDMs, deviance explained was 41% for BATS, 43% for NATT, 50% for TREND and 52% for SWATT. Geographic distance and mean annual rainfall were the top predictors for SWATT and TREND but were out-ranked by mean temperature for NATT, while mean annual rainfall was out-ranked by topographic relief for BATS. The ranking of the remaining variables differed considerably Table 2 Statistics for Generalised Dissimilarity Models (GDM) fitted to plant compositional turnover for four independent bioclimatic transects across Australia (Fig. 1)   among the transects, as detailed in Table 2. Purely spatial structure (i.e. that independent of spatial-environmental covariation) accounted for a relatively small proportion (0.07-0.12) of the deviance explained by 'standard' models (Fig. 3). For NATT and SWATT, spatial-environmental covariation comprised the biggest component of deviance explained, leaving a proportion of 0.26 for the purely environmental component for NATT and just 0.10 for SWATT. In contrast, for BATS and TREND, purely environmental variation and spatial-environmental covariation components were nearly equal (Fig. 3).
Non-linear responses to individual predictors were modelled in 'exploratory' and 'standard' GDMs with a range of sigmoidal, nearlinear or plateauing functions ( Fig. 4; 5). Responses to mean annual rainfall followed a common pattern in which compositional turnover per mm of rainfall was rapid at the more arid extremes of the gradient, with the rate slowing over increasingly wet landscapes and typically plateauing in higher rainfall zones ( Fig. 4c; Fig 5c). The plateauing was more obvious for BATS and NATT, which extend into rainfall zones of over 1200 mm. The two Mediterranean-semi-arid transects, SWATT and TREND, covered similar mean annual rainfall ranges, from approximately 200 to 1100 mm, and had overlapping fitted I-spline basis functions with similar shapes and coefficients (Fig. 5c). Similar responses were also evident to topographic relief, in which there was a consistent transition zone from plains (i.e. zero relief) to more topographically complex landscapes but the response plateauing so that increasing topographical complexity was not associated with additional compositional turnover (Fig. 5d).

Identification of transition zones
Transition zones along gradients were explored in more depth for mean annual rainfall as well as the most important temperature  (Table 1 for codes) predicting plant compositional turnover along the four bioclimatic transects in Australia shown in Fig. 1 (BATS, NATT, SWATT, TREND). Y-axes represent the sum of the I-spline coefficients in fitted Generalised Dissimilarity Models (GDMs), which is equivalent to the height of the partial functions and the amount of species turnover along that gradient. 'Exploratory' refers to model selection via a backwards elimination procedure, while 'Standard' refers to models fitted to a pre-defined subset of variables for comparison. variables in 'exploratory' models: mean temperature for BATS and SWATT; mean maximum temperature of the warmest month for NATT and TREND.
For mean annual rainfall gradients (Fig. 6), spatial transitions (definition 1.) reflected topographic discontinuities for BATS and TREND. For BATS, there was a clear transition associated with mountainous terrain over 1,000 m asl (Appendix S2), whereas for TREND, the transition zone was at the boundary between ranges and lowland plains. For SWATT, the northern arid extreme of the transect was highlighted as a transition zone and, to lesser extent, the southern coastal margin. For NATT, two distinct bands were highlighted as transition zones, at the southern arid extreme of the transect and at the mid-point. However, the rate of change was low for NATT compared to the other transects, while the most rapid rainfall-related transition occurred on TREND.
Mean annual rainfall transitions (subintervals with the highest compositional turnover rate; definition 2.) overlapped for SWATT and TREND, within ranges of 206-311 mm and 238-340 mm, respectively, for 'standard' models (Table 2; Figs. 7 and 8). The rainfall transitional zones for BATS (878-975 mm) and NATT (510-628 mm) differed. Note, however, that these transects occupy different intervals along the rainfall gradient. When mapped spatially, rainfall transition zones mapped for BATS to the margins of the Cumberland Plain and higher elevation areas at the western side of the Dividing Range; for NATT at the arid southern end of the transect; for SWATT throughout the northeastern arid half of the transect; and for TREND along an irregular central strip matching the 'mallee' belt of the Flinders Ranges and Murray Mallee (Fig. 7).
With the exception of NATT, the rate of compositional turnover across temperature transitions was lower than along mean annual rainfall transitions. Transitions coincident with temperature were evident for NATT in the north-eastern corner of the transect, associated with sandstone and higher elevation terrain; for SWATT, along the southern coastal margin and, to a lesser extent, the mid-transect regions of the transect; for BATS along a band through the western end of the transect; and for TREND along the margins of the ranges (Fig. S9 of Appendix S1). The rate of the transition was relatively low for TREND compared to NATT and SWATT.
Temperature transition zones were inconsistent among the transects, although comparisons are limited by the fact that each transect occupies different subintervals along the gradient. For mean maximum temperature of the warmest month, there was some overlap in transition zones, with several models predicting them at 34-36°C or 27-29°C. For temperature, transition zones mapped for BATS to the transition to higher ranges; for NATT the Arnhem Plateau in the far north-eastern corner of the transect; for SWATT to a strip across the arid half of the transect; and for TREND to the northernmost plains (Fig.  S10).

How does the relative importance of spatial structure, rainfall, temperature, edaphic and topographic factors vary along independent bioclimatic transects?
Previous Australian studies have found that rainfall is a key driver of vegetation composition and structure (Williams et al., 1996;Hutley et al., 2011;Fitzpatrick et al., 2013;Guerin et al., 2017). When tested explicitly, the contribution from purely spatial structuring has been found to be relatively small . Our results largely match this conclusion, providing further support to the hypothesis that strong sorting along rainfall gradients has dominated geologically recent plant community assembly throughout Australia. Deterministic environmental controls on present-day species composition are expected to translate to higher levels of predictability under climate change (Ozinga et al., 2005;Van Bodegom et al., 2014), although the inference is still indirect.
The four independent bioclimatic transects we studied in Australia shared the major attributes of strong sorting along mean annual rainfall gradients and low purely spatial structuring. However, specific responses to other variables, such as temperature and edaphic or topographic factors, differed markedly. The unique contribution of space was consistently weak across all transects, as ascertained with variance partitioning (Fig. 3). The higher contribution of pure environment and spatial-environmental covariance points to a high degree of niche-based species sorting along spatially structured gradients. The higher proportion of deviance explained by spatialenvironmental covariance for NATT and SWATT may reflect the relative lack of topographic discontinuities in climate compared to BATS and TREND (Appendix S2), where purely environmental effects were stronger. That is, environmental gradients on less topographically complex transects tend to shift more evenly with increasing geographic distance (Hutley et al., 2011).
Mean annual rainfall was identified as an important environmental predictor of compositional turnover across the transects. However, mean annual rainfall was outranked by topographic relief in two models for the BATS transect, and by mean temperature in one NATT model (Fig. 2). Mean temperature also outranked both rainfall and geographic distance for the 'standard' BATS model (the variable was not statistically significant in the 'exploratory' model). The response to mean annual rainfall was particularly similar for the SWATT and TREND transects, highlighting a common biogeographic transition from Mediterranean to arid climates and floras (Fig. 5). Changes in mean annual rainfall had a more pronounced effect across lower rainfall transition zones, underlining the importance of moisture availability in low rainfall systems (February et al., 2007, Fig. 8).
Landscape factors, such as topographic relief and soil nutrient status, were typically less important for predicting compositional dissimilarity than climate variables and geographic distance. The exception was the finding that topographic relief was more important than Fig. 3. Partitioning of the proportional contribution of purely spatial structure versus environmental variation to deviance explained by Generalised Dissimilarity Models (GDMs) fitted for the four transects in Australia using a 'standard' set of climate, soil and landscape predictors (see Table 2. Fig. 5 for further details). rainfall for BATS, the transect with the highest altitudinal range (Appendix S2). The deliberate orientation of the transects to sample climate, rather than soil, gradients may have influenced this finding. An additional cause of the low importance of soil variables in our models may be the coarse resolution of the analysis of grid cells, in which small-scale edaphic variation is not reflected but that may have particular importance for explaining vegetation patterning at local (Staver et al., 2011) or global (Bruelheide et al., 2018) scales.  Fig. 1 (BATS, NATT, SWATT, TREND). Plots show nonlinear, monotonic functions fitted to variables using Generalised Dissimilarity Models (GDM), in which slope represents the instantaneous rate of compositional turnover and maximum function height represents total compositional turnover and is equivalent to summed spline coefficients. Selected variables are shown based on importance in multiple transect models. Statistics are shown in Table 2

To what extent are zones of especially rapid compositional turnover consistent among independent bioclimatic gradients?
The detection and characterisation of ecotones are useful considerations for climate adaptation planning, because transition zones of more rapid compositional turnover may be more sensitive to climate change (Peñuelas & Boada, 2003;Guerin et al., 2013). The transition zones we detected in geographic and environmental space resulted from nonlinear rates of compositional turnover and environmental gradients that are not spatially uniform. Generally, the transition zones were somewhat spatially contingent rather than fitting a fixed pattern. For example, transition zones were associated with topographic  Fig. 4 for details). 'Standard' models were generated using the same set of pre-defined variables for each transect for comparison. (a) Mean temperature (bio1); (b) Mean maximum temperature of the warmest month (bio5); (c) Mean annual rainfall (bio12); (d) Topographic relief; (e) Soil effective cation exchange capacity (ECE); (f) Soil percent available water capacity (AWC). discontinuities from mountain ranges to adjacent plains (TREND) or onto high country (BATS), but for NATT and SWATT the transitions were less obviously associated with geographic features.
Compositional transition zones were found in the lower rainfall regions of each transect, pointing to the most consistent modelled pattern among transects. Indeed, the rainfall transition zones for TREND and SWATT, the Mediterranean-semi-arid transects, were almost identical in environmental space (Table 2; Fig. 8). Nevertheless, the spatial arrangement of these compositional transitions on the ground is also strongly influenced by environmental transitions. For example, complex topography creates discontinuities in rainfall patterns for BATS and TREND, leading to rapid changes in species composition over small geographic distances ( Fig. 6; Appendix S2).

Implications for climate sensitivity assessment practice
Several studies have assessed the drivers of species distribution and composition using models that encompass large, even continental, areas to infer relative climate change sensitivity (McMahon et al., 2011;Ferrier et al., 2012;Prober et al., 2012). Broadscale approaches provide context to local predictions, enable regional prioritisation and capture the broad climatic tolerances of species and communities (Pearson and Dawson, 2003). However, we found that, despite some important consistencies, compositional turnover responses differed among the four studied transects in different regions of Australia. In particular, identified transition zones of rapid compositional turnover were, in many cases, localised and occupied contrasting regions of climatic space and topography. Given these regional differences, assessments of climate responses or climate change vulnerability across large spatial extents should ideally consider nested spatial extents to account for contrasting responses within and between ecosystems, in addition to continental patterns (Randin et al., 2009;Siegel et al., 2014).

Limitations
There are important limitations to the interpretation of compositional turnover along bioclimatic transects. Compositional dissimilarity models such as GDM represent mainly deterministic variation along a particular set of tested gradients. In doing so, they ignore other sources of variation that may not be fully captured by the decay of compositional similarity with geographic distance, such as disturbance regimes and stochasticity (Guerin et al., 2014a;Baselga et al., 2015). Additionally, inherent biases in sampling regimes for biological data may influence our ability to detect real changes, although GDM appears to be relatively robust in this regard compared to raw data approaches .
The fit and predictions of various species composition models vary partly according to how they are trained (Jones et al., 2008). In our study, a small, standardised selection of spatial and environmental variables had good explanatory power when applied across diverse ecosystems, indeed often better than variables selected for individual regions through a statistical vetting procedure. This counterintuitive result is partly due to the exclusion of variables in exploratory models during significance testing that nonetheless contributed substantially to deviance explained in standard models. Responses to particular variables, including their coefficients and rank importance, were also somewhat influenced by co-variable selection (Zuur et al., 2010;Grueber et al., 2011), notably between exploratory and standard models.
Differences in models and subsequent spatial-environmental variance partitioning can emerge when different datasets, spatial extents, input variables and selection procedures are used for the same general region or ecosystem (Gilbert and Bennett, 2010;Guerin et al., 2018). The phenomenon is demonstrated by differences in variable importance for south-west Western Australia reported here compared to Fitzpatrick et al. (2013) and Jones et al. (2016). The importance of maximising variance explained versus robust detection of real drivers among candidate variables will vary according to the purpose of each study  . 8. Transition zones -environmental space: subintervals of mean annual rainfall (bio12) gradients along which Generalised Dissimilarity Models (GDM) predicted the highest rate of compositional turnover per mm of rainfall for each transect. (Murray and Conner, 2009). Nevertheless, these issues suggest it may be most prudent to consider multiple models or take a model averaging approach (Whittingham et al., 2006) when assessing drivers of compositional turnover or climate sensitivity, or at least to interpret results cautiously.

Conclusions
Our analysis of compositional turnover along independent bioclimatic transects aimed to understand the relative importance of purely spatial structure as compared to rainfall, temperature, edaphic and topographic factors. Mean annual rainfall was found to be a key driver of plant compositional turnover across four diverse Australian regions, with caveats. Purely spatial distance decay effects that may relate to dispersal limitation or ecological drift, accounted for only 7-12% of explained deviance. This suggests that environmental filtering or nichebased processes have had more influence than geographic distance per se on recent community assembly and that moisture limitation is a key consideration across Australia.
The compositional dissimilarity responses to climatic gradients were not uniform within regions, and localised transition zones were identified across bioclimatic gradients and appear to be common, regardless of eco-climate type. However, the nature of transition zones was somewhat spatially contingent, in that they differ from region to region according to local factors such as topography. Transition zones across the more arid ends of rainfall gradients were the most consistent finding, although these only overlapped in climatic space for two transects (SWATT and TREND). The differences in compositional turnover along independent sampling transects (e.g. in terms of rank importance of drivers and nature of transition zones) suggest that climate sensitivity assessment should ideally aim to capture regional variation in ecological responses.