Distinguishing effects of area per se and isolation from the sample‐area effect for true islands and habitat fragments

1051 –––––––––––––––––––––––––––––––––––––––– © 2021 The Authors. Ecography published by John Wiley & Sons Ltd on behalf of Nordic Society Oikos This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. Subject Editor: Henrique Pereira Editor-in-Chief: Miguel Araújo Accepted 21 March 2021 44: 1051–1066, 2021 doi: 10.1111/ecog.05563 44 1 0 5 1 – 1066


Introduction
From true islands to habitat fragments on terrestrial landscapes, positive relationships between species richness and the area of islands or fragments are among the oldest and most widely documented patterns in ecology (Arrhenius 1921, MacArthur and Wilson 1963, Rosenzweig 1995, Gotelli and Graves 1996, He and Legendre 2002, Hanski et al. 2013). These island species-area relationships (ISARs, sensu Triantis et al. 2012, or 'Type IV' curves, sensu Scheiner 2003 have received considerable attention, in part due to their importance to conservation frameworks (reviewed by Shafer 1990, Lomolino 2000, Whittaker and Fernández-Palacios 2007. Although there are several documented divergences between the biogeographies of true islands and terrestrial habitat fragments (Laurance 2008, Mendenhall et al. 2014, Itescu 2019, Farneda et al. 2020, studies addressing trueisland systems can still help resolve what mechanisms underlie ISARs and how fragmentation effects are best measured (Diamond 1975, Simberloff and Abele 1976, 1982, Haila 2002, Haddad et al. 2015, MacDonald et al. 2018a. There are, however, at least two enduring problems with the use of ISARs in conservation that are generalizable to both true islands and habitat fragments: 1) ISARs may emerge from a combination of different mechanisms, each of which potentially informs a different conservation directive (Connor and McCoy 1979, Kadmon and Allouche 2007, MacDonald et al. 2018b); and 2) ISARs are emergent patterns of diversity that can mask differential responses to habitat area among species that may require independent consideration for successful conservation planning (Ewers and Didham 2006, Öckinger et al. 2009, Franzén et al. 2012, Hanski 2015, MacDonald et al. 2018a. Here, we propose an approach to addressing each of these problems within a single modelling framework.

Mechanisms underlying ISARs
Three hypotheses, each with distinct underlying mechanisms and conservation implications, have been proposed to account for ISARs and related spatial patterns of species richness: 1) the passive sampling hypothesis (Connor and McCoy 1979); 2) area per se, as outlined by the theory of island biogeography (MacArthur and Wilson 1963, Wilson andMacArthur 1967); and 3) the habitat diversity hypothesis (Williams 1964). The passive sampling hypothesis, originally developed within the context of oceanic islands, predicts that islands randomly sample individuals from the regional species pool in abundances proportional to their area (Connor and McCoy 1979). As larger islands sample more individuals, they sample more species according to the abundance distribution of the regional species pool (i.e. the 'sample-area effect'). Thus, passive sampling serves as a useful null hypothesis, assuming random assembly of both individuals and species.
Area per se hypothesizes a disproportionate reduction in species richness as island area decreases, steepening the slope of ISARs within archipelagoes or fragmented terrestrial landscapes relative to species-area relationships within landscapes comprised of continuous habitat (Diamond 1972, 1975, Wilson and Willis 1975, Connor and McCoy 1979, Saccheri et al. 1998, Gonzalez 2000, Haila 2002, Haddad et al. 2015, MacDonald et al. 2018a. From a mechanistic perspective, area per se essentially invokes the theory of island biogeography, where species richness arises as a dynamic equilibrium between rates of extinction and colonization, which in turn depend on island area and isolation (MacArthur and Wilson 1963, Wilson andMacArthur 1967). More isolated populations occupying small islands are predicted to be more prone to stochastic extinction and small, isolated islands are less likely to be re-colonized from external source populations than larger, well-connected islands (Levins 1969, Hanski and Gyllenberg 1993, Orrock and Watling 2010. Thus, the theory of island biogeography addresses effects of both island area and isolation, predicting that immigration rates and rescue effects decrease as islands become increasingly isolated from neighboring habitat (i.e. the mainland or other islands), negatively affecting species' probabilities of occurrence and therefore species richness (MacArthur and Wilson 1963, Wilson and MacArthur 1967, Brown and Kodric-Brown 1977, Hanski 1994, 1999.
As an alternative to dynamic balances between demographic rates predicted by the theory of island biogeography, Williams (1964) proposed that ISARs are driven by variation in habitat diversity among islands. The habitat diversity hypothesis predicts that island/fragment area correlates with species richness only insofar as area correlates with the intermediate variable of habitat diversity; larger sample areas generally contain more habitats, which support more species (Williams 1964, Nilsson et al. 1988, Rosenzweig 1995, Gotelli and Graves 1996. It follows that the presence or proportion of suitable habitat within islands/fragments should affect abundances and occurrences of individual species (Buckley 1982, Haila andJärvinen 1983). However, few studies have investigated how habitat associations of single species relate to emergent patterns of species richness in this context (but see Haila et al. 1983). Still, multiple studies addressing species richness support the habitat diversity hypothesis (Nilsson et al. 1988, Kadmon and Allouche 2007, Hortal et al. 2009). It has also been inferred that mechanisms predicted by the passive sampling hypothesis, theory of island biogeography and habitat diversity hypothesis may simultaneously contribute to ISARs (Connor and McCoy 1979, Kadmon and Allouche 2007, MacDonald et al. 2018b, Chase et al. 2019.

Complications in extending ISARs to conservation
Due to its success on true islands, conservation biologists were quick to recognize the potential for the theory of island biogeography to be applied to fragmented habitat on terrestrial landscapes, initially in the design of nature reserves (Diamond 1972, 1975, Wilson and Willis 1975. However, the extension of ISAR-based inferences from true islands to habitat fragments is complicated by differences in their extents of insularity and interactions between habitat and non-habitat (i.e. matrix) areas (Itescu 2019). Whereas edges of true islands clearly delimit suitable habitat from a homogeneous matrix of unsuitable habitat (open water), species occurring on habitat fragments may utilize resources of heterogeneous terrestrial matrices and these matrices may differentially constrain or facilitate colonization rates of species (i.e. 'matrix effects'; Dunning et al. 1992, Ricketts 2001. Thus, understanding habitat fragments as analogous to true islands may be problematic. More specifically, fragmentation effects observed within true-island systems may differ from those typical of fragmented terrestrial landscapes (Laurance 2008, Mendenhall et al. 2014, Farneda et al. 2020.

SLOSS-based inferences
Irrespective of the mechanisms underlying ISARs, observations that species richness generally decreases as island/fragment area decreases and isolation increases have contributed to long-standing inferences that habitat fragmentation poses a major threat to species diversity (Diamond 1972, 1975, Noss 1991, Haila 2002, Hanski 2015, Fletcher et al. 2018. However, many of these inferences are founded on observations or experimental designs that have not sufficiently decoupled the effects of area per se and isolation from the sample-area effect (i.e. decoupled habitat fragmentation from habitat loss; sensu Fahrig 2003, Hadley and Betts 2016. In the majority of studies successfully decoupling habitat fragmentation from habitat loss, single large habitat fragments are generally found to contain an equivalent or lesser number of species than sets of several small habitat fragments summing to an equivalent total area (Quinn and Harrison 1988, Fahrig 2003, Yaacobi et al. 2007, MacDonald et al. 2018a, b, Deane et al. 2020. Such comparisons contribute to the ongoing single-large-or-severalsmall ('SLOSS') debate, addressing how finite conservation efforts should prioritize the area and configuration of fragmented habitat and nature reserves (Diamond 1975, Abele and Connor 1979, Ovaskainen 2002, Tjørve 2010. In light of SLOSS-based observations that species richness is often equal or greater within sets of several small habitat fragments, Fahrig (2013) advanced the habitat amount hypothesis, predicting that the number of species persisting on fragmented landscapes is only a function of total habitat amount at the landscape scale irrespective of its spatial subdivision and configuration. The principal mechanism underlying the habitat amount hypothesis is the sample-area effect, as originally articulated by the passive sampling hypothesis (Connor and McCoy 1979). However, the habitat amount hypothesis extends implications of the sample-area effect to predict that there should also be no detectable effect of fragment isolation on species' abundances, species' occurrences or species richness after the sample-area effect has been accounted for (Fahrig 2013).

The importance of understanding how species-level patterns affect ISARs
While the sample-area effect surely contributes to patterns of species richness within a variety of true-island systems and fragmented terrestrial landscapes, richness-based analyses may obscure important effects of both area per se and isolation on individual species (Ewers and Didham 2006, Öckinger et al. 2009, Franzén et al.2012, Hanski 2015, MacDonald et al. 2018a. Indeed, area per se and isolation effects have been inferred to vary widely among species, even within single landscapes and taxa (Henle et al. 2004, Ewers and Didham 2006, Nowicki et al. 2007, Öckinger et al. 2009, Hanski 2015, Hillebrand et al. 2018, MacDonald et al. 2018a. Functional traits, including body size (Gehring and Swihart 2003, Henle et al. 2004, Larsen et al. 2008, Prugh et al. 2008, Barbaro and Van Halder 2009, Warzecha et al. 2016, mobility/dispersal ability (Roland and Taylor 1997, Lens et al. 2002, Ewers and Didham 2006, Öckinger et al. 2009, MacDonald et al. 2018a, degree of ecological specialization (Tscharntke and Brandl 2004), rarity/conservation status (Ewers and Didham 2006) and trophic position (Tscharntke et al. 2002, Thies et al. 2005, Ewers and Didham 2006 are hypothesized to relate species' sensitivity to area per se and isolation. Still, relatively few empirical studies have investigated how functional traits relate to interspecific variation in responses to area per se and isolation or how this interspecific variation scales to emergent patterns of species richness, such as those reflected in ISARs (Melbourne et al. 2004, but see Barbaro andVan Halder 2009, Öckinger et al. 2009).

Distinguishing mechanisms underlying ISARs
Due to the sample-area effect, observations that species' abundances, species' probabilities of occurrence or species richness positively correlate with island/fragment area cannot be directly interpreted as evidence of effects of area per se (Connor and McCoy 1979, Fahrig 2003, Fletcher et al. 2018. Three established methods may be used to control for the sample-area effect: 1) comparing sets of islands/fragments that sum to equal areas but differ in degree of fragmentation, including the nested-set designs of Yaacobi et al. (2007) and MacDonald et al. (2018a, b) and comparisons of species accumulation curves proposed by Quinn and Harrison (1988); 2) extrapolating a speciesarea regression to the total area of all islands/fragments used to build the regression and comparing predicted and observed species richness (γ-diversity) (Rosenzweig 2004, Yaacobi et al. 2007, Santos et al. 2010, Gavish et al. 2012, MacDonald et al. 2018a; and 3) comparing equal-area sampling plots across islands/fragments (Westman 1983, Stevens 1986, Quinn et al. 1987, Kelly et al. 1989, Fahrig 2013, Watling et al. 2020). However, for methods 1 and 2 (SLOSS-based comparisons), substantial species turnover among several small islands/fragments can increase their aggregate richness relative to single large islands/fragments, such that important effects of area per se on individual species are obscured (sensu Simberloff 1976, MacDonald et al. 2018b, Deane et al. 2020. Assuming species are uniformly distributed within islands/fragments, inferences drawn from method 3 may be robust for sessile taxa (e.g. vascular plants; Westman 1983, Quinn et al. 1987, Kelly et al. 1989), but remain tenuous for vagile species that move frequently within islands/fragments, as individual sampling plots may accumulate all vagile species within single islands/fragments if sampling effort is high. Additionally, if rare species are particularly sensitive to area per se or isolation, these effects are likely to go undetected when sampling at small spatial grain sizes dictated by method 3; rare species and important habitat types within islands/fragments may be missed entirely in comparisons of small sampling plots (Karger et al. 2014). Finally, each of the three established methods only contrast the sample-area effect with those of area per se (methods 1 and 2) or area per se and isolation (method 3); evaluating effects of isolation and habitat diversity requires additional analyses. A fourth method, recently proposed by Chase et al. (2019), employs parameters derived from individual-based rarefaction curves across various spatial scales within islands/fragments to distinguish between the sample-area effect and effects of area per se and habitat diversity. While this framework can effectively resolve mechanisms underlying ISARs, it cannot simultaneously evaluate effects of isolation, assess whether area per se and isolation differentially affect species in relation to their functional traits, or be applied to existing datasets lacking abundance data for subplots stratified within each island/ fragment across diagnosable habitat heterogeneity.
In this paper, we present a novel application of random placement and linear mixed effects models that can simultaneously evaluate: 1) how area per se, isolation, and habitat diversity affect species' abundances, species' occurrences and species richness across true islands or terrestrial habitat fragments; and 2) whether interspecific variation in these responses relates to variation in species' functional traits. This modelling framework is applicable to any dataset for which abundance data were collected for sets of true islands or habitat fragments with sampling effort standardized per unit area. We assess the utility of the framework using a butterfly assemblage persisting on a naturally fragmented landscape of true islands; Lake of the Woods, Canada. Methodological developments and basic inferences presented here are equally applicable to both true islands and terrestrial habitat fragments, so long as the edges of habitat fragments can be consistently delimited.

Overview of the modelling framework
Starting with the assumption that all individuals of each species are randomly distributed across true islands or habitat fragments in abundances proportional to their areas, random placement models can be used to calculate expected species' abundances, expected probabilities of species' occurrences and expected species richness for each island or fragment (Arrhenius 1921, Coleman 1981, Gotelli and Graves 1996, He and Legendre 2002. Resulting random placement values are equivalent to values of species' abundances, species' probabilities of occurrence and species richness predicted by the sample-area effect. Predictions of the passive sampling/ habitat amount hypotheses, theory of island biogeography and habitat diversity hypothesis may be then simultaneously evaluated by modelling relationships between random placement residuals (observed values minus random placement values) and the area, isolation and habitat diversity of individual islands/fragments. Variables measuring species' functional traits may be introduced to abundance and occurrence models via interaction terms with island/fragment area and isolation to evaluate whether effects of area per se and isolation interspecifically vary contingent on the measured traits.

Random placement models
We present random placement models in ascending order of mathematical complexity, from species' abundances, to species' occurrences, to species richness. Within random placement models, a j is the area of the jth island/fragment, A T is the total area of all islands/fragments, n i is the abundance of species i summed across all islands/fragments and S is the total number of species observed. Islands or fragments that were not surveyed are not included in random placement models. According to the sample-area effect, the expected abundance of species i on island/fragment j is simply proportional to j's area (model 1) and the occurrence probability follows the random placement model (model 2). The random placement model for expected richness on any island/fragment is simply the sum of the expected probabilities of occurrence over all species (model 3). This random placement richness model was first proposed a century ago by Arrhenius (1921) and later reinvented by Coleman (1981) with the inclusion of the variance.

Modelling of random placement residuals
Subtracting random placement abundance and occurrence probability values from observed abundance and occurrence values for each species on each island/fragment produces abundance and occurrence residuals. The direction and magnitude of these residuals measure how abundances and occurrences of each species deviate from predictions of the sample-area effect. In the modelling framework described below, all species' abundance residuals and all species' occurrence residuals are concatenated across species for use in single linear mixed effects models; one model addressing all species' abundances and one model addressing all species' occurrences. Abundance residuals require standardization (subtracting the mean and dividing by standard deviation) before concatenation across species, as the range of possible values is greater for common species than rare species. While this is not the case for occurrence residuals (values bound between −1 and 1), standardization is still recommended to generate values that are commensurate among species. Richness residuals are similarly calculated for each island/fragment by subtracting random placement richness values from observed richness values. Standardization is recommended to facilitate comparisons of effect sizes among abundance, occurrence and richness analyses.

Species' abundances and occurrences
Linear mixed effects models are used to relate abundance and occurrence residuals to island/fragment variables while controlling for species identity as a random effect. If area per se affects species' abundances or occurrences beyond variation associated with the sample-area effect, residuals will be positively related to area, indicating a disproportionate concentration of species' abundances or occurrences on larger islands/ fragments. If isolation negatively affects species' abundances or occurrences, residuals will be negatively related to measures of isolation specific to individual islands/fragments, indicating a disproportionate concentration of species' abundances or occurrences on less isolated islands/fragments. Each of these results align with predictions of the theory of island biogeography. Alternatively, the absence of significant relationships between residuals and area or isolation would indicate that the sample-area effect sufficiently accounts for variation in species' abundances or occurrences across islands/ fragments of varying area or isolation. The combination these results would confer support for the passive sampling/habitat amount hypotheses. The proportion of species-specific suitable habitat and presence of species-specific resources within each island/fragment may also be included in linear mixed effects models. Positive relationships between abundance or occurrence residuals and these variables would indicate that availability of specific habitats or resources within islands/ fragments are important considerations that affect species' abundances or occurrences, as predicted by the habitat diversity hypothesis. If data on species' functional traits are available, functional trait variables may be introduced to linear mixed effects models via interaction terms with island/fragment area and isolation. Here, a significant interaction between a functional trait variable and island/fragment area or isolation would indicate that area per se or isolation differentially affects species' abundances or occurrences contingent on the measured trait. Total abundance and number of occurrences (prevalence) for each species are also of interest, as rarity is cited as a predictor of species' sensitivity to fragmentation (Ewers and Didham 2006). A significant positive interaction between total abundance or prevalence and island/fragment area would indicate that rare species are disproportionately concentrated or likely to occur on larger islands/fragments. Similarly, a significant negative interaction between total abundance or prevalence and island/fragment isolation would indicate that rare species are disproportionately concentrated or likely to occur on less isolated islands/fragments.

Species richness
A similar modelling process may be applied to species richness using linear models. Significant relationships between richness residuals and island/fragment area or isolation would indicate that area per se or isolation significantly affects richness after controlling for the sample-area effect. Each of these results align with predictions of the theory of island biogeography. Conversely, the absence of significant relationships between residuals and area and isolation would indicate that the sample-area effect sufficiently accounts for variation in species richness across islands/fragments of varying area or isolation. This combination of results would suggest that only habitat amount at the archipelago-or landscape-scale affects richness, as predicted by the passive sampling/habitat amount hypotheses. Predictions of the habitat diversity hypothesis may also be simultaneously evaluated by including measures of habitat diversity in linear models. Significant relationships between richness residuals and habitat diversity would indicate that, despite correlations between habitat diversity and island/fragment area, variation in habitat diversity among islands/fragments affects species richness beyond variation associated with both the sample-area effect and effects of area per se.

Study area
We assessed the utility of this modelling framework using a butterfly assemblage persisting on a ~1250 km 2 lake-island complex located in Sabaskong Bay, Lake of the Woods, Canada (Fig. 1). Differential isostatic rebound and outlet restriction resulted in the flooding of the study area and isolation of land-bridge islands approximately 3000-4000 YA (Yang and Teller 2005). Given this substantial time-sinceisolation, we infer that species assemblages have relaxed to equilibria (> 1000 generations; based on categories suggested by Fahrig 2020). Butterflies represent a suitable taxon for this investigation, as most species complete their life cycles within relatively small patches of habitat, their detectability is generally high and their diversity correlates with that of many other terrestrial taxa (Thomas 2005, Nowicki et al. 2008, MacDonald et al. 2017, 2018a. Butterflies do not utilize open water at any life stage, meaning the matrix separating islands in this system is entirely unsuitable. This effectively controls for matrix effects (Dunning et al. 1992, Ricketts 2001), but also limits the generalizability of our inferred fragmentation effects to terrestrial landscapes (Laurance 2008, Mendenhall et al. 2014, Farneda et al. 2020. Thirty islands, ranging in area from 0.1 to 8.0 ha, were randomly selected from lists of candidate islands compiled according to the methods of MacDonald et al. (2018a). Islands were only considered as candidates if they were isolated from other landmasses by at least 100 m, beyond the inferred visual ranges of butterflies (Rutowski 2003, MacDonald et al. 2019. The relative isolation of each study island was quantified at multiple scales as the proportion of water within 250-, 500-, 1000-, 1500-, 2000-, 2500-, 3000-, 3500-, 4000-, 4500-and 5000-m buffers. Buffers were generated from island edges, meaning isolation measures are independent from and uncorrelated with island area. We opted for these proportion-based measures because they have been shown to predict immigration rates, rescue effects and related ecological processes more accurately than distance-based measures (Moilanen and Nieminen 2002, Tischendorf et al. 2003, Prugh 2009). Habitat diversity was estimated on each island as the relative proportion of 14 habitat types, defined using structural properties of vegetation and geological features (Supporting information for habitat type descriptions).

Survey protocol
Butterfly abundance data were collected for each of the 30 islands using repeated full-island surveys. Each island was visited four times at intervals between 10 and 14 days during the peak flight season (1 June 2015 to 20 Aug 2015). Sampling effort was standardized to 40 min per ha per survey across all islands, eliminating the need for sampling effort and diversity corrections (e.g. rarefaction or extrapolation, Chao et al. 2014, Fahrig 2020. Care was taken to visit all habitat types during each survey and handheld GPS units were used to ensure uniform coverage of islands. Recording observer tracks and capturing individuals whenever possible (kept as voucher specimens or released at the end of each survey) minimized the possibility of double counts, where individuals are recorded multiple times in single surveys. To ensure optimal and standardized butterfly activity, surveys were restricted to the hours of 10:45 to 15:45 and were not conducted in wind speeds over 15 km h −1 or in temperatures below 13°C. If temperatures were below 17°C, surveys were only conducted in sunny conditions (cloud cover < 40%). Surveys were conducted in temperatures above 17°C, regardless of cloud cover (MacDonald et al. 2017). Figure 1. Map of the study area, located in Sabaskong Bay, Lake of the Woods, Canada. Butterfly abundance, occurrence and species richness data were collected for 30 study islands, varying in area from 0.09 to 8.4 ha, using repeated full island surveys.
The diversities of vascular plants and butterflies have been observed to positively correlate with one another in a variety of systems, primarily due to larval host plant associations (Erhardt 1985, Sparks and Parish 1995, Simonson et al. 2001, Croxton et al. 2005, Nowicki et al. 2007, Kitahara et al. 2008, MacDonald et al. 2018a, Riva et al. 2020. Accordingly, vascular plant species richness was quantified on each island using repeated full-island surveys (four surveys total), standardized to 40 min per ha per survey (MacDonald et al. 2018b for further details on vascular plant surveys). A total survey time of two hr and 40 min per ha is consistent with recent recommendations for boreal plant communities (Zhang et al. 2014). Only presence-absence data were collected for vascular plants, precluding use of our modelling framework (but see Simberloff and Gotelli 1984 for random colonization models applied to presence-absence data).

Data analysis
We calculated values of species' abundances, species' probabilities of occurrence and species richness predicted by random placement for each of the 30 study islands using random placement models (1, 2 and 3, respectively). Abundance, occurrence and richness residuals were estimated as observed values minus random placement values. We next used linear mixed effects models to simultaneously: 1) quantify relationships between either abundance residuals or occurrence residuals and island variables, including area, isolation, proportion of suitable habitat and presence/absence of preferred larval host plants; 2) evaluate whether area per se or isolation differentially affects species' abundances or occurrences contingent on their functional traits. Separate models were built for each isolation buffer size (250, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500 and 5000 m). Each of these models had the same number of parameters (k) and only differed in the buffer size used to measure island isolation. We therefore directly compared models with log likelihood (Lamb et al. 2018), where the model with the maximum log likelihood identified the optimal buffer size. Relationships between log likelihood and buffer sizes were quantified using Pearson's product moment correlation coefficients. The proportion of suitable habitat for each species was estimated as the total area of suitable habitat types divided by the total area of the island. Habitat types were classified as suitable for a species if we observed at least one individual within them during the repeated full-island surveys. The presence/absence of preferred larval host plants (compiled from Hall et al. 2014, Acorn andSheldon 2017) for each species on each island was included as a binary variable. We included species' wingspan (mm; as reported in Burke et al. 2011) in models as a functional trait variable, serving as a measure of both body size and mobility/dispersal ability (Roland and Taylor 1997, Lens et al. 2002, Ewers and Didham 2006, Öckinger et al. 2009, MacDonald et al. 2018a. Other functional traits predicted to relate to species' sensitivity to fragmentation (e.g. degree of ecological specialization, trophic position) were either not measured or did not vary substantially among butterfly species and so were not investigated. As an inverse measure of species' rarity, we included species' total abundance and prevalence in abundance and occurrence models, respectively. To evaluate whether area per se or isolation differentially affected species' abundances or occurrences contingent on their functional traits, we included the following interaction terms: wingspan:area, wingspan:isolation, rarity:area and rarity:isolation. All non-binary predictor variables were standardized (subtracting the mean and dividing by standard deviation), permitting comparisons of effect sizes. The structure for both abundance and occurrence linear mixed effects models was as follows, where 'habitat' is the proportion of suitable habitat for each species on each island and 'plants' is the presence/absence of each species' preferred larval host plants: Linear models were fitted for species richness following the same basic protocol as abundance and occurrence linear mixed effects models. Predictor variables included island area, isolation, habitat diversity and vascular plant species richness. The most supported isolation buffer size was again assessed using log likelihood comparisons among models differing only in buffer size. Habitat diversity was estimated as the number of habitat types on each island. All predictor variables were standardized. The structure for richness linear models was as follows, where 'habitat' is the total number of habitat types recorded on each island and 'plants' is vascular plant species richness:

Results
A total of 869 butterflies belonging to 34 species were observed during repeated full island surveys. Butterfly abundance data are reported in Supporting information. One species, Feniseca tarquinius, uniquely feeds on woolly aphids in its larval stage (Hall et al 2014, Acorn andSheldon 2017).
Only one individual of this species was observed across all surveys. We excluded it from abundance and occurrence models, which included presence/absence of preferred larval host plants as a predictor variable. All individuals belonging to all species were included in richness models.

Species' abundances
Comparing log likelihood among linear mixed effects models differing only in isolation buffer size resolved that the proportion of water within 250 m (smallest buffer size) was most supported (Table 1, Fig. 2a). Model support significantly declined across increasing buffer sizes (r = −0.709; p = 0.015), indicating that the amount of habitat immediately surrounding individual islands better predicted species' abundances than the amount of habitat at broader spatial scales. Within the most supported model, abundance residuals were significantly related to both island area and isolation (Table 2, Fig. 3). Area per se had a significant positive effect on species' abundances, while isolation had a significant negative effect, in accordance with mechanisms predicted by the theory of island biogeography. The absolute magnitude of the effects of area per se and isolation, inferred from standardized regression coefficients, were approximately equivalent. Together, these results indicate that the sample-area effect cannot account for variation in species' abundances across islands of varying area and that habitat configuration, and not just total area, has important effects on species' abundances in this system. Other island variables, including the proportion of suitable habitat and presence/absence of preferred larval host plants, were not related to species' abundances. Coefficients of the wingspan:area and rarity:area interaction terms were significantly negative, indicating that effects of area per se systematically varied across species in respect to these functional traits. Causality behind the wingspan:area interaction is clear; effects of area per se on abundance were greater for smaller, less mobile butterfly species. However, for rarity:area, it cannot be resolved whether effects of area per se on abundance were greater for rare species, or whether these species were rare within the dataset because they experience greater effects of area per se. Comparisons of the relative abundances of species between the mainland (continuous habitat) and islands (fragmented habitat) would help resolve the causal direction of this relationship; however, this was beyond the scope of this study. Relationships between abundance residuals and functional trait variables were not significant. This result is expected, as abundance residuals were standardized for each species individually before they were concatenated for use in linear mixed effects models.

Species' occurrences
As with abundance linear mixed effects models, the most supported isolation buffer size for predicting occurrence residuals was 250 m (Table 1, Fig. 2a). Model support significantly declined across increasing buffer sizes (r = −0.746; p = 0.008). Within the most supported model, occurrence residuals showed no relationship to island area, suggesting that the sample-area effect sufficiently accounts for variation in species' occurrences across islands that vary in area (Table 2, Fig. 3). This result confers support for the passive sampling/habitat amount hypotheses. However, occurrence residuals were significantly negatively related to island isolation, suggesting that island configuration has important effects on species' occurrences, as predicted by the theory of island biogeography. Other island variables, including the proportion of suitable habitat and presence/absence of preferred larval host plants, were not significantly related to occurrence residuals.
Effects of functional trait variables (wingspan and rarity) and their interaction with island variables were not significant at α = 0.05. However, the coefficient of the wingspan:area interaction term was marginally significant at p = 0.069, suggesting that smaller, less mobile species were less likely than larger, more mobile species to occur on small islands. In other words, butterfly species richness on small islands may be disproportionately comprised of large butterfly species with high mobility.

Species richness
The most supported isolation buffer size for predicting species richness residuals was 1500 m, which was only marginally more supported than the smallest (250 m) buffer size (Table 1, Fig. 2a). Notwithstanding, model support still significantly declined across increasing buffer sizes (r = −0.833; p = 0.001). Within the most supported model, residuals were not significantly related to island area, isolation, habitat richness or vascular plant species richness (Table 2, Fig. 3). It should be noted, however, that the directionality of the island isolation coefficient aligned with predictions of the theory of island biogeography, and was significant at α = 0.10. Failure to resolve a significant effect at α = 0.05 suggests that richness-based analyses confer less analytical power than abundance-and occurrence-based analyses. If we adopt the most commonly used significance threshold of α = 0.05, our linear model would suggest that the sample-area effect sufficiently accounts for variation in richness observed across our study islands, conferring support for the passive sampling/habitat amount hypotheses and random assembly of species.   Table 2 and Figure 3. Separate models were built using different island isolation buffer sizes, quantifying using the proportion of open water within 250, 500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, and 5000 m of island shores. To permit comparisons of relationships between log likelihood values and buffer sizes among abundance, occurrence, and species richness model sets, log likelihood scores were standardized for each model set by subtracting the mean and dividing by standard deviation. Model support significantly declined across increasing isolation buffer sizes in all instances. b) Log-log island species-area relationship (ISAR) for butterflies occurring on 30 study islands. The shaded blue polygon represents 95% confidence intervals for the loglog ISAR linear regression (solid blue line), parameterized using observed richness values for all 30 islands. Random placement richness values were calculated using a random placement model (model 3; see Materials and Methods). The shaded green polygon represents 95% interpolated confidence intervals for random placement values, calculated using Coleman's (1981) formula for variance. Table 2. Standardized regression coefficient estimates ('coeff.'), standard errors ('SE') and p-values ('p') from linear models fitting random placement residuals for species' abundances, species' occurrences and species richness. Included in all models were island area, measured in m 2 and island isolation, measured as the proportion of water within the most supported buffer size (250 m for abundance and occurrence; 1500 m for species richness). Within abundance and occurrence models, the proportion suitable habitat ('habitat') was measured for each species as the area of suitable habitat on each island divided by the area of the island. Presence/absence of each species' preferred larval host plants ('plants') was included as a binary variable. Wingspan was included as a measure of species' body size and as a proxy of dispersal ability. Each species' total abundance and prevalence were used as an inverse measure of rarity in abundance and occurrence models, respectively. Species identity was included as a random effect in abundance and occurrence models. Within the species richness model, habitat diversity ('habitat') was estimated as the total number of habitat types recorded on each island.

Discussion
Distinguishing mechanisms that underlie variation in species' abundances, species' occurrences and emergent patterns of species richness is not only of great interest within the context of fundamental ecology, but is also of paramount importance to applied ecology and understanding how habitat fragmentation affects species diversity (Diamond 1975, Simberloff and Abele 1976, 1982, Haila 2002, Haddad et al. 2015, Chase et al. 2019). Here, we detail a novel modelling framework to test hypotheses on which conservation directives are contingent (Fahrig 2003, Haddad et al. 2015, Hanski 2015, Fletcher et al. 2018. Applying this modelling framework to a butterfly assemblage persisting on a naturally fragmented landscape of true islands, we were able to resolve that: 1) island area per se and isolation significantly affect species' abundances and occurrences contingent on their functional traits; and 2) important effects of area per se and isolation are not always apparent in aggregate diversity measures, such as those reflected in ISARs. Although there are several documented divergences between the biogeographies of true-island systems and fragmented habitat on terrestrial landscapes, findings from our study clearly demonstrate that fragmentation effects should not be inferred from richness-based analyses, but rather evaluated on a species-by-species basis.

Inferences from the ISAR
Our modelling framework resolved that spatial patterns in butterfly species richness (i.e. the ISAR) did not significantly deviate from random placement in relation to island area, isolation, habitat diversity or vascular plant diversity. These ISAR-based inferences align with those of MacDonald et al. Figure 3. Standardized regression coefficients and 95% confidence intervals from linear models relating random placement residuals to island characteristics and species' functional traits for (a) species' abundances, (b) species' occurrences and (c) species richness. Included in all models were island area, measured in m 2 and island isolation, measured as the proportion of water within the most supported buffer size (250 m for abundance and occurrence; 1500 m for species richness). Within abundance and occurrence models, the proportion suitable habitat ('suitable habitat') was measured for each species as the area of suitable habitat on each island divided by the area of the island. Presence/absence of each species' preferred larval host plants was included as a binary variable. Wingspan was included as a measure of species' body size and as a proxy of dispersal ability. Each species' total abundance and prevalence were used as an inverse measure of rarity in abundance and occurrence models, respectively. Species identity was included as a random effect in abundance and occurrence models.
Within the species richness model, habitat diversity was estimated as the total number of habitat types recorded on each island. Plant diversity was measured as vascular plant species richness. The shading of each variable's point estimate (coefficient) and confidence interval is proportional to its p-value, with darker shades indicating greater significance. Coefficients with 95% confidence intervals not overlapping zero were inferred to be significant at α = 0.05.
(2018a), who used a series of SLOSS-based analyses to infer that butterfly species richness in this naturally fragmented landscape approximately conform to predictions of the sample-area effect, with no significant effects of area per se or isolation. Therefore, all analyses addressing effects of area per se and isolation on butterfly species richness in this system support the passive sampling/habitat amount hypotheses and Fahrig's (2003Fahrig's ( , 2013Fahrig's ( , 2017 general conclusion that effects of fragmentation (i.e. area per se and isolation) are generally negligible after controlling for deleterious effects of habitat loss; only habitat amount at the landscape-scale affects species richness via its influence on regional species pools.
While there seems to be considerable support for the passive sampling/habitat amount hypotheses in the literature (Fahrig 2003, Martin 2018, Watling et al. 2020), a recent meta-analysis by Fahrig (2020), including 157 SLOSS comparisons from 58 studies, found that several small islands/ fragments contained more species than single large islands/ fragments in 72% of comparisons, equivalent numbers of species in 22% of comparisons and fewer species in 6% of comparisons. Removing studies with biased sampling effort in relation to island/fragment area shifted these figures to 58, 37 and 5%. Regardless, these results suggest that species richness varies with degree of fragmentation more often than it does not. Thus, the sample-area effect implicated in the passive sampling/habitat amount hypotheses cannot consistently account for SLOSS-based richness patterns. Furthermore, methods employed by Fahrig (2020) -specifically, comparisons of Quinn and Harrison (1988) species accumulation curves -suffer from an important limitation: substantial species turnover among several small islands/fragments can inflate their aggregate richness, such that important deviations from random placement (e.g. effects of area per se) are obscured (sensu Simberloff 1976, MacDonald et al. 2018b, Deane et al. 2020. This relationship may explain why several small islands/fragments are generally found to contain more species than single large fragments in the majority of SLOSSbased studies. While results of Fahrig's (2020) meta-analysis are of great interest to both fundamental and applied ecology, they cannot necessarily be used to distinguish effects of area per se from the sample-area effect and cannot resolve additional effects of isolation or habitat diversity. Future investigations focusing on species richness would benefit from the inclusion of ISAR-based analyses that assess deviations from random placement on an island-by-island or fragment-byfragment basis, such as those included within the modelling framework presented here.
Additional deviations from random placement may be resolved by visually examining the ISAR and random placement richness values (Fig. 2b). In this study, observed richness values were generally less than random placement richness values (all islands except four). This pattern cannot be attributed to effects of area per se, because the direction and magnitude of richness residuals was relatively consistent across islands of varying area (Fig. 2b), as indicated by the insignificant coefficient of island area within the species richness linear model (Table 2). Rather, this relationship is best explained by spatial species aggregation, wherein conspecific individuals are more likely to co-occur on islands than what is predicted by random placement, reducing the species richness of individual islands (He and Legendre 2002). Therefore, although spatial patterns of species richness did not significantly deviate from random placement in respect to either island area or isolation, they did deviate from random placement in respect to spatial species aggregation. It is therefore clear that failure to resolve significant effects of area per se and isolation on species richness cannot be taken as direct evidence for random assembly of individuals and species, the fundamental prediction of the habitat amount/passive sampling hypotheses (c.f. Fahrig 2003Fahrig , 2013Fahrig , 2017Fahrig , 2020. Rigorous evaluation of random assembly requires comparison of observed richness to expected richness predicted by null models, such as the random placement models presented here.

Inferences from species' abundances and occurrences
Inferring whether fragmentation is 'good' or 'bad' (sensu Fahrig 2017, Fletcher et al. 2018) based on emergent patterns of species richness is potentially susceptible to 'ecological fallacy' (sensu Robinson 1950), which describes biases that may arise when observed effects on aggregated variables (e.g. species richness) differ from causal relationships at more reductive and informative levels of organization (e.g. species' abundances and occurrences). Indeed, our abundance and occurrence models resolved important effects of area per se and isolation that were not apparent in either ISARor SLOSS-based analyses (this study and MacDonald et al. 2018a, respectively). This discrepancy among inferences suggests that conflating responses of all species into a single aggregate measure (e.g. species richness) reduces our power to detect important relationships on which conservation directives should be contingent. Two such relationships were resolved when considering the entire species assemblage in abundance and occurrence models: 1) there was a disproportionate concentration of individuals on larger and less-isolated islands relative to what was predicted by the sample-area effect (passive sampling/habitat amount hypotheses); and 2) species were more likely to occur on less-isolated islands than what was predicted by the sample-area effect. It is therefore clear that the sample-area effect described by the passive sampling/habitat amount hypotheses cannot adequately account for spatial patterns in butterfly abundances and occurrences in this naturally fragmented landscape, which are better predicted by mechanisms outlined by the theory of island biogeography. It should be noted that the directionality of area per se and isolation effects were consistent among abundance, occurrence and richness models, but only statistically significant (α = 0.05) in the first two analyses, wherein species' responses were not aggregated into a single measure.
Most interestingly, our modelling framework simultaneously resolved that effects of area per se on species' abundances and occurrences varied significantly with species-specific functional traits, suggesting that mechanisms outlined by the theory of island biogeography are not neutral with respect to species identity. Whether species' sensitivity to fragmentation varies predictably with functional traits is a long-standing and pertinent question in conservation biology (Roland and Taylor 1997, Lens et al. 2002, Gehring and Swihart 2003, Henle et al. 2004, Tscharntke and Brandl 2004, Thies et al. 2005, Ewers and Didham 2006, Prugh et al. 2008, Barbaro and Van Halder 2009, Öckinger et al. 2009, Hanski 2015, Warzecha et al. 2016, MacDonald et al. 2018a. In this study, effects of area per se were significantly greater for butterfly species with smaller wingspans. For Canadian butterflies, wingspan is one of the strongest correlates of estimated mobility and dispersal ability (Burke et al. 2011); thus, we can infer that effects of area per se are greatest for small species of limited mobility (c.f. Larsen et al. 2008). This relationship may be explained by larger and more mobile butterfly species having the ability to move among multiple small islands to meet their resource requirements. These 'transient' species may thereby exhibit patterns of abundance and occurrence that approximate random placement (Wilson and MacArthur 1967: Chapter 2;Rosenzweig 2004, MacDonald et al. 2018a. Such patterns would be predicted by ideal free distribution theory (sensu Dreisig 1995) if two conditions are met: 1) islands of varying area contain equivalent densities of resources; and 2) the mobility/dispersal ability of individuals is sufficient to render costs of inter-island movements negligible. By contrast, island edges may be perceived as impassible barriers for smaller and less mobile species, with energy expenditures and mortality risks associated with movement through the open-water matrix being too high for regular inter-island movements. Island edges may therefore delimit populations of smaller and less mobile species, which are generally restricted to larger islands that contain all resources required for mate location, reproduction, resting, roosting, predator escape and feeding (i.e. the functional resource-based habitat concept; Dennis et al. 2003). This hypothesis is supported by analyses of MacDonald et al. (2018a), who resolved that the probability of butterfly species occurring on islands without their preferred larval host plants was positively related to their wingspan and estimated mobility. Considered together, these results suggest that functional traits may be used to predict species' sensitivity to fragmentation and that species identity should not be ignored when investigating mechanisms that underlie ISARs or in conservation planning. It is, however, important to recognize that the open-water matrix of this study landscape may be more unsuitable and less permeable than those of many fragmented terrestrial landscapes (Dunning et al. 1992, Ricketts 2001, Laurance 2008, Mendenhall et al. 2014, Itescu 2019, Farneda et al. 2020. It is therefore unclear the degree to which these relationships between species' functional traits and effects of area per se and isolation are generalizable to conservation efforts addressing terrestrial landscapes fragmented through anthropogenic activities. The abundance and occurrence modelling framework proposed here may also be implemented on a species-by-species basis by regressing island/fragment variables (area, isolation, presence/amount of suitable habitat or specific resources, etc.) on abundance and occurrence residuals for each species in separate linear models. This method of analysis precludes the simultaneous integration of functional trait analyses within models, but has the added advantage of identifying single species that are particularly sensitive to fragmentation (area per se and isolation) or other island/fragment variables of interest. This simple decomposition of our modelling framework may be used to resolve whether particular species require independent consideration within conservation frameworks.

Isolation versus habitat amount
The relative isolation of islands addressed in this study was quantified across multiple scales as the proportion of water within 11 buffer sizes, ranging from 250 to 5000 m. These measures are equal to 1 minus the amount of landmass (habitat) within each buffer distance. Fahrig (2013) suggests that the habitat amount hypothesis would be supported by species' abundances, species' occurrences or species richness of equal-area sampling plots (stratified across fragments of varying area and isolation) correlating with the amount of habitat on the surrounding landscape more strongly than with the area of the individual fragments on which the sampling plots are located. This is because landscapes containing less habitat should contain fewer individuals (belonging to fewer species) due to the sample-area effect, meaning fragments within such landscapes will have smaller species pools from which their own diversities are randomly sampled. However, the 'appropriate distance' for quantifying the amount of habitat surrounding equal-area sampling plots is undefined and, most problematically, the area of individual fragments on which sampling plots are located becomes increasingly correlated with habitat amount as this distance is reduced. Thus, it may not be possible to decouple fragment area from habitat amount using Fahrig's (2013) proposed method; particularly, for taxa that respond to habitat amount and configuration at fine spatial scales, including butterflies (Thomas and Abery 1995, MacDonald et al. 2017, 2018a, Saura 2020. Within our modelling framework, variation in species' abundances, species' occurrences and species richness associated with full-island surveys and the sample-area effect is nullified in the calculation of random placement residuals. Abundance, occurrence and richness residuals can therefore be correlated with island/fragment area and the amount of surrounding habitat in a fashion similar to Fahrig's (2013) proposed method of using equal-area sampling plots. However, because the area of individual islands/fragments is not included in our measures of the amount of surrounding habitat (our buffers are generated from island/fragment edges), the problem of island/fragment area becoming increasingly correlated with habitat amount at fine spatial scales is avoided. After controlling for effects of area per se, abundance and occurrence residuals were significantly related to island isolation. Model support significantly declined across increasing isolation buffer sizes for species' abundances, species' occurrences and species richness, suggesting that that the amount of habitat immediately surrounding islands, rather than the amount of habitat at broader landscape scales, has the greatest effect on butterflies in this system. Importantly, island area and isolation were effectively decoupled, as indicated by the absence of correlation between island area and isolation (e.g. 250 m buffer; r = 0.082; p = 0.668). Thus, in contrast with mechanisms outlined by the habitat amount hypothesis, isolation effects observed in this study are better attributed to the fact that individual butterflies moving through the open-water matrix are less likely to encounter more isolated islands (Andrén 1994), reducing species' abundances and probabilities of occurrence, as predicted by the theory of island biogeography. This inference is further corroborated by analyses of MacDonald et al. (2018a), which showed that butterfly species turnover among islands of equal areas -a proxy for variation in species pools if islands indeed randomly sample species -was unrelated to Euclidean distance between islands. This result suggests a uniform species pool throughout the study landscape. Again, the degree to which these findings apply to fragmented terrestrial landscapes, wherein the suitability and permeability of matrices may vary, is unclear. For studies addressing fragmented terrestrial landscapes, isolation measures should not only account for the proportion of suitable habitat within various buffer distances, but also include measures of matrix suitability and permeability, if they are available (MacDonald et al. 2020).

Habitat fragmentation and the narcissus effect
It is important to recognize a bias within this case study, and potentially other study designs addressing spatial patterns of species' abundances, species' occurrences or species richness across true islands or terrestrial habitat fragments. Here, we investigated effects of area per se and isolation on a naturally fragmented landscape of true islands with substantial time-since-isolation (3000-4000 YA; Yang and Teller 2005). Therefore, species that are particularly sensitive to fragmentation are unlikely to occur on islands at all. Although our modelling framework resolved significant effects of area per se and isolation on butterfly species' abundances and occurrences, effects of area per se and isolation on the regional species pool were likely underestimated, as the species assemblage of adjacent continuous habitat was not quantified. This bias may be described as the 'narcissus effect', which addresses situations wherein a null model or study design unintentionally accounts for or excludes effects that are of interest (sensu Colwell and Winkler 1984). It is possible that this bias contributed to results of Fahrig's (2017) review, where 68% (158/232) of studies addressing single species reported positive fragmentation effects; species that are particularly sensitive to fragmentation may be completely missed in many studies. We therefore suggest caution in interpreting results from study designs that are susceptible to the narcissus effect and encourage future studies to compare the identities and functional traits of species between islands/fragments and adjacent continuous habitat to assess potential biases resulting from the historic exclusion of fragmentation-sensitive species. This may be accomplished using our proposed modelling framework by surveying continuous habitat equal in area to the sum of all surveyed islands/fragments. Abundance data from continuous habitat may then be used in place of total abundances across all islands/fragments (n i ) to calculate abundance, occurrence and richness random placement values for each island/fragment using the random placement models described above. Subtracting these random placement values from observed abundance, occurrence and richness values for each island/fragment will result in abundance, occurrence and richness residuals that may be used in our linear mixed effects models (abundance and occurrence) and linear models (richness) to resolve whether there are additional effects of area per se and isolation on the regional species pool.

Conservation implications
Considerable uncertainty exists in the literature regarding the influence of area per se and isolation (i.e. habitat fragmentation) on populations and communities of wildlife (Fahrig 2003, Haddad et al. 2015, Hanski 2015, Fletcher et al. 2018. There is an immediate need to resolve this debate, as habitat loss and fragmentation are widespread and increasing (Hanski et al. 2013, Ibisch et al. 2016, Chase et al. 2020, Deane et al. 2020. We demonstrate here that ISAR-and SLOSS-based inferences, founded on emergent patterns of species richness, have the potential to obscure important interspecific variation in responses to area per se and isolation. To infer support for the passive sampling, habitat amount and related hypotheses from emergent patterns of species richness that spuriously conform to predictions of the sample-area effect is to simplistically cut rather than carefully untie the Gordian knot of ecological complexity. We suggest that, in addition to emergent patterns of species richness, information at more reductive and informative levels of organization (e.g. species' abundances and occurrences) should be included in studies aiming to measure and understand effects of habitat fragmentation.
Acknowledgements -Vascular plant surveys and identifications were completed by Iraleigh Anderson and we thank him for valuable discussions and insights. We also thank Federico Riva for his review of a previous version of this manuscript. Permits -Permission to survey butterflies and vascular plants, including collection of voucher specimens, was granted by Ontario Parks.