Species-energy relationships of indigenous and invasive species may arise in different ways – a demonstration using springtails

Although the relationship between species richness and available energy is well established for a range of spatial scales, exploration of the plausible underlying explanations for this relationship is less common. Speciation, extinction, dispersal and environmental filters all play a role. Here we make use of replicated elevational transects and the insights offered by comparing indigenous and invasive species to test four proximal mechanisms that have been offered to explain relationships between energy availability, abundance and species richness: the sampling mechanism (a null expectation), and the more individuals, dynamic equilibrium and range limitation mechanisms. We also briefly consider the time for speciation mechanism. We do so for springtails on sub-Antarctic Marion Island. Relationships between energy availability and species richness are stronger for invasive than indigenous species, with geometric constraints and area variation playing minor roles. We reject the sampling and more individuals mechanisms, but show that dynamic equilibrium and range limitation are plausible mechanisms underlying these gradients, especially for invasive species. Time for speciation cannot be ruled out as contributing to richness variation in the indigenous species. Differences between the indigenous and invasive species highlight the ways in which deconstruction of richness gradients may usefully inform investigations of the mechanisms underlying them. They also point to the importance of population size-related mechanisms in accounting for such variation. In the context of the sub-Antarctic our findings suggest that warming climates may favour invasive over indigenous species in the context of changes to elevational distributions, a situation found for vascular plants, and predicted for springtails on the basis of smaller-scale manipulative field experiments.


Mechanism Synopsis
Time for speciation Longer time periods provide more opportunity for speciation.
Diversification rate Increased energy produces faster speciation or slower extinction rates.
Niche breadth Higher energy results in greater abundance of preferred resources, a switch away from non-preferred ones, reduction in niche overlap, lower competition, and thus greater richness.
Niche position Higher energy increases the abundance of rare resources and niche position resource specialists, leading to higher richness.
More trophic levels Increased energy enables additional trophic levels to occur that are occupied by new consumer species so increasing richness.
Consumer pressure As a consequence of other mechanisms, consumers are more abundant or diverse, so reducing prey populations and promoting co-existence, resulting in higher richness.

Sampling
Higher energy results in greater numbers of individuals, and random selection from a regional species pool with larger numbers of individuals results in an increased number of novel species in a focal assemblage.
Increased population size/ more individuals Higher energy areas support more individuals, leading to lower extinction rates, and thus greater numbers of species.
Dynamic equilibrium Increased energy enables faster recovery rates from disturbance, reducing the time during which small population size-associated stochastic extinction is likely to occur, hence elevating richness.
Range limitation As solar energy increases, climatic conditions are within the physiological tolerance range of more species. Table 1. Mechanisms underlying relationships between energy and species richness based on two recent theoretical treatments 12,20 . even surface estimates are unavailable at the resolution required, especially for understanding seasonal variation. Hence, we have used temperature as a proxy for energy availability. Because both geometric constraints 36 and the species-area relationship 37 should be considered a priori in the context of elevational variation in richness, we also test for these effects, examining the latter together with temperature variation as an explanation for variation in richness. We then test explicitly four of the eight proximal explanations for richness-energy relationships set out by Evans et al. 12 and specifically the sampling, more individuals, dynamic equilibrium, and range limitation hypotheses. The niche breadth, niche position, more trophic levels, and consumer pressure explanations were not considered. No evidence for interspecific competition exists for the group on the island 38 , and none for specialization -reflected by the ease with which both indigenous and invasive species can be reared using food sources from elsewhere than on the island 39 and by the broad habitats of the invasive species in their home ranges 33 , thus excluding the niche position and breadth mechanisms. The more trophic levels mechanism was excluded because all of the species belong to the same trophic level 39 . The consumer pressure hypothesis was not considered because spiders are the only major predators of springtails on the island, and the density of spiders is higher in areas with greater thermal energy availability 40,41 -areas where springtail abundance is also higher 38 , contrary to the mechanism.

Methods
Study area, species and sampling approach. Sub-Antarctic Marion Island (MI) (46°54′S, 37°45′E) is the larger (300 km 2 ) of two biologically very similar islands in the Prince Edward Island group. This sub-Antarctic, Indian Ocean island has a cool, wet, windy climate, which varies considerably with elevation (see Chown and Froneman 32 for review). Two major biomes are present: sub-Antarctic tundra, which predominates in lowland areas, and sub-Antarctic polar desert, restricted to high elevations. Sixteen species of springtails have been recorded from the island, of which six are introduced and invasive 33,38 . Here, invasive species mean those alien species which have colonised the entire lowland extent of the island, and thus are not dispersal limited (see evidence from Gabriel et al. 38 ; Hugo 42 ). The invasive species are European and in their native ranges occur across the full range of temperatures found on the island 32,43,44 . Species composition and abundance were investigated along two altitudinal transects, one on the eastern side of the island and one on the west (Fig. 1). Because some springtail species on Marion Island show pronounced seasonality in abundance 45 , both transects were sampled twice: once in winter (June-July 2008) and once in summer (February 2009). The habitat complexes see 32 sampled are representative of those found at the specific altitudes: biotically influenced vegetation (10 m a.s.l.), mires (50 m a.s.l. and 200 m a.s.l.), fellfield (400 m a.s.l. and 600 m a.s.l.) and polar desert (750 m a.s.l., 850 m a.s.l. and 1000 m a.s.l.), resulting in eight sites per transect ( Fig. 1). At each site, the GPS coordinates and the habitat type and/or plant species composition were recorded.
At each site, 20 samples were taken either using a 34 mm diameter soil corer or, for the polar desert, by removing the top 5 cm of soil and stones from a 10 cm × 10 cm area. The sampling design varied between the lower five sites and the upper three sites because in the upper sites no vegetation is present and the scoria (loose gravel) substrate cannot be sampled using a standard soil corer (full sampling descriptions are provided in the Supporting Information Appendix S1). Previous work has found comparable results using these approaches 38 .
For each of the samples taken at the five lower elevation sites, springtails were extracted from the top 10 cm of cores using a Macfadyen high-gradient extractor into 40% ethanol (extraction protocol: 2 days at 25 °C followed by 2 days at 30 °C 45 ) for the east transect. Samples from the west transect were extracted using Tullgren funnels in a field hut. The difference in extraction method was not considered a source of bias because densities from the western samples were well within those found by Hugo et al. 46 and Hugo 42 for comparable sites after high gradient extraction. For each of the samples from the three higher elevation sites, the protocol followed Barendse and Chown 45 . Material was washed (using water previously sieved through 125 µm mesh) and sieved through 125 µm mesh three times. A separate mesh was used for each sample and after sample washing each mesh was placed in a separate 35 ml plastic container with 40% ethanol. All samples were further sorted at the research  energy availability. We use temperature as a metric of energy availability. In systems which are not water limited, as is the case for this cool island which typically receives in excess of 1900 mm of precipitation per year 32 , plant productivity is constrained by the influence of temperature on both growing season length and on rates of leaf metabolism 35,47 . Thus, we consider temperature a reasonable proxy for productivity in such systems 35 . While there are other proxies for energy availability, such as solar radiation or potential evapotranspiration 35 , these have not been measured across the island. This proxy approach should be kept in mind in considering the outcomes of the tests we present here. Moreover, variation in temperature also plays a direct physiological role in ectotherm physiology 48 , making it applicable to both tests of the purely productive energy forms of species-energy mechanisms, and those that have to do with temperature variation directly 12,35 .
Soil temperature, appropriate for springtails, was measured using Thermochron iButton temperature dataloggers (Model DS1921, accurate to ±0.5 °C, Dallas Semiconductors, Dallas, TX, USA) placed ca. 2 cm below the ground surface at each site. Hourly temperatures were recorded for three winter months (mid-May 2008 to mid-August 2008) and three summer months (mid-November 2008 to mid-February 2009). Temperature data were processed to obtain mean daily minimum, maximum, mean and temperature range for each site in R (version 2.12.0) 49 (hereafter: short-term data).
A longer-term data series (hereafter: long-term data), processed for the same variables described above, was obtained by recording soil temperatures on an hourly basis at c. 100 m intervals from sea level to 750 m, across the eastern slope of Marion Island between 2002 and 2009 using iButton temperature dataloggers 50 . Missing data (due to datalogger loss, exposure, or damage) were interpolated using a sinusoidal function (i.e. to approximate daily temperature cycles, written in R2.12.0), with the initial and final interpolated values (as well as the amplitude of the sinusoidal curve) calculated from the temperature records of the 48 hours before and after the missing data.
Sampling adequacy and species richness. Raw abundance data for indigenous and invasive species were converted to density (individuals.m −2 ) to enable comparison of data from the two different sampling approaches (though we use the term abundance hereafter). To determine sampling adequacy for each site, sample-based (i.e. per core or 10 cm × 10 cm sample) rarefaction curves were calculated using the Mao Tau moment-based interpolation method in EstimateS 51 .

Richness-energy relationships, geometric constraints and area effects.
Because the level of interest here is variation among sites with different energy availability, and because multiple samples were taken to ensure adequate site sampling 52 , the site level was used as the level for investigation.
To test for geometric constraints, species richness data were compared with null model predictions using a Monte Carlo simulation procedure, Mid-Domain Null, in Visual Basic for Excel 53 . Randomisation techniques were used, based on 50 000 simulations sampled without replacement from empirical range sizes, and a regression of the empirical values on predicted values provided r 2 , slope and intercept as estimates of the fit of the null model 53 . We used the coastal and high elevation end-points for each transect as the limits to the domains, to avoid biasing assessments against a mid-domain effect, which would have been the case had we used a single domain (i.e. a coast to coast approach). Springtails have not been found above 1000 m, largely because higher elevations (the island rises to 1230 m) are represented by just a few isolated peaks 32 . www.nature.com/scientificreports www.nature.com/scientificreports/ To examine the influence of surface area on richness, we used generalised linear models (GLMs, assuming a Quasipoisson distribution to correct for overdispersion, and with a log link function, implemented using the glm function in R2.12.0 54 ) with species richness as the response variable, and surface area and mean temperature (from the short-term data) as the predictor variables to establish the role of surface area. Sites were allocated to altitudinal bands and the surface area of each of these bands extracted from Meiklejohn and Smith 55 . Full models included energy and log transformed surface area as variables, and single term models included either temperature or log(area) individually as variables. We re-assessed the relationships using generalised estimating equation models to take the potential effects of spatial autocorrelation into account 56 .
Following this assessment, which revealed relationships between temperature and richness, we compared the estimates for the richness-mean temperature relationship for the indigenous and invasive species to determine whether they differ. We used GLMs (Quasipoisson distribution, log link function) with species richness as the response variable, and mean temperature (from the short-term data) and species category (indigenous/invasive) as the predictor variables, specifically focussing on whether the interaction term is significant, indicating different forms of the relationships for indigenous and invasive species.
tests of the potential mechanisms. Following Evans et al. 's 12 recommendations, to distinguish the sampling and increased population size hypotheses, the relationships between species abundance and richness (using S obs for the sampling explanation and the Jacknife2 estimator 57 for the increased population size explanation, calculated for each of the sites using the vegan package in R2.12.0), species richness and mean temperature, and abundance and mean temperature were investigated for both indigenous and invasive species as well as the combined assemblages using GLMs (Quasipoisson distribution, log link function). To test for decelerating relationships, the GLMs were rerun using a squared term for the predictors. In all cases the units were sites on each transect, resulting in n = 16. If the relationships are not positive and decelerating using the estimator, but are for the raw data, the increased population size explanation in this specific form can be rejected 12 .
The dynamic equilibrium explanation posits that elevated energy enables populations to recover faster from disturbance. On Marion Island, large-scale disturbances are mostly from low temperature events 32,58 (Supporting Information Appendix S2). Temperature disturbances in the form of two thresholds were considered. These are 0 °C (the freezing point of water), and the mean lower development threshold (LDT) of the springtail eggs (the most sensitive developmental stage 39 ). Both the number of times the thresholds were crossed and the maximum duration spent below the thresholds at any one time (hereafter -longest duration) were assessed as disturbance events for each site using the short-term and long-term temperature data and empirical values of the springtail egg LDTs (Table S1) 39 . The effects of thresholds (number of events and longest duration of sub-zero events and those below the LDT) on the abundance (summed density at each site) of indigenous and invasive  www.nature.com/scientificreports www.nature.com/scientificreports/ species were investigated using generalised linear models (in R3.1.0, Quasipoisson distribution, log link function; Deviance Explained (DE) was calculated using the BiodiversityR package). Minimum adequate model selection was adopted 54 .
The range limitation explanation proposes that species occur in areas where they can meet their physiological requirements 12 , and assumes that more species can tolerate warmer than cooler conditions 9 . To test this idea and to distinguish possible outcomes from the dynamic equilibrium explanation, the summed abundances of indigenous and invasive species groups were compared to the minimum and maximum temperatures recorded at each site (i.e. winter mean daily minimum and summer mean daily maximum), and to the number of generations possible at each elevation. Conditions that do not exceed springtail tolerances and enable multiple generations are more likely to result in positive population growth, than those that limit populations either via development or thermal tolerances 59,60 . Both the short-term and long-term temperature datasets were used (see above). The number of generations possible was calculated using the sum of effective temperatures (SET) 61 . Mean SET values were calculated for the indigenous and invasive groups of species at each site (detailed information in Appendix S2). The effects of minimum and maximum temperatures and generations possible on the abundances of indigenous and invasive species were investigated using generalised linear models (as above). Minimum adequate model selection was adopted 54 .

Results
Strong elevational declines in temperature were found ( Fig. 2; Table S2). The low temperature values, together with site-specific rainfall records indicating rainfall above 1000 mm over the year (Fig. S1), substantiated the assumption that the soil ecosystem is not water limited. Eight indigenous and four invasive springtail species were sampled on the west transect, and seven indigenous and four invasive species were found on the east. When summer and winter data were pooled, rarefaction curves indicated that sampling had reached an asymptote for most sites (Fig. S2), and was therefore considered adequate for the summed data, which was used for all further analyses (Table S3).
Fit to the geometric constraints null model was poor for observed species richness for both transects (east: r 2 = 0.003; west: r 2 = 0.064) and the relationships were not significant (east: p = 0.66, slope = 0.11, intercept = 6.92; Although both surface area and mean temperature were significant predictor variables for indigenous and invasive species richness in the single term GLMs, surface area was not significant in the full models ( Table 2) and was therefore not considered any further (though bearing in mind that these variables are collinear, r 2 = 0.63, p < 0.001). Similar outcomes were found using the GEE models, which explicitly consider spatial relationships (Table S4), though in this case significance values could not be fully relied on because Poisson rather than Quasipoisson distributions were required in the models. Positive species-energy relationships were found for the indigenous and invasive species (Table 3; Fig. 3), with a significant difference in the form of the relationship ( Table 4).
The strong positive relationships between energy and richness tended not to be present with the richness estimator and a similar effect was found for the relationship between richness and abundance, although the effect for richness was larger for the indigenous than the invasive species (Table 3). Overall, the variable positive and decelerating relationships between richness and energy, richness and abundance, and between abundance and temperature suggested that neither the sample size nor sampling explanations could account for the energy-richness relationship.
Significant relationships were found between species abundance and disturbance thresholds for the long-term, but not for the short-term datasets (Table 5). In the case of the latter, abundance of indigenous species tended to decline with the number of events below 0 °C, but weakly so, and a weak positive relationship between duration below LDT and abundance was found, with less than 50% of the deviance being explained by the model overall. For the invasive species 93% of the deviance was explained by a single negative relationship between the number of events below the LDT and abundance.
In the context of potential range limitation, no variation in abundance was explained for the indigenous species (Table 6), for either the short-term or long-term datasets, although for invasive species abundance increased with minimum temperature and accounted for much of the deviance in the data (DE = 64 to 92%) ( Table 6).

Discussion
Species richness declined monotonically in both transects, in keeping with patterns found for some, but not all elevational gradients in richness 4,37 . No evidence was present for the influence of geometric constraints, in keeping with the relatively small scale of the work 62 . Although richness was positively related both to energy and area, area tended not to be significant in the models including both variables (bearing in mind the variables are collinear). Thus, positive species-energy relationships exist for this system. Moreover, the variation in richness accounted for by mean temperature (ca. 50%) was in keeping with (or perhaps slightly higher than) that found for other organisms at this spatial extent 6,63 .
Differences between the indigenous and invasive species in the form of the species-energy relationship were significant, with a steeper relationship for the invasive species. The outcome, and especially the shallower slope of the richness-energy relationship for the indigenous species suggests that a time for speciation mechanism 4,21 might be important. In other words, occupation of indigenous springtails of multiple habitats over a prolonged period of the island's history would have enabled speciation in all of them (see discussion in Myburgh et al. 64 ), thus reducing the slope of the species-energy relationship. The lack of springtail clades which have diversified on the island suggests, however, that the mechanism is unlikely to have played a large role. Nonetheless, at least two species are endemic or suspected to be so 33,65 , indicating that local speciation has taken place.
By contrast, much of the evidence points to variation in population size as an important explanation for the richness-energy relationships, but in a manner different to that proposed directly by the more individuals hypothesis. Indeed, of the mechanisms for the species-energy relationship compiled by Evans et al. 12 that we examined, neither the more individuals hypothesis, nor some form of sampling artefact were supported. By contrast, the generalised linear models suggested that both the range limitation and dynamic equilibrium mechanisms are plausible for the richness-energy relationships, especially in the case of the invasive species. Although Evans et al. 12 rightly considered the dynamic equilibrium and range limitation mechanisms distinct from the more individuals hypothesis, given that the former have to do with productive energy and the latter may incorporate direct effects of temperature too, it is clear that they constitute a subset of a broader population size-related set of mechanisms. Indeed, they both invoke variation in population size or abundance associated with an environmental filter 4 , supporting a more general emphasis on population size variation as a factor explaining energy-richness relationships 13 .
In the case of the dynamic equilibrium mechanism, small populations are unable to recover from disturbance, and are more prone to negative effects of environmental stochasticity 66,67 which ultimately leads to local extirpation. How this might play out for springtails on the island is clear. Invasive species typically have higher lower  www.nature.com/scientificreports www.nature.com/scientificreports/ developmental thresholds for development than do their indigenous counterparts 39 , which would limit population growth at low temperature. Similarly, differences in the ability of the two groups of species to cope with temperatures below freezing have been recorded 68 , with the indigenous species capable of dealing with the most extreme events found on the island. Moreover, variation in temperatures along the altitudinal gradient shows that low temperature disturbance events are more common at higher rather than lower altitudes (see also discussion in Lee et al. 58 ).
In the case of range limitation, lack of physiological tolerance accounts for either an inability to exist at a site, or low abundances 9,15 , which in turn result in higher extinction probability. These processes mean that even if a site is within the dispersal range of a species, extinction-related processes preclude viable populations 58,69,70 . For the springtails considered here, the indigenous species show no limitation by low or high temperatures, perhaps unsurprising given their tolerance limits and the relationship between egg development and temperature 39,68 . By contrast the strong relationship between minimum temperature and invasive abundance suggests that their low temperature sensitivity 39,68 may be limiting.
Differences in the likely mechanisms underlying species-energy relationships among the indigenous and invasive species illustrate that while the patterns for these two groups of species appear similar in broad terms   www.nature.com/scientificreports www.nature.com/scientificreports/ (richness declines with decreasing energy), the explanations for them may be quite different. In this case, dynamic equilibrium mechanisms apply to both groups, while range limitation plays a larger role for the invasive than the indigenous species. Moreover, the lower variation explained in the models for the indigenous species, yet the significant, though shallow, species-energy relationship suggests that speciation-associated mechanisms may also be important for them, in keeping with studies suggesting that time for speciation effects are relevant over small spatial extents 21 . Importantly, a strong relationship between abundance and temperature for the invasive species across the gradient, but little effect on the indigenous species, adds support to the idea that environmental constraints are more important for the former group 71,72 , despite their origins in cold-temperate Europe.
While these differences might at first appear to constrain lessons from comparisons among the two groups of species, they may prove insightful by enabling various likely mechanisms in an area to be disentangled by relying on such a priori expectations of differences (see also Marquet et al. 73 ; Hawkins et al. 74 ). In the context of the specific systems of the sub-Antarctic islands, the strong influence of temperature on richness and on abundance patterns of the invasive compared with the indigenous species suggests also that rising temperatures and declining rainfall on many of the islands e.g. [75][76][77][78] will enable the spread of invasive species to higher elevations, as has already been documented for plants 79 . Field manipulations have supported an assumption of greater population-level success for invasive species at local scales under warming and drying 80 , with this broader scale analysis providing further substantiation for this idea. Nonetheless, differentiation of temperature effects from energy availability per se is still required, given that we used temperature as a proxy for the latter 35 . Only with well-developed, fine resolution data on either net primary productivity or proxies such as potential evapotranspiration or solar radiation will this be possible. While estimates of primary productivity have been made for Marion Island 81,82 , these have not been spatially explicit and remain challenging for polar desert environments that either have an interstitial flora, or high cloud cover(precluding remote-sensing estimates), or both.

Data Availability
The core level site data collected in this study are publicly available from the Monash Figshare repository (https:// doi.org/10.26180/5cf4f64629f6a).