Determining the distributions of plant communities in subantarctic vegetation using species distribution models

Keywords: feldmark, grassland, herbaceous plants, herbfield, megaherb, plant communities, random forest model, tundra, vegetation boundaries.

The subantarctic islands have a distinct environment, reflecting their extreme isolation and hyper-oceanic climate. They are characterised by tundra vegetation (short woody vegetation being absent in most cases) with low vascular plant species diversity and a complete lack of native mammalian herbivores (Smith and Lewis Smith 1987; Bergstrom and Chown 1999; Van der Putten et al. 2012). Invasive species have had severe impacts in many cases, for example, rabbits caused massive changes to the vegetation of Macquarie Island prior to their eradication in 2011 (Whinam et al. 2014). Aside from localised long-term monitoring sites (Scott and Kirkpatrick 2013; Whinam et al. 2014), there are limited means of examining how the vegetation of the entire island has responded to herbivore removal.

The vegetation map for Macquarie Island (see Selkirk and Adamson 1998 and subsequently Department of Primary Industries Parks Water and Environment 2013; Supplementary Fig. S1) was produced by manual classification of vegetation height and groundcover, largely from a 1994 satellite image. This method did not consider floristics and therefore the structural classes often include a variety of floristic assemblages (Selkirk and Adamson 1998). Modelling the distribution of individual key plant species could help understand the co-occurrence patterns of species and consequently provide a method for high resolution repeat monitoring of vegetation change at the island scale.

Species niche or distribution models (SDMs) have been widely used in diverse contexts to estimate the potential niches of species (Elith and Leathwick 2009; Franklin 2009). Numerous methodological approaches have been employed to develop species models depending on the purpose of the modelling, the available data and the background of the researcher (Austin 2007; Elith and Leathwick 2009). Species distribution modelling typically determines climatic limits for a species or probabilities of occurrence but does not indicate the optimum habitat or consider interspecific competition (Heikkinen et al. 2007; Anderson 2017). Typical SDMs, therefore, provide limited inferences about plant species assemblages and vegetation types. More complex modelling approaches are required to capture ecological interactions and dynamics (Kissling et al. 2012; Chardon et al. 2020), such as including species interactions as variables in SDMs (le Roux et al. 2013). SDM approaches have also been applied to the mapping of plant species communities or vegetation types (Ferrier and Guisan 2006; Norberg et al. 2019).

SDM projects mostly use presence only or presence or absence data. Following a previous study, which applied a three-class classification of Azorella macquariensis cover in feldmark vegetation on Macquarie Island (Bricher et al. 2013), we applied this approach to modelling a range of key species on the island. This repeatable method has prospective application for monitoring temporal and spatial vegetation change in a more nuanced way than simple presence or absence of species in an environment where vegetation type (e.g. short herbfield, tussock grassland) reflects the local dominance of widespread plant species.

We use a dominance category that is easily determined during field data collection or from existing floristic plot data to represent the core range of a species. We define the core range as the environmental envelope in which a plant species is a dominant species, with a minimum of 25% projected foliage cover. Biomass or leaf area index are more accurate means of determining dominance, but projected foliage cover is easier to measure and is closely correlated to these other metrics in herbaceous tundra vegetation (Chen et al. 2009).

This study aims to compare inferences about vegetation communities from SDMs, by determining the extent of assemblages of co-occurring species, using both the entire range of a species and its core range. We expect that the widespread dominant vascular species will each have similar environmental envelopes, with distributions limited by extremes of climatic and edaphic conditions (total range) but will be more differentiated in their core ranges. We also determine the relative importance of predictor variables for each of the nine species.

Macquarie Island (54.50°S, 158.94°E) emerged from the Southern Ocean less than 1 × 10⁶ years ago and is the only exposed portion of the extensive submarine Macquarie Ridge (Adamson et al. 1996). The 128-km² island is largely composed of igneous rocks such as basalts. Escarpment slopes rise steeply from the ocean on the east coast and from a flat coastal terrace in the west, forming an undulating plateau mostly over 200 m above sea level (ASL), which comprises most of the landmass, with peaks reaching 420 m ASL (Fig. 1).

Fig. 1.  Location of species observation points on Macquarie Island used for species distribution modelling. Inset shows location of Macquarie island (red outline) in the Southern Ocean.

Macquarie Island is subject to constant geostrophic winds and frequent cold fronts, in a region with the highest annually averaged sea level wind speeds on Earth (Hande et al. 2012). The hyper-oceanic climate is characterised by very low diurnal and annual temperature ranges (January mean monthly temperature range, 5.3–8.8°C; July, 1.6–4.9°C), high relative humidity (mean 86%) and almost daily precipitation. Annual mean rainfall (1033 mm for 1980–2010) has increased by ∼20% since 1948 (Jovanovic et al. 2012). Low cloud (fog) is frequent at higher elevations (Fitzgerald and Kirkpatrick 2020).

The vegetation is dominated by bryophytes, ferns, forbs and graminoids. Megaherbs and tussock grasses are the largest lifeforms. Vegetation types recognised on the island include tussock grassland, short grassland, herbfield, fernbrake, mire and feldmark (Selkirk et al. 1990). Feral rabbits (Oryctolagus cuniculus) have had a profound effect on the vegetation of Macquarie Island by intense grazing of select plant species, promoting vegetation dominated by less palatable or more resilient plants (Scott and Kirkpatrick 2013; Whinam et al. 2014). There are signs of vegetation changes following eradication of rabbits in 2011 (Shaw et al. 2011; Williams et al. 2016; Fitzgerald et al. 2021).

Cover data for 9 vascular plant taxa were compiled from 770 plots (Fig. 1) spanning the period 2009–2017 (Bricher et al. 2013; Department of Primary Industries, Parks, Water and Environment, unpubl. data, 2017; Jenny Scott, unpubl. data 2014; Fitzgerald, unpubl. data 2015–2017).

The 9 taxa comprise 10 species (with 2 Acaena species combined into 1 class because they are ecologically similar and sometimes difficult to distinguish) representing most of the dominant vascular plants of Macquarie Island vegetation (Selkirk et al. 1990): Acaena spp., Agrostis magellanica, Azorella macquariensis, Deschampsia cespitosa, Festuca contracta, Luzula crinita, Pleurophyllum hookeri, Poa foliosa, Stilbocarpa polaris (syn. Azorella polaris). Lack of data precluded the inclusion of species largely confined to localised habitats such as mires, coastline and fernbrakes.

A taxon was considered dominant if projected foliage cover exceeded 25%. This threshold was chosen to capture the most extensive species in each plot, including instances where total vegetation cover is as low as 25% or where more than one species has similar coverage (co-dominance). Plots varied in size from 1 × 1 m (in which case several nearby plots were averaged as a single point) to 2.5 × 3.5 m, to a circle of 2.5-m radius. The 410 observations from Bricher et al. (2013) were sampled using a geostratified design to capture environmental gradients on Macquarie Island. Frequency of occurrence of taxa in the plots is shown in Table S1.

We used 11 terrain and topoclimatic predictor layers derived from a 5-m resolution digital elevation model for Macquarie Island (Table 1). Correlations between pairs of predictors did not exceed 0.7 (Pearson correlation coefficient) so we consider each predictor to potentially contribute to the models. The most frequent and strongest winds on Macquarie Island are from the northwesterly and westerly sectors, whereas southwesterly winds appear to be the most limiting to plant growth (Fitzgerald and Kirkpatrick 2017). We therefore modelled local wind speed across Macquarie Island based on the 95th percentile wind speed (from Bureau of Meteorology data) for each of three wind directions (225, 270, 315°) using WindNinja software (Forthofer 2007). The resulting wind exposure layers were highly correlated so we used only the westerly (270°) model for the species modelling. These wind speed models do not incorporate turbulence so may have reduced accuracy at high wind speeds in complex or steep terrain, as is typical on Macquarie Island.

Table 1.  Topographic layers used as predictors for the SDMs, with some examples of how they might influence the distribution of plant species either positively or negatively.

To investigate whether species models could be improved by the addition of spectral reflectance data we also performed the modelling with a set of eight spectral layers. Spectral layers were obtained from a cloud-free 2016 Sentinel-2 image of Macquarie Island. We trialled each of the three visible wavelength bands and the near infra-red band as individual layers plus four spectral vegetation indices derived from these bands using ENVI software (ver. 5.4, Harris Geospatial, Broomfield, CO, USA, see http://www.harrisgeospatial.com/docs/VegetationIndices.html). Test models using spectral variables with topographic variables were not significantly better than topographic variables alone, using the AUC metric (area under the curve of the receiver operating characteristic, a measure of the classification accuracy of the model). Further, modelling with spectral layers was limited to a subset of the species data that was from the same year as the satellite imagery and the spectral layers were highly correlated (>0.7 Pearson correlation coefficient), which can be problematic for modelling purposes (Dormann et al. 2013). Spectral variables were therefore not used in the final models.

Modelling was performed using the biomod2 package for R (ver. 3.5.1, see https://biomodhub.github.io/biomod2/index.html; Thuiller et al. 2009). Exploratory analyses with four types of models (GLM, GAM, Artificial Neural Networks, Random Forest) showed consistently good results from random forest models (RFMs) with our data. Random forest is a machine-learning algorithm for classification or regression (Breiman 2001; Liaw and Wiener 2002). We employed random sampling of the dataset to select 65% of records for training and set aside 35% for evaluation for each decision tree in each RFM. Trial runs of RFMs with a subset of species showed that the default parameters of ntree = 500 (trees per RFM) and mtry = 4 (number of variables randomly selected for each decision tree split) were effective at minimising out-of-bag (OOB) error rates. Unlike some other SDMs, random forest models are not adversely affected by spatial autocorrelation of species observations (Marmion et al. 2009) so we did not investigate spatial autocorrelation in our data. Modelling involved running 100 RFMs and selecting the best model based on the highest AUC score. The species dataset was split into two datasets: presence or absence (PA); and dominant or not-dominant (DN), which comprises all presence records but not absences (i.e. not-dominant is defined as present with <25% cover). For each species two RFM models were produced: PA, which models the total range (presence) of the species; and DN, which models the dominant or core range as a subset of the total range. The resulting RFM regression layers for each species were converted to binary layers using a threshold that maximises the true skill statistic (Hanssen and Kuipers discriminant) (Allouche et al. 2006). These two layers were overlaid to produce a final layer with three classes: absent, non-core range, core range.

To determine which predictor variables were related to the distribution of each species we calculated the mean and standard deviation of the variable importance value from the 100 RFMs and plotted response plots for all variables for each species. Variable importance is calculated in the biomod2 package by randomising the values for a single variable and comparing the model results to the reference model to determine the influence of that variable on the model.

To predict the vegetation from the species projections we employed a ruleset (Table S2) based on observations of plant community assemblages on Macquarie Island and largely following the classification of Selkirk et al. (1990) (Table S3). Essentially the tallest plants with largest leaf area were expected to dominate smaller plants. An assemblage comprising megaherbs and tussock grass could have a single dominant or two co-dominant species. The short grass species and the perennial herb Acaena can occur as three or more co-dominant species. Sites with no modelled dominant species (i.e. all species <25% cover) were classified as feldmark if within the modelled range of A. macquariensis. We applied this ruleset to both modelled core range and total range (i.e. core plus non-core range) models to see how they differ.

Model performance was acceptable with AUC values exceeding 0.78 for the best models in all cases (Fig. 2, Table S4). The mean AUC value for the core range and total range models was not significantly different (two-sample t-test, P = 0.89).

Fig. 2.  Model performance of 100 random forest models was assessed using the area under the curve of a receiver operating characteristic plot (AUC). Bars are mean AUC of 100 models, error bars are one standard deviation. Dark grey bars indicate the presence or absence model, pale grey bars indicate the dominant or not-dominant model. Acaena, Acaena spp.; Agrmag, Agrostis magellanica; Azomac, Azorella macquariensis; Desch, Deschampsia cespitosa; Festuca, Festuca contracta; Luzula, Luzula crinita; Pleuro, Pleurophyllum hookeri; Poafol, Poa foliosa; Stilbo, Stilbocarpa polaris.

The short graminoid species Agrostis magellanica, Festuca contracta and Luzula crinita have ranges covering almost the entire island but their core ranges are limited to low to mid elevations, particularly in the east (Fig. 3). Poa foliosa is mostly confined to the coastal terraces and slopes, whereas S. polaris and P. hookeri have ranges covering a wide range of elevation but mostly in the northern half of the island. The core range of Azorella macquariensis covers the entire high elevation plateau habitat. Deschampsia cespitosa has the smallest and most fragmented range of the species considered here (Fig. 3).

Fig. 3.  Modelled core range (black) and non-core (total) range (pale grey) for nine dominant Macquarie Island plant taxa, classified from random forest models. Acaena, Acaena spp.; Agrmag, Agrostis magellanica; Azomac, Azorella macquariensis; Desch, Deschampsia cespitosa; Festuca, Festuca contracta; Luzula, Luzula crinita; Pleuro, Pleurophyllum hookeri; Poafol, Poa foliosa; Stilbo, Stilbocarpa polaris.

The most useful variables for modelling distributions depended on the taxon. Elevation and distance from the coast are important for several species including A. magellanica, A. macquariensis, F. contracta, L. crinita, P. foliosa and S. polaris (Fig. 4, Table S5). Several species show a geographic bias on longitudinal or latitudinal axes. For example, Pleurophyllum hookeri has a distribution biased to the north, whereas A. macquariensis displays the opposite pattern. The short grasses D. cespitosa and F. contracta are biased toward the east (lee) side of the island. Proximity to the coast and topographic position had some influence on the core range of S. polaris. Aspect (eastness and northness) featured strongly in the models for D. cespitosa but had low importance for other species. Wind exposure had low importance for all models (Fig. 4, Table S5).

Fig. 4.  Variable importance for species core range models. Mean variable importance with standard deviation from 100 random forest models is shown. Coast dist, distance from coast; Solarad, solar radiation; tpi180, Topographic Position Index within 180-m radius.

Core range as a proportion of total range varied greatly, from 12.5% for S. polaris to 74% for A. macquariensis (Table 2). Most of the 10 dominant species considered here potentially co-occur across large parts of Macquarie Island (Fig. S2a). However, the extent of overlapping core ranges is much lower, ranging from eight of the nine taxa in some coastal situations and declining with elevation (Fig. S2b). The core ranges of some species coincide strongly, for example S. polaris is almost always associated with P. foliosa (Table S6).

Table 2.  Modelled habitat area of key Macquarie Island plant species.

Combining range maps into a species assemblage classification produced an original map with nine classes plus 3.4% of the island in an unclassified category (Fig. 5a, Table 3). Applying the assemblage ruleset to total range maps instead of core ranges resulted in a very different spatial pattern with the taller vegetation (three assemblages) and short grassland–herbfield complex covering most of the island, with dramatic reductions in feldmark, P. hookeri herbfield and Acaena herbfield (Fig. 5b, Table 3).

Fig. 5.  Species assemblages: (a) based on core ranges of key species; (b) based on total ranges of key species.

Table 3.  Total area (ha) classified in each floristic assemblage using either modelled core ranges or total ranges.

By using modelled core range distributions our research provides the ability to map the current species assemblages across most of Macquarie Island. Generalist species on Macquarie Island had much narrower core ranges than full ranges, consistent with the model of Austin and Smith (1989). Core ranges may be further narrowed by interspecific competition or past grazing impacts. The low number of co-occurring core ranges across most of the island suggests that these effects may be strong. The assemblage ruleset gives precedence to the taller growing species. When applied to the core range models, where these species have more than 25% cover they were assumed to dominate over co-occurring short plants which also have more than 25% cover.

The vegetation distribution modelled on core ranges indicates tussock grassland dominating at lower elevations, with P. hookeri herbfield and feldmark at higher elevations. This pattern approximates descriptions from the 19th century when tall vegetation dominated by tussocks and megaherbs extended to the lower altitudinal limit of the feldmark (Kirkpatrick 2009).

A more complex pattern with extensive areas of short grassland and grassland–herbfield and very limited extent of feldmark results from assemblages derived from the total range models. Applying the assemblage ruleset to the total range of all species assumes that the tallest species will potentially be dominant across their entire range. This outcome could happen where a species has an abrupt range boundary due to interspecific competition or a sharp environmental discontinuity.

Although both vegetation maps indicate the dominance of the large megaherb and tussock species at lower elevations, the assemblages differ in composition with varying combinations of P. foliosa, P. hookeri and S. polaris. The very limited core range of S. polaris likely reflects the relative scarcity of this species at the time of sampling due to rabbit impacts (Whinam et al. 2014). The large fern Polystichum vestitum was not modelled because of lack of data, yet historically this species formed a distinct fernbrake community in places and was sometimes co-dominant with the mergaherb and tussock species (Selkirk et al. 1990). Therefore, we expect that the mapped lowland megaherb and tussock vegetation may be complex and changeable in composition. This could be examined by further sampling of these environments.

Shorter vegetation, largely confined to middle to high elevations, indicates either abiotic limits on the dominance of larger plants or the influence of disturbance favouring particular species such as small graminoids and herbs. Acaena herbfield is possibly a grazing-induced disclimax. On Grande Terre in the Kerguelen Islands, extensive Acaena magellanica herbfields are attributed to rabbits (Lebouvier and Frenot 2007). Without grazing this community may be replaced by taller vegetation, as is indicated by our total range map (Fig. 5b), or by a more diverse assemblage of short grasses and herbs.

Similarly, short grassland and short grassland and herbfield complex probably represent a grazing-induced disclimax in many places. However, the gradient theory of Austin and Smith (1989) suggests they may also have a natural niche in the transition between highly productive and highly stressed environments. Two such ecotones might occur on Macquarie Island: permanently wet soils between the well-drained tussock grass–tall herbfield and waterlogged mires, and the transition from lowland closed vegetation to feldmark.

In the latter example, this topographic partitioning of vegetation types is consistent with the stress gradient hypothesis, where interaction between plant species is increasingly positive (facilitative) as abiotic stress increases, as demonstrated in a similar environment by vegetation patterns along an altitudinal gradient on subantarctic Marion Island (le Roux and McGeoch 2008). Azorella grows slowly and does not tolerate shade and therefore is excluded by competition except where harsh conditions preclude taller growing plants (Bergstrom et al. 1997).

Some areas were undefined by the assemblage ruleset, suggesting that these areas either do not fit into the community classification or the species models are underestimating the core ranges. The main gap in our analysis is the exclusion of mire species such as Montia fontana, Juncus scheuchzerioides and Isolepis aucklandica because of insufficient data. Mires occur on the coastal terrace and inland valleys where the watertable is constantly high. Pleurophyllum can occur in mires (Rich 1996), so the core map is not entirely inaccurate in mapping these mires as herbfield. Nevertheless, extensive areas of coastal terrace mires were wrongly classified as being dominated by S. polaris and P. foliosa in the present study. Bricher (2012) noted the same problem with a similar spatial modelling approach which neglected mire species. An additional problem is the microtopographic scale at which herbfield and mire are differentiated on the coastal terraces (Rich 1996).

Well-constructed SDMs with good predictive accuracy can discriminate variable importance in most cases and provide inferences about the relative significance of environmental factors to explore in experimental studies. This is despite complications due to small sample sizes, variable collinearity and variation in species niche width (Smith and Santos 2020). A key consideration is the interpretation of causation (as opposed to correlation), particularly when using indirect gradients (e.g. elevation) as surrogate measures for direct drivers of species distributions (Elith and Leathwick 2009; Mod et al. 2016).

Elevation is an important predictor variable for several species in this study, yet given the correlation between elevation and several direct gradients that pertains in this and other studies, it is difficult to interpret which direct environmental factors or combination of factors are involved in species responses to the elevation gradient. For example, relative humidity and windspeed increase with elevation on Macquarie Island, whereas soil temperature and air temperature decrease (Tweedie 2000; Fitzgerald and Kirkpatrick 2020). Soil freeze–thaw cycles influence plant growth at higher elevations (Boelhouwers et al. 2003; Selkirk-Bell and Selkirk 2013) and frequency of fog increases with elevation (Fitzgerald and Kirkpatrick 2020).

Wind exposure is an important climatic influence on vegetation in the subantarctic, particularly in determining the distribution of feldmark (Fitzgerald and Kirkpatrick 2017; Momberg et al. 2021a), yet wind exposure was not important for core or total range models in our analysis. This contrasts with plant species modelling for subantarctic Marion Island, where wind stress was an important factor in species distribution and cover (Momberg et al. 2021b), perhaps reflecting different approaches to modelling of wind stress.

Indirect gradients can involve biotic interactions, such as the latitudinal variation in Azorella on Macquarie Island, which is likely to be driven by lower frequency of freezing days in the north and consequent increase in pathogen-induced mortality compared to the south (Dickson et al. 2021). Shifts between interspecific competition and facilitation on gradients of environmental stress (le Roux and McGeoch 2008; Chardon et al. 2020) may also influence species distribution patterns.

Environmental factors not included in the RFMs include geology and soils. Geology does not appear to influence plant species dominance on Macquarie Island (Adamson et al. 1993). Although spatial data for soils arenot available, the topographic variables used in our models may provide some surrogacy for soil characteristics (e.g. topographic position and wetness). Satellite derived variables, particularly texture metrics (Murray et al. 2010; Bricher et al. 2013), may be useful for detecting some vegetation types, although obtaining suitable satellite imagery covering Macquarie Island has proven problematic.

We demonstrate a simple variation of the typical binary presence or absence SDMs by adding a third class representing the core range of a species. This approach is useful on Macquarie Island where most of the dominant plant species co-occur across a wide range of environments, yet the vegetation varies considerably in individual species abundance and, consequently, in community structure. Mapping of structural vegetation types can be derived from core range models, whereas absolute ranges are not useful for this purpose.

The simple field data collection and GIS-derived topographic variables make this method repeatable and data-driven. Although models could be improved with more data points from more extensive fieldwork, for most species the model performance is sufficient for repeat modelling to monitor changes in the extent of communities. Topographic variables are readily available and produce effective models.

The distribution of plant species assemblages on Macquarie Island reflects environmental gradients and may be locally modified by other factors such as grazing history. We expect that plant species distributions on Macquarie Island are still shifting as species approach their climatic limits following release from grazing pressure. Data points capturing recent expansion in species distributions would allow the models to be updated to monitor this change.

Supplementary material is available online.

Datasets are archived with the Australian Antarctic Data Centre (http://dx.doi.org/doi:10.26179/m528-md85).

The authors declare that they have no conflicts of interest.

This research was supported by an Australian Government Australian Postgraduate Award Scholarship to N. B. Fitzgerald. This research did not receive any additional specific funding.