Supporting Information Appendix for

The Area-Heterogeneity Trade-off and the Diversity of Ecological Communities

Omri Allouchea*, Michael Kalyuzhnya, Gregorio Moreno-Ruedab, Manuel Pizarrob, Ronen Kadmona

a

Department of Ecology, Evolution & Behavior, Institute of Life Sciences, The Hebrew University of Jerusalem, 91904, Jerusalem, Israel.

b

Departamento de Zoología, Facultad de Ciencias, Universidad de Granada, 18071, Granada, Spain.

* To whom correspondence should be addressed. E-mail: [email protected]

SUPPORTING INFORMATION APPENDIX METHODS A. Empirical analyses A.1

General considerations

A.2

Study system

A.3

Quantification of the response variables

A.4

Data analysis A.4.1 Community responses to elevation range A.4.2 Effects of niche width

B. Theoretical analyses B.1

Objectives

B.2

Description of the model

B.3

Simulation experiments

C. Meta-analysis C.1

Objectives

C.2

Data source

C.3

Statistical analysis

C.4

Criteria used for functional relationship determination

D. Additional References

A. Empirical analysis A.1 General considerations A comprehensive test of the area-heterogeneity trade-off hypothesis requires data on species richness, extinction rates, mean species abundance, and an appropriate measure of environmental heterogeneity. While data on species richness, species abundance, and environmental heterogeneity can be obtained for a single point in time, determining extinction rates requires data on species composition in a fixed set of sites for (at least) two distinct survey periods. The interval between the two surveys should be larger than the typical life span of the relevant species to allow detectable levels of extinction, but not too large to minimize opportunities for recolonization following extinction events. The spatial extent of the observation unit is also important – it should be small enough to capture the footprints of local stochastic extinctions, but large enough to minimize the impact of immigration and emigration on the results and to obtain sufficient spatial variation in environmental conditions. Choosing an appropriate measure of environmental heterogeneity is also not trivial. Ideally, environmental heterogeneity should be defined with respect to an environmental factor that is known to affect the distribution of the target species in the relevant system. Taking into account that the strength of the areaheterogeneity trade-off depends on niche width (Fig. 1), it should also be possible to quantify niche width of the target species with respect to the selected factor, and there should be evidence that the study species differ in both the position and the width of their niches along the selected environmental gradient.

A.2 Study system Unfortunately, data satisfying all of the above requirements are extremely rare. In this study we used data on breeding bird distributions in Catalonia (NE Spain) as a basis for our analysis. This system was ideal for our purposes because (1) data were collected using the same set of observation units (UTM grid cells of 10x10 km) in two distinct surveys, the first from 1975-1982 with most data collected in 1980-82 (17) and the second from 19992002 (18); (2) the interval between the two surveys (two decades) was larger than the typical lifespan of most bird species (34); (3) the size of the observation unit (10x10km) was small enough to detect extinction events (35) but large enough to show a wide range of within-cell heterogeneity in environmental conditions; and (4) the data provided by the second survey included estimates of species abundance that were generated using logistic regression models based on a stratified random sampling of smaller units (1x1km) within the 10x10km cells. These estimates were available for species that were observed in at least 10 squares and for which cross validation revealed highly accurate predictions. The second atlas (18) also provided estimates of the effective number of hours invested in bird surveying within individual grid cells in each survey period (five categories: 160h), thus allowing to evaluate the potential effect of sampling bias on the results. Analyses of the data collected in the second survey has shown that elevation is a major determinant of breeding bird distributions in this area and that different species show considerable differences in both the position and the width of their niches along the elevation gradient (18). We therefore used elevation range as an effective measure of environmental heterogeneity in our analysis. A total of 372 grid cells for which appropriate data were available were selected for the analysis. These grid cells covered 99% of the terrestrial area of Catalonia.

Since many grid cells (179 out of 372) were border cells or coastal cells in which only part of the grid cell was surveyed, we calculated the terrestrial area of Catalonia in each grid cell by intersecting the UTM grid used to design the two field surveys with a vector map of Catalonia (downloaded from ArcGIS online world administrative divisions, http://www.esri.com/data/free-data/index.html). The 372 polygons were then intersected with a ~1 km resolution DTM of Catalonia (downloaded from http://dds.cr.usgs.gov/srtm/version1/Eurasia/) and the mean elevation of each polygon was determined. All GIS operations were performed using ESRI ArcGIS Desktop 10 (Environmental Systems Research Institute, Redlands, California, USA). Additional data on mean annual rainfall, mean annual temperature, and human population density of each grid cell were available from a previous analysis (36). Determining the appropriate scale for quantifying elevation range was more complicated due to differences in area among the peripheral polygons and the fact that the center of such polygons did not coincide with the center of the corresponding UTM cells. We therefore determined elevation range at three alternative scales using buffer zones with radii of 7km, 10km, and 14km around the center of each polygon. The three measures were highly correlated (Pearson correlation = 0.93-0.97) and had similar effects on species richness. We therefore used the smallest scale (7km) which is approximately one-half of the diagonal of a 10x10km cell, as a basis for quantifying an area-independent measure of elevation range. This radius ensured that all parts of all grid cells were included within the buffer. It also fits the finding of Sutherland et al. (37) that median distances of natal dispersal in birds range from 0.03 km to 10 km.

A.3 Quantification of the response variables Our predictions focused on three community-level characteristics: mean species abundance (Fig. 2A), extinction rates (Fig. 2B), and species richness (Fig. 2C). Species richness was determined for each grid cell as the average number of species that were present in that cell in the two survey periods. Estimates of species abundance were averaged for individual species in each UTM cell and were standardized prior to the analysis in order to control for differences in absolute abundance among species. Two standardization methods were applied, by maximum value and by the mean and standard deviation. In the first method, mean species abundance in a given grid cell was estimated as:

Pj

Pi , j

1 Nj

P i 1

max i

where Pj is the mean (standardized) abundance of all species present in grid cell j, Nj is the number of species for which estimates of abundance were available in grid cell j, Pi,j is the mean abundance of species i in grid cell j (i = 1,...,N), and Pmax i is the abundance of species i in the UTM cell with its highest abundance in Catalonia. In the second method mean species abundance in each cell was estimated as:

Pj

1 Nj

i

Pi , j i

i

where i is the mean abundance of species i in the study area and i is the corresponding standard deviation. We report the results obtained for the first (maximum-based) standardization because individual species abundances were not normally distributed and our analyses showed that both methods gave similar results. Yet, it should be noted that estimates of species abundance were based on predictive models of site occupancy (GLM regressions predicting probabilities of occurrence at a spatial resolution of 1km, 18) rather than direct observations. Thus, although occupancy is usually correlated with abundance

(38) and only species with accurate predictions were included in our analysis, these results should be interpreted with some cautious. Extinction rates were estimated for each grid cell by counting the number of species that were present in that cell in the first survey but not in the second survey (note that this definition includes local 'extinctions' caused by emigration to other cells). Although survey data like those analyzed in this study are known to suffer from imperfect detection, we believe that our estimates of extinction rates were not very sensitive to such bias because of the intensive sampling effort of the second survey (18). Spatial differences in sampling efforts might cause bias in breeding bird surveys (39), but this source of bias was minimized in the second survey by using a standardized sampling protocol, originally designed to minimize such bias. Moreover, a detailed analysis of the overall data collected from 1983-2002 showed that sampling effort in this period was positively and significantly correlated with elevation range (39). All expected consequences of such bias (higher species abundances, lower extinction rates and higher species richness at grid cells characterized by high elevation range) are exactly opposite to our predictions. We therefore believe that patterns of abundance and species richness supporting our predictions cannot be attributed to underlying bias in sampling effort. This assumption was confirmed by complementary analyses (see below).

A.4 Data analysis A.4.1 Community responses to elevation range The effect of elevation range on mean species abundance, extinction rates, and species richness was analyzed using four types of regression models. The first set of models tested the combined effects of elevation range (H) and its squared term (H2) on each of the three response variables using ordinary least squares (OLS) regressions. The squared value of

elevation range was included in the models to account for possible unimodal effects (though we do not expect any relationship to be explicit parabolic). These and all other OLS models were performed using SPSS for Windows, Ver. 15.0 (SPSS, Chicago, IL, U.S.A). The results of these analyses are presented in the first column of SI Appendix Tables 1-3. In the analysis of mean species abundance (SI Appendix Table 1) the quadratic term of elevation range was negative and highly significant but the linear term was not significant, indicating a non-linear negative response (Fig. 2D). In the analysis of extinction rates (SI Appendix Table 2) the linear term of elevation range was positive and statistically significant but the squared term was not significant, indicating a positive linear response (Fig. 2E). In the analysis of species richness (SI Appendix Table 3) the linear term was positive, the quadratic term was negative, both terms were statistically significant, and the inflection point was within the range of the data, indicating a unimodal response (Fig. 2F). These results are consistent with our theoretical predictions (Fig. 2AC). To evaluate the robustness of these patterns to other, potentially confounding factors that might influence the dynamics and structure of the study species, we constructed more complex OLS models in which we incorporated the effects of area (A), mean absolute elevation (E), mean annual rainfall (R), mean annual temperature (T), human population density (P), and broad-scale spatial effects (X coordinate, Y coordinate, X2, Y2, and XY) as additional predictors in the models. Analyses were performed using the backward variable selection method of SPSS with significance levels to enter and to leave the model being set at α = 0.01 and 0.05, respectively. The results were similar to those obtained in the simpler models, indicating a negative effect of elevation range on mean species abundance, positive effect on extinction rates, and unimodal effect on species

richness (SI Appendix Tables 1-3). Similar results were obtained in models incorporating the simultaneous effects of all 12 independent variables and in the 'best' models selected from all possible models using the AIC criterion. We also repeated the analysis replacing area with the effective surface area, which corrects for area rugosity, and obtained similar results, as expected due to the small variation of area introduced by considering rugosity (less than 5%). Although broad-scale geographic trends in the response variables were partially accounted for by the OLS models (SI Appendix Tables 1-3), the residuals of all models were positively autocorrelated (Moran's I < 0.05) at spatial scales lower than 40km (40). We therefore performed two additional sets of analyses in which all variables that were found statistically significant in the OLS models were reanalyzed for their effects on the response variables using Conditional Autoregressive Regression (CAR) and Simultaneous Autoregressive Regression (SAR) models (41). Both types of models assume that errors of the respective non-spatial models are spatially autocorrelated. As spatial autocorrelation was limited to small spatial scales (40km), the spatial models were performed using binary connectivity matrices with a distance threshold of 40km (42). All spatial analyses were performed with SAM Ver. 4.0 (Spatial Analysis in Macroecology (43)). The results obtained for the spatial models were qualitatively similar to those obtained for the non-spatial models (SI Appendix Tables 1-3). A comparison of the four regression models constructed for each response variable based on the coefficient of determination (R2) and Akaike's Information Criteria (44) adjusted for small sample size (AICC) indicated that the SAR models had the best fit to the data (highest R2 and lowest AICC) in all analyses (SI Appendix Tables 1-3).

Two additional analyses were performed to verify that the observed responses to elevation range were not artifacts of spatial or temporal variation in sampling effort. First, we repeated the original analyses using only grid cells that had a similar sampling effort for both surveys. This procedure reduced sample size by 50% (from 372 to 189) but did not affect the results of the original analyses (SI Appendix Fig. 1, SI Appendix Table 4). Second, the combined effects of elevation range (both the linear and quadratic terms) and sampling effort (160h) on species richness were tested separately for each survey period using GLM. The results indicated that sampling effort had a significantly positive effect on species richness but for both survey periods the effect of elevation range was similar to that in the original analysis (SI Appendix Fig. 2, SI Appendix Tables 5-6). We therefore feel quite confident that the unimodal pattern of species richness was not a sampling artifact.

A.4.2 Effects of niche width Further analyses were performed to test the prediction that increasing niche width shifts the level of elevation range that maximizes species richness to higher levels. In addition to its theoretical significance, this analysis was methodologically important for two reasons. First, it has often been suggested that unimodal responses of species richness to environmental factors in general and elevation in particular can arise through a simple middomain effect (45). In our system elevation range was correlated with mean elevation (r = 0.75, P < 0.001). It was therefore important to reject the possibility that the observed effect of elevation range on species richness is not a statistical artifact of this correlation. According to mid-domain theory, differences among species in niche width should affect the height of the response to the relevant gradient, but not the position of the inflection point (45). A statistically significant effect of niche width on the position of the inflection

point can therefore be interpreted as evidence against a mid-domain explanation of the unimodal pattern observed at the level of the whole community. Second, although the unimodal effect of elevation range on species richness was robust to variation in a number of potentially confounding factors (SI Appendix Tables 13), it is possible that the observed response was caused by other confounding factors for which we did not have data. A consistent and directional shift in the position of the inflection point with increasing niche width is not expected to arise from other simple mechanisms and can therefore be interpreted as a further support for the area-heterogeneity trade-off. Niche width of the study species was determined by ranking the records of each species with respect to elevation, removing outliers (highest and lowest 5% of the records), and calculating the difference between the remaining maximum and minimum elevation values. Based on the calculated niche widths we categorized the species into four 'niche quartiles' and regressed the number of species belonging to each niche quartile against elevation range and its squared value.

B. Theoretical analyses B.1 Objectives The theoretical part of the study had three main objectives. The first objective was to confirm the logic and validity of our predictions (Fig. 2A-C) by deriving the expected relationships between environmental heterogeneity, mean species abundance, extinction rates, and species richness, from the fundamental processes of birth, death, and migration of individuals. A second objective was to get a sense of the 'noise' caused by pure demographic stochasticity. Such an evaluation was important because the empirically observed data were very noisy (Fig. 2D-F). The third objective was to evaluate whether and

how community responses to environmental heterogeneity depend on the demographic rates of the modeled species, their niche width, the area, and the length of time over which the data are collected. Answering all of these questions is crucial for a proper interpretation of empirically observed heterogeneity-diversity relationships. Evaluating the magnitude of 'noise' expected in empirical tests of our predictions is also important because ecologists studying patterns of species diversity tend to ignore the inherent noise expected from demographic stochasticity and often attribute low values of R2 to the effect of ‘hidden’ variables not included in the analysis.

B.2 Description of the model Our theoretical analyses were based on a stochastic, individual based, spatially implicit model of a local community that receives immigrants from a regional species pool with N species. Individuals go through the fundamental demographic processes of reproduction, mortality and migration, all modeled with explicit rates, that are constant among species. The landscape consists of A sites, where each site can be occupied by at most one individual. Each site is conceptualized as having a unit area, so that A also indicates the total area available for the local community. Environmental heterogeneity is introduced by assigning each site an environmental value Ej that is randomly drawn from a uniform distribution between Emin and Emax. The range of E available for the local community (ER = Emax-Emin) represents the degree of environmental heterogeneity. This range is embedded within a larger range (ERR = ERmin to ERmax) representing the range of E in the area available for the regional species pool (SI Appendix Fig. 4). In this formulation, the environment is reduced into a single environmental variable representing a major niche axis (19,46) or a complex gradient such as elevation (47).

The niche of each species along the environmental gradient is described by a Gaussian function indicating the probability of establishment as a function of E (SI Appendix Fig. 4): ( E j μi )2

Pri , j Z i e

2 σi 2

,

where Pri,j(E) is the per-capita probability of establishment of species i in site j ( i {1,..., N}, j {1,..., A} ), µi is the optimum level of E for species i, and σi is a parameter representing the corresponding niche width (assumed to be equal for all species in most analyses). Values of µi are randomly selected from a uniform distribution between ERmin and ERmax. Zj is a normalization constant ensuring that all species are equivalent in their overall probability of establishment within the range ERmin – ERmax: E R max

Zj ( R E

[

e

(E j μi )2 2 σi 2

]

dE)1

min

Similar results were obtained for equal distribution of environmental conditions, created by dividing the local environmental range ERR into A equally-spaced values and assigning each to one local site. Each individual dies at an annual rate M and gives birth to one offspring at rate R. A new offspring is immediately dispersed into a random site. Immigration from the regional species pool to each site in the landscape occurs at a rate I, with equal immigration probability for each of the N species in the regional species pool. Individuals arriving at empty sites establish with probability Prij and individuals arriving in occupied sites die without being able to establish. The dynamics of the local community is modeled as a continuous-time Markov process where, in a small enough time interval, only a single demographic event (birth, death or immigration) can occur (48). This model can be considered as an individual-based, stochastic extension of Levins' model (49) in which all species are identical except for the position of their niche along an

environmental gradient. It can also be considered as an extension of the model of Kadmon and Allouche (14). While in Kadmon and Allouche's model habitats were discrete and unordered and each species could only persist in a single habitat, in the current model the environmental axis is continuous and species can persist in multiple environmental conditions with varying success. We simulated the system using Poisson thinning (50,51), implemented in C and MATLABTM (MATLAB version 7.10.0.499, The MathWorks Inc., Natick, MA)

B.3 The simulation experiments The effect of environmental heterogeneity on community characteristics was evaluated by varying the range of environmental conditions available for the local community (ER = Emax-Emin) while keeping all other parameters constant. Two sets of simulations were conducted. The first set was designed to enable evaluation of the 'noise' expected in empirical tests of our predictions due to demographic stochasticity. In these simulations each response variable was sampled only once in each simulation run (in order to ensure independence) and no averaging was applied to the results. The second set was designed to evaluate the robustness of the predicted patterns to variations in model parameters. In these simulations the results were averaged within and among different runs in order to remove the noise caused by demographic stochasticity from the observed patterns (SI Appendix Figs. 5-8). Preliminary analyses showed that all variables reach their steady state distribution after less than 5,000 years. We therefore started all simulations by running the model for 20,000 years with initial conditions of zero local abundance of all species. Similar results were obtained for simulations of a shorter time period of 50 years, with initial conditions of zero local abundance of all species, in which the community did not reach it steady state

distribution (SI Appendix Fig. 9). Sampling started from year 20,000 and differed between the two sets of simulations. In the first set (designed to generate independent realizations) the modeled community was sampled using only two snapshots, at the end of year 20,000 and at the end of year 20,001. The first snapshot was used to determine mean species abundance (average number of individuals per species) and species richness (number of species) in the local community. The second snapshot was used to determine extinction rates (number of species extinct between year 20,000 and 20,001). This procedure was repeated using 20 independent simulations for each of 15 levels of environmental range (3, 6, 9,…,45) under a fixed parameter regime (N = 500, R = 45, M = 1, I = 0.2, σ = 2, A = 1,000, ERR = 65). The qualitative results of these simulations (Fig. 2, G to I) were consistent with our empirical findings, indicating that environmental heterogeneity reduces mean species abundance, increases extinction rates, and has a unimodal effect on species richness. It is also important to note that all community responses were extremely ‘noisy’, suggesting that empirical analyses of these patterns are expected to show very low values of R2 even at the absence of any other deterministic or stochastic source of environmental variation. The only difference between the simulations and the empirical data was the functional form of the abundance-heterogeneity relationship (Fig. 2E,H) which may reflect restrictive assumptions of our model (e.g., ignorance of environmental stochasticity) or differences in the manner by which abundance was determined in the field and in our simulations. In the second set of simulations (designed to evaluate the robustness of the predicted patterns to variations in model parameters) each run lasted 120,000 years and snapshots of the community were sampled at the end of each year starting from year 20,000. The 100,001 snapshots obtained for a given run were used to calculate mean values of the three response variables (mean species abundance, extinction rates, and

species richness). This procedure was repeated using 10 independent runs for each level of environmental range under each parameter regime. We used the second type of simulations to test the manner by which the effects of environmental heterogeneity depend on reproduction rates, immigration rates, area, and temporal scale of data collection. The effect of each of these factors was tested under a wide range of values (reproduction: 30120; immigration: 0.01-0.5; area: 100-50,000, temporal scale of data collection: 1-20 years). Further simulations were performed to test the hypothesis that increasing niche width shifts the level of heterogeneity that maximizes species richness to higher levels. In these simulations the regional species pool consisted of a mixture of species representing different niche widths and the number of independent runs was 50 due to higher stochasticity. Typical results of our simulations are presented in SI Appendix Figs. 5-8. As expected from the area-heterogeneity trade-off, the general relationship between species richness and environmental heterogeneity was unimodal, but the level of heterogeneity that maximized species richness increased with increasing reproduction, immigration, area, temporal scale of data collection, and niche width. All of these factors reduce the likelihood of stochastic extinction, thereby increasing the range of environmental heterogeneity under which the relationship is positive.

C. Meta-analysis C.1 Objectives The purpose of this part of the study was to test our hypothesis that evidence of positive heterogeneity-diversity relationships obtained for insular systems might be biased due to underlying positive correlation between island area and heterogeneity. To test this hypothesis we surveyed the literature for studies providing estimates of both

environmental heterogeneity and area in insular systems, and compared the effects of environmental heterogeneity on species richness before and after accounting for the effect of area in the analysis. In contrast to most previous analyses of this sort, we included the quadratic term of environmental heterogeneity in the models. We predicted that accounting for the effect of area would significantly reduce the frequency of positive relationships and increase the frequency of non-significant and unimodal relationships in the data.

C.2 Data source To be compatible with the analyses performed in this study we searched for studies in which environmental heterogeneity was expressed by elevation range. We found 40 such studies (20,22,52-89). These studies included 66 species-island datasets, but some datasets had partial overlap in their species composition (e.g., endemic species vs. all species, woody species vs. all plant species, native vs. all species, etc.). We therefore used only a single dataset for each group of species in each study (the one that was not nested within any other dataset) for our analysis. This selection resulted in 54 different datasets.

C.3 Statistical analysis In order to test our hypothesis we first identified all datasets for which linear regression analysis revealed statistically significant effect of elevation range on species richness (a total of 43 datasets, SI Appendix Table 7). Each of these datasets was then analyzed in three steps. First, elevation range was regressed against area and the residuals of the model were determined. Second, log species richness was regressed against log area and the predicted values of log species richness were back transformed into units of species richness in order to determine the residuals of species richness over area. In the third step,

the residuals of species richness over area were regressed against the residuals of elevation range over area using both linear and quadratic regression models (note that this analysis is equivalent to the calculation of partial regression coefficients in multiple regression analysis and was necessary to account for logarithmic species-area relationships). A comparison of the results obtained for the quadratic vs. linear models using the sample size-corrected Akaike Information Criterion (AICC) enabled us to determine which of the two models better fits each dataset. The results (SI Appendix Table 7, Fig. 4A) showed that only 6 out of the 43 datasets remained significantly positive after the effect of area was adjusted for. The remaining datasets become significantly unimodal (14 patterns), non-significant (20 patterns), negative (1 pattern) or U-shaped (2 patterns). As a more formal test for the prevalence of unimodal relationships in the data we calculated the average effect size of all area-corrected quadratic models based on their Fvalues. We first transformed the F-statistic into rF, an estimate of effect size, using the formula rF F ( F df ) . The sign of the rF-values was determined by the sign of the quadratic regression coefficients. Values of rF were then transformed to Fisher's z-values for determining the mean and 95% confidence intervals of the effect size (90). The results (mean: -0.1785; confidence intervals: -0.1129 to -0.2441) indicated that the average effect size was significantly smaller than zero, supporting the hypothesis that the general relationship between species richness and elevation range was unimodal. Eighteen of the datasets we found also provided data on habitat diversity. We therefore analyzed these datasets using the same analysis performed for elevation range but with habitat diversity as a measure of environmental heterogeneity. The results were very similar to those obtained for elevation range (Fig. 4B). These overall findings support our hypothesis that previously reported evidence for positive effects of environmental heterogeneity on species richness in insular systems might be biased.

C.4 Criteria used to determine the functional relationship between species richness and the area-corrected elevation range A dataset was classified as unimodal if it satisfied the following requirements: (1) the quadratic model was statistically significant (P < 0.05), (2) the quadratic term was negative, (3) the inflection point was within the range of the data, and (4) the quadratic model showed a smaller value of the sample size-corrected Akaike Information Criterion (AICC) than the corresponding linear model. If requirements 1,3 and 4 were satisfied and the quadratic term was positive, the dataset was classified as U-shaped. Datasets were classified as positive (negative) if: (1) the dataset was not classified as unimodal, (2) the linear model was statistically significant (P < 0.05), (3) the linear term in the linear model was positive (negative) or: (1) the quadratic model was statistically significant (P < 0.05), (2) the quadratic term was negative, (3) the inflection point was to the right (left) of the range of the data or: (1) the quadratic model was statistically significant (P < 0.05), (2) the quadratic term was positive, (3) the inflection point was to the left (right) of the range of the data If the dataset was not classified as either unimodal, U-shaped, positive or negative, it was classified as non-significant.

D. References 1.

Rosenzweig ML (1995) Species diversity in space and time (University of Cambridge press, New York).

2.

Hanski I (1999) Metapopulation ecology (Oxford University Press).

3.

MacArthur RH & Wilson EO (1967) The theory of island biogeography (Princeton University Press, Princeton, NJ).

4.

Bell G (2000) The distribution of abundance in neutral communities. Am. Nat. 155(5):606-617.

5.

S. P. Hubbell, The unified neutral theory of biodiversity and biogeography. Monographs in Population Biology (Princeton University Press, Princeton, NJ, 2001).

6.

Hutchinson GE (1957) Concluding remarks, Cold Spring Harbor Symposium. Quantitative Biology 22:415-427.

7.

Chase JM & Leibold MA (2003) Ecological niches: linking classical and contemporary approaches (University of Chicago Press, Chicago, IL).

8.

Chase JM (2005) Towards a really unified theory for metacommunities. Funct. Ecol. 19(1):182-186.

9.

Haegeman B & Loreau M (2011) A mathematical synthesis of niche and neutral theories in community ecology. J. Theor. Biol. 269(1):150-165.

10.

Diamond JM (1976) Island biogeography and conservation - strategy and limitations. Science 193(4257).

11.

Simberloff D (1988) The contribution of population and community biology to conservation science. Annu. Rev. Ecol. Syst. 19:473-511.

12.

Benton TG, Vickery JA, & Wilson JD (2003) Farmland biodiversity: is habitat heterogeneity the key? Trends Ecol. Evol. 18(4):182-188.

13.

Flather CH, Hayward GD, Beissinger SR, & Stephens PA (2011) Minimum viable populations: is there a 'magic number' for conservation practitioners? Trends Ecol. Evol. 26(6):307-316.

14.

Kadmon R & Allouche O (2007) Integrating the effects of area, isolation, and habitat heterogeneity on species diversity: A unification of island biogeography and niche theory. Am. Nat. 170(3):443-454.

15.

Allouche O & Kadmon R (2009) A general framework for neutral models of community dynamics. Ecol. Lett. 12(12):1287-1297.

16.

McKinney ML (1997) Extinction vulnerability and selectivity: Combining ecological and paleontological views. Annu. Rev. Ecol. Syst. 28:495-516.

17.

Muntaner J, X. Ferrer, and A. Martı´nez-Vilalta (1984) Atles dels ocells nidificants de Catalunya i Andorra [Breeding Birds Atlas of Catalonia and Andorra] (Ketres, Barcelona, Spain).

18.

Estrada J, Pedrocchi, V., Brotons, L. & Herrando, S. (2004) Atles dels ocells nidificants de Catalunya 1999–2002 [Breeding Birds Atlas of Catalonia] (Institut Català d’Ornitologia (ICO)/Lynx Editions, Barcelona, Spain).

19.

Gravel D, Canham CD, Beaudet M, & Messier C (2006) Reconciling niche and neutrality: the continuum hypothesis. Ecol. Lett. 9(4):399-409.

20.

Ricklefs RE & Lovette IJ (1999) The roles of island area per se and habitat diversity in the species-area relationships of four Lesser Antillean faunal groups. J. Anim. Ecol. 68(6):1142-1160.

21.

Moreno Saiz JC & Lobo JM (2008) Iberian-Balearic fern regions and their explanatory variables. Plant Ecol. 198(2):149-167.

22.

Veech JA & Crist TO (2007) Habitat and climate heterogeneity maintain betadiversity of birds among landscapes within ecoregions. Glob. Ecol. Biogeogr. 16(5):650-656.

23.

Tamme R, Hiiesalu I, Laanisto L, Szava-Kovats R, & Partel M (2010) Environmental heterogeneity, species diversity and co-existence at different spatial scales. J. Veg. Sci. 21(4):796-801.

24.

Kerr JT & Packer L (1997) Habitat heterogeneity as a determinant of mammal species richness in high-energy regions. Nature 385(6613):252-254.

25.

Costa GC, Nogueira C, Machado RB, & Colli GR (2007) Squamate richness in the Brazilian Cerrado and its environmental-climatic associations. Divers. Distrib. 13(6):714-724.

26.

Huang Y, et al. (2011) Lizard species richness patterns in China and its environmental associations. Biodivers. Conserv. 20(7):1399-1414.

27.

Rahbek C & Graves GR (2001) Multiscale assessment of patterns of avian species richness. Proc. Natl. Acad. Sci. U. S. A. 98(8):4534-4539.

28.

Jimenez I, Distler T, & Jorgensen PM (2009) Estimated plant richness pattern across northwest South America provides similar support for the species-energy and spatial heterogeneity hypotheses. Ecography 32(3):433-448.

29.

Montigny MK & MacLean DA (2005) Using heterogeneity and representation of ecosite criteria to select forest reserves in an intensively managed industrial forest. Biol. Conserv. 125(2):237-248.

30.

Kati V, et al. (2010) Towards the use of ecological heterogeneity to design reserve networks: a case study from Dadia National Park, Greece. Biodivers. Conserv. 19(6):1585-1597.

31.

Fuhlendorf SD, et al. (2006) Should heterogeneity be the basis for conservation? Grassland bird response to fire and grazing. Ecol. Appl. 16(5):1706-1716.

32.

Tilman D (2004) Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community assembly. Proc. Natl. Acad. Sci. U. S. A. 101(30):10854-10861.

33.

Ruokolainen L, Ranta E, Kaitala V, & Fowler MS (2009) When can we distinguish between neutral and non-neutral processes in community dynamics under ecological drift? Ecol. Lett. 12(9):909-919.

34.

Wasser, D. E. & Sherman, P. W. Avian longevities and their interpretation under evolutionary theories of senescence. Journal of Zoology 280, 103-155 (2010).

35.

Carnicer, J., Brotons, L., Sol, D. & Jordano, P. Community-based processes behind species richness gradients: contrasting abundance-extinction dynamics and sampling effects in areas of low and high productivity. Glob. Ecol. Biogeogr. 16, 709-719 (2007).

36.

Moreno-Rueda, G. & Pizarro, M. Relative influence of habitat heterogeneity, climate, human disturbance, and spatial structure on vertebrate species richness in Spain. Ecol. Res. 24, 335-344 (2009).

37.

Sutherland, G. D., Harestad, A. S., Price, K. & Lertzman, K. P. Scaling of Natal dispersal distance in terrestrial birds and mammals. Conservation Ecology (online) 4, 16 (2000).

38.

Kunin, W. E. Biodiversity at the edge: A test, of the importance of spatial "mass effects" in the Rothamsted Park Grass experiments. Proc. Natl. Acad. Sci. U. S. A. 95, 207-212 (1998).

39.

Ferrer, X., Carrascal, L. M., Gordo, O. & Pino, J. Bias in avian sampling effort due to human preferences: An analysis with Catalonian birds (1900-2002). Ardeola 53, 213-227 (2006).

40.

Legendre, P. & Legendre, L. Numerical Ecology. (Elsevier Science B.V., 1998).

41.

Dormann, C. F. et al. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. Ecography 30, 609-628 (2007).

42.

Sokal, R. R. Testing statistical significance of geographic-variation patterns. Syst. Zool. 28, 227-232 (1979).

43.

Rangel, T. F., Diniz-Filho, J. A. F. & Bini, L. M. SAM: a comprehensive application for Spatial Analysis in Macroecology. Ecography 33, 46-50 (2010).

44.

Burnham, K. P. & D.R.Anderson. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. (Springer-Verlag, 2002).

45.

McCain, C. M. Global analysis of bird elevational diversity. Glob. Ecol. Biogeogr. 18, 346-360 (2009).

46.

Pulliam, H. R. On the relationship between niche and distribution. Ecol. Lett. 3, 349-361 (2000).

47.

Whittaker, R. H. Vegetation of the Great Smoky Mountains. Ecol. Monogr. 26, 180 (1956).

48.

Allouche, O. & Kadmon, R. A general framework for neutral models of community dynamics. Ecol. Lett. 12, 1287-1297 (2009).

49.

Levins, R. Some demographic and genetic consequences of environmental heterogeneity for biological control. Bull. Entomol. Soc. Am. 15, 227-240 (1969).

50.

Durrett, R. & Levin, S. The importance of being discrete (and spatial). Theor. Popul. Biol. 46, 363-394 (1994).

51.

Berec, L. Techniques of spatially explicit individual-based models: construction, simulation, and mean-field analysis. Ecol. Model. 150, 55-81 (2002).

52.

Abbott, I. Numbers of plant, insect and land bird species on nineteen remote islands in southern-hemisphere. Biol. J. Linn. Soc. 6, 143-152 (1974).

53.

Ackerman, J. D., Trejo-Torres, J. C. & Crespo-Chuy, Y. Orchids of the West Indies: predictability of diversity and endemism. J. Biogeogr. 34, 779-786 (2007).

54.

Cardoso, P., Arnedo, M. A., Triantis, K. A. & Borges, P. A. V. Drivers of diversity in Macaronesian spiders and the role of species extinctions. J. Biogeogr. 37, 10341046 (2010).

55.

Case, T. J. Species numbers, density compensation, and colonizing ability of lizards on islands in gulf of california. Ecology 56, 3-18 (1975).

56.

Case, T. J. & Bolger, D. T. The role of introduced species in shaping the distribution and abundance of island reptiles. Evol. Ecol. 5, 272-290 (1991).

57.

Chen, Y.-H. Combining the Species-Area-Habitat Relationship and Environmental Cluster Analysis to Set Conservation Priorities: a Study in the Zhoushan Archipelago, China. Conserv. Biol. 23, 537-545 (2009).

58.

Duarte, M. C., Rego, F., Romeiras, M. M. & Moreira, I. Plant species richness in the Cape Verde islands - eco-geographical determinants. Biodivers. Conserv. 17, 453-466 (2008).

59.

Dueser, R. D. & Brown, W. C. Ecological correlates of insular rodent diversity. Ecology 61, 50-56 (1980).

60.

Fattorini, S. Biogeography of tenebrionid beetles (Coleoptera: Tenebrionidae) in the circum-Sicilian islands (Italy, Sicily): Multiple biogeographical patterns require multiple explanations. Eur. J. Entomol. 108, 659-672 (2002).

61.

Hamilton, T. H., Rubinoff, I., Barth, R. H., Jr. & Bush, G. L. Species Abundance: Natural Regulation of Insular Variation. Science 142, 1575-1577 (1963).

62.

Hannus, J.-J. & von Numers, M. Vascular plant species richness in relation to habitat diversity and island area in the Finnish Archipelago. J. Biogeogr. 35, 10771086 (2008).

63.

Harris, M. P. Galapagos avifauna. Condor 75, 265-278 (1973).

64.

Johnson, M. P., Simberloff D. S. Environmental Determinants of Island Species Numbers in the British Isles. J. Biogeogr. 1, 149-154 (1974).

65.

Johnson, M. P., Mason, L. G. & Raven, P. H. Ecological parameters and plant species diversity. Am. Nat. 102, 297-& (1968).

66.

Johnson, M. P. & Raven, P. H. Species number and endemism - galapagos archipelago revisited. Science 179, 893-895 (1973).

67.

Kratter, A. W. Montane avian biogeography in southern california and bajacalifornia. J. Biogeogr. 19, 269-283 (1992).

68.

Losos, J. B. Island biogeography of day geckos (phelsuma) in the indian-ocean. Oecologia 68, 338-343 (1986).

69.

McMaster, R. T. Factors influencing vascular plant diversity on 22 islands off the coast of eastern North America. J. Biogeogr. 32, 475-492 (2005).

70.

Millien-Parra, V. & Jaeger, J. J. Island biogeography of the Japanese terrestrial mammal assemblages: an example of a relict fauna. J. Biogeogr. 26, 959-972 (1999).

71.

Nakamura, K., Suwa, R., Denda, T. & Yokota, M. Geohistorical and current environmental influences on floristic differentiation in the Ryukyu Archipelago, Japan. J. Biogeogr. 36, 919-928 (2009).

72.

Newmark, W. D. Species area relationship and its determinants for mammals in western north-american national-parks. Biol. J. Linn. Soc. 28, 83-98 (1986).

73.

Panitsa, M., Trigas, P., Iatrou, G. & Sfenthourakis, S. Factors affecting plant species richness and endemism on land-bridge islands - An example from the East Aegean archipelago. Acta Oecologica-International Journal of Ecology 36, 431437 (2010).

74.

Panitsa, M., Tzanoudakis, D., Triantis, K. A. & Sfenthourakis, S. Patterns of species richness on very small islands: the plants of the Aegean archipelago. J. Biogeogr. 33, 1223-1234 (2006).

75.

Perry, G., Rodda, G. H., Fritts, T. H. & Sharp, T. R. The lizard fauna of Guam's fringing islets: island biogeography, phylogenetic history, and conservation implications. Glob. Ecol. Biogeogr. 7, 353-365 (1998).

76.

Power, D. M. Numbers of bird species on california islands. Evolution 26, 451-& (1972).

77.

Power, D. M. Avifauna richness on california channel islands. Condor 78, 394-398 (1976).

78.

Price, J. P. Floristic biogeography of the Hawaiian Islands: influences of area, environment and paleogeography. J. Biogeogr. 31, 487-500 (2004).

79.

Reed, T. The number of breeding land-bird species on british islands. J. Anim. Ecol. 50, 613-624 (1981).

80.

Rogers, G. & Overton, J. Regional patterns of plant species richness in southern New Zealand. N. Z. J. Bot. 38, 609-627 (2000).

81.

Sfenthourakis, S. The species-area relationship of terrestrial isopods (Isopoda; Oniscidea) from the Aegean archipelago (Greece): A comparative study. Global Ecology and Biogeography Letters 5, 149-157 (1996).

82.

Sfenthourakis, S. & Triantis, K. A. Habitat diversity, ecological requirements of species and the Small Island Effect. Divers. Distrib. 15, 131-140 (2009).

83.

Vanderwerff, H. Species number, area and habitat diversity in the galapagosislands. Vegetatio 54, 167-175 (1983).

84.

Vuilleumier, F. Insular biogeography in continental regions .1. northern andes of south-america. Am. Nat. 104, 373-& (1970).

85.

Weissman, D. B. & Rentz, D. C. Zoogeography of the Grasshoppers and Their Relatives (Orthoptera) on the California Channel Islands. J. Biogeogr. 3, 105-114 (1976).

86.

Welter-Schultes, F. W. & Williams, M. R. History, island area and habitat availability determine land snail species richness of Aegean islands. J. Biogeogr. 26, 239-249 (1999).

87.

White, P. S., Miller, R. I. & Ramseur, G. S. The species-area relationship of the southern Appalachian high peaks - vascular plant richness and rare plantdistributions. Castanea 49, 47-61 (1984).

88.

Wilcox, B. A., Murphy, D. D., Ehrlich, P. R. & Austin, G. T. Insular biogeography of the montane butterfly faunas in the great-basin - comparison with birds and mammals. Oecologia 69, 188-194 (1986).

89.

Yeakley, J. A. & Weishampel, J. F. Multiple source pools and dispersal barriers for Galapagos plant species distribution. Ecology 81, 893-898 (2000).

90.

Borenstein, M., Hedges, L. V., Higgins, J. P. T. & Rothstein, H. R. Introduction to Meta-Analysis. (John Wiley & Sons, 2009).

SI APPENDIX TABLES SI Appendix Table 1. Results of regression models testing the effects of elevation range (H), the squared value of elevation range (H2), area (A), geographical coordinates (X, Y, X2, Y2, XY), mean absolute elevation (E), mean annual temperature (T), mean annual rainfall (R), and human population size (P) on mean species abundance of breeding birds in Catalonia. Analyses were performed using four different models: OLS incorporating the simultaneous effects of the linear term and squared term of elevation range (OLS – basic), OLS incorporating all independent variables with backward variable selection (OLS – extended), Conditional Autoregressive models (CAR) incorporating those variables that were found statistically significant in the OLS models, and corresponding Simultaneous Autoregressive models (SAR). For each model we report the t-values of each variable, the respective significance level (bold values indicate significant effects), the R2 value of the model, the AICC value, and the difference between AICC value of the relevant model and the 'best' model (the one with the lowest value).

H H2 A X Y X2 Y2 XY E R T

OLS (Basic) OLS (Extended) 0.69 -----4.13***

CAR -----

SAR -----

-4.30*** 4.70*** -----4.65*** -7.83*** ----7.66***

-4.57*** 3.53*** -----1.15 3.41*** ----3.27***

-3.91*** 2.07 ----0.29 -1.11 ----0.95

-9.83*** ---------

-10.30*** ---------

-10.89*** ---------

P -4.05*** -4.15*** -4.24*** ----------------------------------------------------------------------------------------------R2 0.32 0.59 0.54 0.63 AICC -1668.6 -1848 -1802.7 -1880.6 ΔAICC 212 32.6 77.9 0 ----------------------------------------------------------------------------------------------* P < 0.05, ** P < 0.01, *** P < 0.001

SI Appendix Table 2. Results of regression models testing the effects of elevation range (H), the squared value of elevation range (H2), area (A), geographical coordinates (X, Y, X2, Y2, XY), mean absolute elevation (E), mean annual temperature (T), mean annual rainfall (R), and human population size (P) on extinction rates of breeding birds in Catalonia. Analyses were performed using four different models: OLS incorporating the simultaneous effects of the linear term and squared term of elevation range (OLS – basic), OLS incorporating all independent variables with backward variable selection (OLS – extended), Conditional Autoregressive models incorporating those variables that were found statistically significant in the OLS models (CAR), and corresponding Simultaneous Autoregressive models (SAR). For each model we report the t-values of each variable, the respective significance level (bold values indicate significant effects), the R2 value of the model, the AICC value, and the difference between AICC values of the relevant model and the 'best' model (the one with the lowest value).

H H2 A

OLS (Basic) OLS (Extended) 3.08** 4.63*** -1.54 -----4.09***

CAR 4.64*** -----2.12*

SAR 3.54*** -----1.36

X

-----

-----

-----

Y X2 Y2 XY E R

7.38*** ---------------------

3.10** ---------------------

3.81*** ---------------------

T 4.56*** 2.47* 2.21* P -------------------------------------------------------------------------------------------------------------------R2 0.32 0.23 0.31 0.40 AICC 2587.4 2528.9 2493.4 2442.0 ΔAICC 145.4 86.9 51.4 0 ----------------------------------------------------------------------------------------------* P < 0.05, ** P < 0.01, *** P < 0.001

SI Appendix Table 3. Results of regression models testing the effects of elevation range (H), the squared value of elevation range (H2), area (A), geographical coordinates (X, Y, X2, Y2, XY), mean absolute elevation (E), mean annual temperature (T), mean annual rainfall (R), and human population size (P) on species richness of breeding birds in Catalonia. Analyses were performed using four different models: OLS incorporating the simultaneous effects of the linear term and squared term of elevation range (OLS – basic), OLS incorporating all independent variables with backward variable selection (OLS – extended), Conditional Autoregressive models incorporating those variables that were found statistically significant in the OLS models (CAR), and corresponding Simultaneous Autoregressive models (SAR). For each model we report the t-values of each variable, the respective significance level (bold values indicate significant effects), the R2 value of the model, the AICC value, and the difference between AICC values of the relevant model and the 'best' model (the one with the lowest value).

H H2 A X

OLS (Basic) OLS (Extended) 6.27*** 5.89*** -5.07*** -3.63*** 7.59*** -----

Y X2 Y2 XY

5.70*** ---------3.18**

CAR 4.86*** -2.63** 6.77*** -----

SAR 5.08*** -3.13*** 6.86*** -----

3.06** ---------2.09**

0.53 ---------1.43

E -4.53*** -3.24*** -3.26*** R ------------T ------------P -------------------------------------------------------------------------------------------------------------------R2

0.13

0.31

0.35

0.36

AICC 2985.2 2904.6 2884.3 2882.6 ΔAICC 102.6 22.0 1.7 0 ----------------------------------------------------------------------------------------------* P < 0.05, ** P < 0.01, *** P < 0.001

SI Appendix Table 4. Results of regression models testing the effects of elevation range (H) and its squared value (H2) on mean species abundance, extinction rates, and species richness, based on 189 grid cells for which sampling effort was similar in the two survey periods (1983 and 2002). Values reported are t-values of each variable and their significance levels (bold values indicate significant effects).

H H2

Mean abundance 0.616 -3.430***

Extinction rates 2.647** -1.708

Species richness 5.452*** -4.735***

----------------------------------------------------------------------------------------------* P < 0.05, ** P < 0.01, *** P < 0.001

SI Appendix Table 5. Results of GLM testing the effects of elevation range (H) and its squared value (H2) on species richness in the first survey (1983) with (left) and without (right) including the effect of sampling effort (SE, five categories based on the number of sampling hours) in the model. Note that sampling effort had a positive effect on species richness but did not affect the qualitative results obtained for the effect of elevation range. See SI Appendix Fig. 2 for a graphical presentation of the relationship between species richness and elevation range.

H H2 SE < 40h 40h < SE < 80h 80h < SE < 120h 120h < SE < 160h SE > 160h

Without correction Coefficient t 43.8 7.21*** -16.8 -5.31***

With correction Coefficient t 29.2 5.66*** -11.66 -4.35*** -36.6 -4.12*** -34.7 -3.95**** -19.7 -2.22* -11.8 -1.32 (Set to zero)

SI Appendix Table 6. Results of GLM testing the effects of elevation range (H) and its squared value (H2) on species richness in the second survey with (left) and without (right) including the effect of sampling effort (SE, five categories based on the number of sampling hours) in the model. Note that sampling effort had a positive effect on species richness but did not affect the qualitative results obtained for the effect of elevation range. See SI Appendix Fig. 2 for a graphical presentation of the relationship between species richness and elevation range.

H H2 SE < 40h 40h < SE < 80h 80h < SE < 120h 120h < SE < 160h SE > 160h

Without correction Coefficient t

With correction Coefficient t

23.8 -11.7

14.2 -6.4 -31.1 -26.3 -16.5 -10.0 (Set to zero)

3.88*** -3.65***

2.70** -2.34* -8.50*** -7.41*** -4.75*** -2.73**

SI Appendix Table 7. Results of linear and quadratic regression analyses for the area-corrected effect of elevation range on species richness in 43 datasets of insular systems. In both models the response variable is the residuals of species richness over area and the independent variable is the residuals of elevation range over area (including its squared term in the quadratic model). The model that better fits the data was determined as the one with the smaller value of the sample size-corrected Akaike Information Criterion (AICC) provided that the pattern was statistically significant (see SI Methods section C.4 for full details).

Linear model

Quadratic model

Area Group Birds

R2

[hectare]

AICC

R2

P-value

AICC

P-value Pattern

Source

9,000-151,000

0.125

-57.3

0.0758

0.323

-60.8

0.0172

Unimodal

Ricklefs & Lovette 1999

9,000-151,000

0.229

-59.7

0.0221

0.322

-60.8

0.0175

Unimodal

Ricklefs & Lovette 1999

Butterflies

9,000-151,000

0.447

-66.4

0.0038

0.428

-64.7

0.0139

Positive

Ricklefs & Lovette 1999

Birds

0.8-10,692

0.0758

-318

0.0105

0.141

-322

0.0018

Unimodal

Reed 1981

Plants

3,500-99,100

0.37

-26.3

0.0367

0.76

-34.1

0.0028

Unimodal

Duarte et al. 2008

Spiders

300-205,800

0.172

-62.5

0.0394

0.184

-61.5

0.069

Positive

Cardoso et al. 2010

Birds

215-46,870

-0.0769

-38.3

0.995

-0.166

-35.6

0.998

NS

Chen 2009

0.05-68.3

0.0147

-362

0.141

0.0044

-360

0.313

NS

Hannus & von Numers 2008

Reptiles and amphibians

Plants, Brunskar archipelago

Plants, Getskar

0.004-14.314

0.0248

-341

0.0897

0.0599

-343

0.0368

Unimodal

Hannus & von Numers 2008

Mammals

1.22-30.51

0.0329

-51.5

0.227

0.0070

-49.6

0.371

NS

Wilcox et al. 1986

Butterflies

1.22-30.51

-0.0351

-50.2

0.524

-0.066

-48.3

0.632

NS

Wilcox et al. 1986

Birds

1.22-30.51

0.525

-41.8

0.0030

0.478

-39

0.0156

Positive

Wilcox et al. 1986

Birds

2,500-348,700

-0.075

-38.3

0.88

-0.119

-36.2

0.779

NS

Vuilleumier 1970

Plants

0.05-50

0.115

-393

0.0008

0.19

-399

5.78E-05

Unimodal

Panitsa et al. 2006

0.0942

-77.5

0.079

0.109

-76.6

0.115

NS

Newmark 1986

archipelago

Mammals

8,3002,073,600

Plants

189-274,063

0.234

-239

9.31E-05

0.239

-239

0.0003

Positive

Rogers and Overton 2000

Terrestrial isopods

0.23-826,120

0.12

-415

0.0005

0.168

-419

0.0001

U-shaped

Sfenthourakis & Triantis 2009

0.225

-132

0.0024

0.204

-130

0.0098

Positive

Veech & Crist 2007

0.556

-218

1.12E-09

0.605

-222

5.11E-10

U-shaped

Ackerman et al. 2007

Birds

Orchids

960,00075,290,000 2,10010,500,700

Lizards

60-143,000

-0.0047

-75

0.355

0.177

-78.6

0.0497

Unimodal

Case 1975

Plants

60-143,000

0.109

-77.9

0.064

0.702

-103

1.17E-06

Unimodal

Case 1975

Rodents

29-2,197

0.337

-22.1

0.0593

0.23

-18.7

0.193

NS

Dueser & Brown 1980

Land birds

259-582,488

0.139

-41.7

0.0944

0.287

-43

0.052

NS

Harris 1973

Plants

259-582,488

-0.0769

-38.3

0.995

0.216

-41.6

0.0918

NS

Harris 1973

Plants

0.5-229,850

-0.0236

-155

0.819

-0.0385

-153

0.788

NS

Johnson & Simberloff 1974

Day geckos

0.4-552,500

0.0781

-189

0.0305

0.192

-194

0.0031

Unimodal

Losos 1986

Plants

52-34,706

0.298

-48.8

0.0169

0.281

-47

0.0462

Negative

Power 1972

Birds

2,600-24,900

0.214

-17.1

0.139

0.209

-14.8

0.24

NS

Power 1976

Plants

2,600-24,900

0.335

-18.5

0.0775

0.212

-14.8

0.237

NS

Power 1976

Plants

52-582,488

-0.0621

-46

0.804

0.299

-51.6

0.0325

Unimodal

Hamilton 1963

Orthoptera

2.6-249

-0.145

-14.1

0.747

-0.368

-10.4

0.943

NS

Weissman & Rentz 1976

Plants

1-466,930

-0.0461

-70

0.863

0.34

-79.3

0.0060

Unimodal

Johnson & Raven 1973

Plants

3-26,668

-0.0498

-65.8

0.953

-0.0893

-63.7

0.871

NS

McMaster 2005

Plants

1-1,636

0.157

-62.2

0.0473

0.144

-60.6

0.103

Positive

Panitsa et al. 2010

Snails

0.58-826,000

-0.0149

-269

0.806

0.0222

-271

0.186

NS

Welter-Schultes & Williams 1999

Plants

1-466,930

-0.0072

-96.4

0.379

0.333

-107

0.0020

Unimodal

Yeakley & Weishampel 2000

Tenebrionid beetles

380-826,000

-0.0312

-109

0.806

-0.0537

-107

0.812

NS

Fattorini 2002

Plants

5-34,706

0.268

-40.1

0.0337

0.2050

-37.4

0.1130

NS

Johnson et al. 1968

Mammals

400-23,051,000

-0.0484

-8.81

0.43

-0.295

-4.26

0.685

NS

Millien-Parra & Jaeger 1999

Plants

2-1,043,300

0.176

-54.4

0.0468

0.272

-55.2

0.0361

Unimodal

Price 2004

Terrestrial Isopods

3-47,620

-0.0239

-160

0.885

-0.0212

-159

0.574

NS

Sfenthourakis 1996

Plants

23-467,000

-0.026

-50.4

0.462

0.345

-57.1

0.0163

Unimodal

Vanderwerff 1983

Plants

3-4,900

0.0251

-22

0.299

-0.0015

-19.8

0.417

NS

White 1984

Extinction rates

Mean abundance

Species richness

SI APPENDIX FIGURES

Environmental heterogeneity SI Appendix Figure 1. Effect of elevation range on mean species abundance, extinction rates, and species richness, in 189 grid cells that had similar sampling efforts in the two survey periods (1983 and 2002). All patterns show the same trends as in the analysis of the full dataset (Fig. 2) and in spite of a 50% reduction in sample size are highly significant (SI Appendix Table 4)

Species richness

• 1983 • 2002

Environmental heterogeneity

SI Appendix Figure 2. Effect of elevation range on species richness in the first (1983) and second (2002) surveys based on the complete dataset (372 grid cells). For both surveys the relationship was significantly unimodal both before and after correction for spatial variation in sampling effort (SI Appendix Tables 5-6).

SI Appendix Figure 3. Effect of environmental heterogeneity on extinction probability of species with different niche widths. Data refer to breeding birds in Catalonia with elevation range as a measure of heterogeneity. Each color represents a different quartile of niche width. Increasing niche width decreases the slope of the regression line. P-values indicate significance levels of the corresponding linear regression models.

Environmental range (regional)

Environmental range (local)

SI Appendix Figure 4. Schematic presentation of species responses to variation in environmental variation in our model. The environment is reduced to a single variable E and each species i is characterized by a Gaussian response with optimum µi and standard deviation ('niche width') σi indicating its establishment probability as a function of E. The range of E available for the local community (ER = Emax-Emin) is embedded within the regional range ERR = ERmax-ERmin.

Species richness Species abundance Extinction rates

R = 30 R = 45 R = 90

Environmental heterogeneity (ER)

SI Appendix Figure 5. Effect of reproduction rates on the relationship between environmental heterogeneity and mean species abundance, extinction rates, and species richness, based on the stochastic model of community dynamics. Each dot represents the mean of 10 independent runs of 100,001 years. Parameters: N = 500, A = 1,000, ERR = 65, σ = 2, M = 1, I = 0.2.

Species richness Species abundance Extinction rates

I = 0.05 I = 0.20 I = 0.50

Environmental heterogeneity (ER)

SI Appendix Figure 6. Effect of immigration rates on the relationship between environmental heterogeneity and mean species abundance, extinction rates, and species richness, based on the stochastic model of community dynamics. Each dot represents the mean of 10 independent runs of 100,001 years. Parameters: N = 500, A = 1,000, ERR = 65, σ = 2, M = 1, R = 45.

Log Species richness

A A = 100 A = 1000 A =50000

Log Species richness

B TE = 1 TE = 5 TE = 20

Environmental heterogeneity (ER) SI Appendix Figure 7. Effect of area (A = number of sites) and temporal extent of data collection (TE = number of years over which the data are recorded) on the relationship between environmental heterogeneity and species richness, based on the stochastic model of community dynamics. Note that extending the spatial and temporal scale of data collection may turn a unimodal pattern into a positive one. Each dot represents the mean of 10 independent runs of 100,001 years. Parameters: N = 500, A = 1,000, R = 45, M = 1, I = 0.2, σ = 2, ERR = 65.

Niche width = 4

I = 0.01 I = 0.05 I = 0.10

Niche width = 12

Niche width = 30

Environmental heterogeneity (ER)

Environmental heterogeneity (ER)

Log Species richness

Log Species richness

Niche width = 1

SI Appendix Figure 8. Effect of niche width (σ) and immigration (I) on the relationship between environmental heterogeneity and species richness based on the stochastic model of community dynamics. Note that increasing niche width may turn a predominantly negative pattern into a positive one with unimodal patterns occurring at intermediate niche widths. Each dot represents the mean of 50 independent runs of 100,001 years. Parameters: N = 1,500, A = 5,000, R = 100, M = 1, ERR = 150.

Species richness Species abundance Extinction rates Environmental heterogeneity (ER)

SI Appendix Figure 9. Simulated responses of ecological communities to environmental heterogeneity after different time intervals. Samples taken after 20,000 years (as in the manuscript, allowing system to reach steady state, blue) and 50 years (orange). Each dot represents the average of 10 independent simulation runs for the 20,000 years period, and 2500 independent runs for the 50 years period. Parameter values are identical to those used in Fig. 3 of the manuscript: A (area) = 1000, R (reproduction) = 45, M (mortality) = 1, I (immigration) = 0.2, and σ (niche width) = 2.

The Area-Heterogeneity Trade-off and the Diversity of Ecological Communities

Omri Allouchea*, Michael Kalyuzhnya, Gregorio Moreno-Ruedab, Manuel Pizarrob, Ronen Kadmona

a

Department of Ecology, Evolution & Behavior, Institute of Life Sciences, The Hebrew University of Jerusalem, 91904, Jerusalem, Israel.

b

Departamento de Zoología, Facultad de Ciencias, Universidad de Granada, 18071, Granada, Spain.

* To whom correspondence should be addressed. E-mail: [email protected]

SUPPORTING INFORMATION APPENDIX METHODS A. Empirical analyses A.1

General considerations

A.2

Study system

A.3

Quantification of the response variables

A.4

Data analysis A.4.1 Community responses to elevation range A.4.2 Effects of niche width

B. Theoretical analyses B.1

Objectives

B.2

Description of the model

B.3

Simulation experiments

C. Meta-analysis C.1

Objectives

C.2

Data source

C.3

Statistical analysis

C.4

Criteria used for functional relationship determination

D. Additional References

A. Empirical analysis A.1 General considerations A comprehensive test of the area-heterogeneity trade-off hypothesis requires data on species richness, extinction rates, mean species abundance, and an appropriate measure of environmental heterogeneity. While data on species richness, species abundance, and environmental heterogeneity can be obtained for a single point in time, determining extinction rates requires data on species composition in a fixed set of sites for (at least) two distinct survey periods. The interval between the two surveys should be larger than the typical life span of the relevant species to allow detectable levels of extinction, but not too large to minimize opportunities for recolonization following extinction events. The spatial extent of the observation unit is also important – it should be small enough to capture the footprints of local stochastic extinctions, but large enough to minimize the impact of immigration and emigration on the results and to obtain sufficient spatial variation in environmental conditions. Choosing an appropriate measure of environmental heterogeneity is also not trivial. Ideally, environmental heterogeneity should be defined with respect to an environmental factor that is known to affect the distribution of the target species in the relevant system. Taking into account that the strength of the areaheterogeneity trade-off depends on niche width (Fig. 1), it should also be possible to quantify niche width of the target species with respect to the selected factor, and there should be evidence that the study species differ in both the position and the width of their niches along the selected environmental gradient.

A.2 Study system Unfortunately, data satisfying all of the above requirements are extremely rare. In this study we used data on breeding bird distributions in Catalonia (NE Spain) as a basis for our analysis. This system was ideal for our purposes because (1) data were collected using the same set of observation units (UTM grid cells of 10x10 km) in two distinct surveys, the first from 1975-1982 with most data collected in 1980-82 (17) and the second from 19992002 (18); (2) the interval between the two surveys (two decades) was larger than the typical lifespan of most bird species (34); (3) the size of the observation unit (10x10km) was small enough to detect extinction events (35) but large enough to show a wide range of within-cell heterogeneity in environmental conditions; and (4) the data provided by the second survey included estimates of species abundance that were generated using logistic regression models based on a stratified random sampling of smaller units (1x1km) within the 10x10km cells. These estimates were available for species that were observed in at least 10 squares and for which cross validation revealed highly accurate predictions. The second atlas (18) also provided estimates of the effective number of hours invested in bird surveying within individual grid cells in each survey period (five categories: 160h), thus allowing to evaluate the potential effect of sampling bias on the results. Analyses of the data collected in the second survey has shown that elevation is a major determinant of breeding bird distributions in this area and that different species show considerable differences in both the position and the width of their niches along the elevation gradient (18). We therefore used elevation range as an effective measure of environmental heterogeneity in our analysis. A total of 372 grid cells for which appropriate data were available were selected for the analysis. These grid cells covered 99% of the terrestrial area of Catalonia.

Since many grid cells (179 out of 372) were border cells or coastal cells in which only part of the grid cell was surveyed, we calculated the terrestrial area of Catalonia in each grid cell by intersecting the UTM grid used to design the two field surveys with a vector map of Catalonia (downloaded from ArcGIS online world administrative divisions, http://www.esri.com/data/free-data/index.html). The 372 polygons were then intersected with a ~1 km resolution DTM of Catalonia (downloaded from http://dds.cr.usgs.gov/srtm/version1/Eurasia/) and the mean elevation of each polygon was determined. All GIS operations were performed using ESRI ArcGIS Desktop 10 (Environmental Systems Research Institute, Redlands, California, USA). Additional data on mean annual rainfall, mean annual temperature, and human population density of each grid cell were available from a previous analysis (36). Determining the appropriate scale for quantifying elevation range was more complicated due to differences in area among the peripheral polygons and the fact that the center of such polygons did not coincide with the center of the corresponding UTM cells. We therefore determined elevation range at three alternative scales using buffer zones with radii of 7km, 10km, and 14km around the center of each polygon. The three measures were highly correlated (Pearson correlation = 0.93-0.97) and had similar effects on species richness. We therefore used the smallest scale (7km) which is approximately one-half of the diagonal of a 10x10km cell, as a basis for quantifying an area-independent measure of elevation range. This radius ensured that all parts of all grid cells were included within the buffer. It also fits the finding of Sutherland et al. (37) that median distances of natal dispersal in birds range from 0.03 km to 10 km.

A.3 Quantification of the response variables Our predictions focused on three community-level characteristics: mean species abundance (Fig. 2A), extinction rates (Fig. 2B), and species richness (Fig. 2C). Species richness was determined for each grid cell as the average number of species that were present in that cell in the two survey periods. Estimates of species abundance were averaged for individual species in each UTM cell and were standardized prior to the analysis in order to control for differences in absolute abundance among species. Two standardization methods were applied, by maximum value and by the mean and standard deviation. In the first method, mean species abundance in a given grid cell was estimated as:

Pj

Pi , j

1 Nj

P i 1

max i

where Pj is the mean (standardized) abundance of all species present in grid cell j, Nj is the number of species for which estimates of abundance were available in grid cell j, Pi,j is the mean abundance of species i in grid cell j (i = 1,...,N), and Pmax i is the abundance of species i in the UTM cell with its highest abundance in Catalonia. In the second method mean species abundance in each cell was estimated as:

Pj

1 Nj

i

Pi , j i

i

where i is the mean abundance of species i in the study area and i is the corresponding standard deviation. We report the results obtained for the first (maximum-based) standardization because individual species abundances were not normally distributed and our analyses showed that both methods gave similar results. Yet, it should be noted that estimates of species abundance were based on predictive models of site occupancy (GLM regressions predicting probabilities of occurrence at a spatial resolution of 1km, 18) rather than direct observations. Thus, although occupancy is usually correlated with abundance

(38) and only species with accurate predictions were included in our analysis, these results should be interpreted with some cautious. Extinction rates were estimated for each grid cell by counting the number of species that were present in that cell in the first survey but not in the second survey (note that this definition includes local 'extinctions' caused by emigration to other cells). Although survey data like those analyzed in this study are known to suffer from imperfect detection, we believe that our estimates of extinction rates were not very sensitive to such bias because of the intensive sampling effort of the second survey (18). Spatial differences in sampling efforts might cause bias in breeding bird surveys (39), but this source of bias was minimized in the second survey by using a standardized sampling protocol, originally designed to minimize such bias. Moreover, a detailed analysis of the overall data collected from 1983-2002 showed that sampling effort in this period was positively and significantly correlated with elevation range (39). All expected consequences of such bias (higher species abundances, lower extinction rates and higher species richness at grid cells characterized by high elevation range) are exactly opposite to our predictions. We therefore believe that patterns of abundance and species richness supporting our predictions cannot be attributed to underlying bias in sampling effort. This assumption was confirmed by complementary analyses (see below).

A.4 Data analysis A.4.1 Community responses to elevation range The effect of elevation range on mean species abundance, extinction rates, and species richness was analyzed using four types of regression models. The first set of models tested the combined effects of elevation range (H) and its squared term (H2) on each of the three response variables using ordinary least squares (OLS) regressions. The squared value of

elevation range was included in the models to account for possible unimodal effects (though we do not expect any relationship to be explicit parabolic). These and all other OLS models were performed using SPSS for Windows, Ver. 15.0 (SPSS, Chicago, IL, U.S.A). The results of these analyses are presented in the first column of SI Appendix Tables 1-3. In the analysis of mean species abundance (SI Appendix Table 1) the quadratic term of elevation range was negative and highly significant but the linear term was not significant, indicating a non-linear negative response (Fig. 2D). In the analysis of extinction rates (SI Appendix Table 2) the linear term of elevation range was positive and statistically significant but the squared term was not significant, indicating a positive linear response (Fig. 2E). In the analysis of species richness (SI Appendix Table 3) the linear term was positive, the quadratic term was negative, both terms were statistically significant, and the inflection point was within the range of the data, indicating a unimodal response (Fig. 2F). These results are consistent with our theoretical predictions (Fig. 2AC). To evaluate the robustness of these patterns to other, potentially confounding factors that might influence the dynamics and structure of the study species, we constructed more complex OLS models in which we incorporated the effects of area (A), mean absolute elevation (E), mean annual rainfall (R), mean annual temperature (T), human population density (P), and broad-scale spatial effects (X coordinate, Y coordinate, X2, Y2, and XY) as additional predictors in the models. Analyses were performed using the backward variable selection method of SPSS with significance levels to enter and to leave the model being set at α = 0.01 and 0.05, respectively. The results were similar to those obtained in the simpler models, indicating a negative effect of elevation range on mean species abundance, positive effect on extinction rates, and unimodal effect on species

richness (SI Appendix Tables 1-3). Similar results were obtained in models incorporating the simultaneous effects of all 12 independent variables and in the 'best' models selected from all possible models using the AIC criterion. We also repeated the analysis replacing area with the effective surface area, which corrects for area rugosity, and obtained similar results, as expected due to the small variation of area introduced by considering rugosity (less than 5%). Although broad-scale geographic trends in the response variables were partially accounted for by the OLS models (SI Appendix Tables 1-3), the residuals of all models were positively autocorrelated (Moran's I < 0.05) at spatial scales lower than 40km (40). We therefore performed two additional sets of analyses in which all variables that were found statistically significant in the OLS models were reanalyzed for their effects on the response variables using Conditional Autoregressive Regression (CAR) and Simultaneous Autoregressive Regression (SAR) models (41). Both types of models assume that errors of the respective non-spatial models are spatially autocorrelated. As spatial autocorrelation was limited to small spatial scales (40km), the spatial models were performed using binary connectivity matrices with a distance threshold of 40km (42). All spatial analyses were performed with SAM Ver. 4.0 (Spatial Analysis in Macroecology (43)). The results obtained for the spatial models were qualitatively similar to those obtained for the non-spatial models (SI Appendix Tables 1-3). A comparison of the four regression models constructed for each response variable based on the coefficient of determination (R2) and Akaike's Information Criteria (44) adjusted for small sample size (AICC) indicated that the SAR models had the best fit to the data (highest R2 and lowest AICC) in all analyses (SI Appendix Tables 1-3).

Two additional analyses were performed to verify that the observed responses to elevation range were not artifacts of spatial or temporal variation in sampling effort. First, we repeated the original analyses using only grid cells that had a similar sampling effort for both surveys. This procedure reduced sample size by 50% (from 372 to 189) but did not affect the results of the original analyses (SI Appendix Fig. 1, SI Appendix Table 4). Second, the combined effects of elevation range (both the linear and quadratic terms) and sampling effort (160h) on species richness were tested separately for each survey period using GLM. The results indicated that sampling effort had a significantly positive effect on species richness but for both survey periods the effect of elevation range was similar to that in the original analysis (SI Appendix Fig. 2, SI Appendix Tables 5-6). We therefore feel quite confident that the unimodal pattern of species richness was not a sampling artifact.

A.4.2 Effects of niche width Further analyses were performed to test the prediction that increasing niche width shifts the level of elevation range that maximizes species richness to higher levels. In addition to its theoretical significance, this analysis was methodologically important for two reasons. First, it has often been suggested that unimodal responses of species richness to environmental factors in general and elevation in particular can arise through a simple middomain effect (45). In our system elevation range was correlated with mean elevation (r = 0.75, P < 0.001). It was therefore important to reject the possibility that the observed effect of elevation range on species richness is not a statistical artifact of this correlation. According to mid-domain theory, differences among species in niche width should affect the height of the response to the relevant gradient, but not the position of the inflection point (45). A statistically significant effect of niche width on the position of the inflection

point can therefore be interpreted as evidence against a mid-domain explanation of the unimodal pattern observed at the level of the whole community. Second, although the unimodal effect of elevation range on species richness was robust to variation in a number of potentially confounding factors (SI Appendix Tables 13), it is possible that the observed response was caused by other confounding factors for which we did not have data. A consistent and directional shift in the position of the inflection point with increasing niche width is not expected to arise from other simple mechanisms and can therefore be interpreted as a further support for the area-heterogeneity trade-off. Niche width of the study species was determined by ranking the records of each species with respect to elevation, removing outliers (highest and lowest 5% of the records), and calculating the difference between the remaining maximum and minimum elevation values. Based on the calculated niche widths we categorized the species into four 'niche quartiles' and regressed the number of species belonging to each niche quartile against elevation range and its squared value.

B. Theoretical analyses B.1 Objectives The theoretical part of the study had three main objectives. The first objective was to confirm the logic and validity of our predictions (Fig. 2A-C) by deriving the expected relationships between environmental heterogeneity, mean species abundance, extinction rates, and species richness, from the fundamental processes of birth, death, and migration of individuals. A second objective was to get a sense of the 'noise' caused by pure demographic stochasticity. Such an evaluation was important because the empirically observed data were very noisy (Fig. 2D-F). The third objective was to evaluate whether and

how community responses to environmental heterogeneity depend on the demographic rates of the modeled species, their niche width, the area, and the length of time over which the data are collected. Answering all of these questions is crucial for a proper interpretation of empirically observed heterogeneity-diversity relationships. Evaluating the magnitude of 'noise' expected in empirical tests of our predictions is also important because ecologists studying patterns of species diversity tend to ignore the inherent noise expected from demographic stochasticity and often attribute low values of R2 to the effect of ‘hidden’ variables not included in the analysis.

B.2 Description of the model Our theoretical analyses were based on a stochastic, individual based, spatially implicit model of a local community that receives immigrants from a regional species pool with N species. Individuals go through the fundamental demographic processes of reproduction, mortality and migration, all modeled with explicit rates, that are constant among species. The landscape consists of A sites, where each site can be occupied by at most one individual. Each site is conceptualized as having a unit area, so that A also indicates the total area available for the local community. Environmental heterogeneity is introduced by assigning each site an environmental value Ej that is randomly drawn from a uniform distribution between Emin and Emax. The range of E available for the local community (ER = Emax-Emin) represents the degree of environmental heterogeneity. This range is embedded within a larger range (ERR = ERmin to ERmax) representing the range of E in the area available for the regional species pool (SI Appendix Fig. 4). In this formulation, the environment is reduced into a single environmental variable representing a major niche axis (19,46) or a complex gradient such as elevation (47).

The niche of each species along the environmental gradient is described by a Gaussian function indicating the probability of establishment as a function of E (SI Appendix Fig. 4): ( E j μi )2

Pri , j Z i e

2 σi 2

,

where Pri,j(E) is the per-capita probability of establishment of species i in site j ( i {1,..., N}, j {1,..., A} ), µi is the optimum level of E for species i, and σi is a parameter representing the corresponding niche width (assumed to be equal for all species in most analyses). Values of µi are randomly selected from a uniform distribution between ERmin and ERmax. Zj is a normalization constant ensuring that all species are equivalent in their overall probability of establishment within the range ERmin – ERmax: E R max

Zj ( R E

[

e

(E j μi )2 2 σi 2

]

dE)1

min

Similar results were obtained for equal distribution of environmental conditions, created by dividing the local environmental range ERR into A equally-spaced values and assigning each to one local site. Each individual dies at an annual rate M and gives birth to one offspring at rate R. A new offspring is immediately dispersed into a random site. Immigration from the regional species pool to each site in the landscape occurs at a rate I, with equal immigration probability for each of the N species in the regional species pool. Individuals arriving at empty sites establish with probability Prij and individuals arriving in occupied sites die without being able to establish. The dynamics of the local community is modeled as a continuous-time Markov process where, in a small enough time interval, only a single demographic event (birth, death or immigration) can occur (48). This model can be considered as an individual-based, stochastic extension of Levins' model (49) in which all species are identical except for the position of their niche along an

environmental gradient. It can also be considered as an extension of the model of Kadmon and Allouche (14). While in Kadmon and Allouche's model habitats were discrete and unordered and each species could only persist in a single habitat, in the current model the environmental axis is continuous and species can persist in multiple environmental conditions with varying success. We simulated the system using Poisson thinning (50,51), implemented in C and MATLABTM (MATLAB version 7.10.0.499, The MathWorks Inc., Natick, MA)

B.3 The simulation experiments The effect of environmental heterogeneity on community characteristics was evaluated by varying the range of environmental conditions available for the local community (ER = Emax-Emin) while keeping all other parameters constant. Two sets of simulations were conducted. The first set was designed to enable evaluation of the 'noise' expected in empirical tests of our predictions due to demographic stochasticity. In these simulations each response variable was sampled only once in each simulation run (in order to ensure independence) and no averaging was applied to the results. The second set was designed to evaluate the robustness of the predicted patterns to variations in model parameters. In these simulations the results were averaged within and among different runs in order to remove the noise caused by demographic stochasticity from the observed patterns (SI Appendix Figs. 5-8). Preliminary analyses showed that all variables reach their steady state distribution after less than 5,000 years. We therefore started all simulations by running the model for 20,000 years with initial conditions of zero local abundance of all species. Similar results were obtained for simulations of a shorter time period of 50 years, with initial conditions of zero local abundance of all species, in which the community did not reach it steady state

distribution (SI Appendix Fig. 9). Sampling started from year 20,000 and differed between the two sets of simulations. In the first set (designed to generate independent realizations) the modeled community was sampled using only two snapshots, at the end of year 20,000 and at the end of year 20,001. The first snapshot was used to determine mean species abundance (average number of individuals per species) and species richness (number of species) in the local community. The second snapshot was used to determine extinction rates (number of species extinct between year 20,000 and 20,001). This procedure was repeated using 20 independent simulations for each of 15 levels of environmental range (3, 6, 9,…,45) under a fixed parameter regime (N = 500, R = 45, M = 1, I = 0.2, σ = 2, A = 1,000, ERR = 65). The qualitative results of these simulations (Fig. 2, G to I) were consistent with our empirical findings, indicating that environmental heterogeneity reduces mean species abundance, increases extinction rates, and has a unimodal effect on species richness. It is also important to note that all community responses were extremely ‘noisy’, suggesting that empirical analyses of these patterns are expected to show very low values of R2 even at the absence of any other deterministic or stochastic source of environmental variation. The only difference between the simulations and the empirical data was the functional form of the abundance-heterogeneity relationship (Fig. 2E,H) which may reflect restrictive assumptions of our model (e.g., ignorance of environmental stochasticity) or differences in the manner by which abundance was determined in the field and in our simulations. In the second set of simulations (designed to evaluate the robustness of the predicted patterns to variations in model parameters) each run lasted 120,000 years and snapshots of the community were sampled at the end of each year starting from year 20,000. The 100,001 snapshots obtained for a given run were used to calculate mean values of the three response variables (mean species abundance, extinction rates, and

species richness). This procedure was repeated using 10 independent runs for each level of environmental range under each parameter regime. We used the second type of simulations to test the manner by which the effects of environmental heterogeneity depend on reproduction rates, immigration rates, area, and temporal scale of data collection. The effect of each of these factors was tested under a wide range of values (reproduction: 30120; immigration: 0.01-0.5; area: 100-50,000, temporal scale of data collection: 1-20 years). Further simulations were performed to test the hypothesis that increasing niche width shifts the level of heterogeneity that maximizes species richness to higher levels. In these simulations the regional species pool consisted of a mixture of species representing different niche widths and the number of independent runs was 50 due to higher stochasticity. Typical results of our simulations are presented in SI Appendix Figs. 5-8. As expected from the area-heterogeneity trade-off, the general relationship between species richness and environmental heterogeneity was unimodal, but the level of heterogeneity that maximized species richness increased with increasing reproduction, immigration, area, temporal scale of data collection, and niche width. All of these factors reduce the likelihood of stochastic extinction, thereby increasing the range of environmental heterogeneity under which the relationship is positive.

C. Meta-analysis C.1 Objectives The purpose of this part of the study was to test our hypothesis that evidence of positive heterogeneity-diversity relationships obtained for insular systems might be biased due to underlying positive correlation between island area and heterogeneity. To test this hypothesis we surveyed the literature for studies providing estimates of both

environmental heterogeneity and area in insular systems, and compared the effects of environmental heterogeneity on species richness before and after accounting for the effect of area in the analysis. In contrast to most previous analyses of this sort, we included the quadratic term of environmental heterogeneity in the models. We predicted that accounting for the effect of area would significantly reduce the frequency of positive relationships and increase the frequency of non-significant and unimodal relationships in the data.

C.2 Data source To be compatible with the analyses performed in this study we searched for studies in which environmental heterogeneity was expressed by elevation range. We found 40 such studies (20,22,52-89). These studies included 66 species-island datasets, but some datasets had partial overlap in their species composition (e.g., endemic species vs. all species, woody species vs. all plant species, native vs. all species, etc.). We therefore used only a single dataset for each group of species in each study (the one that was not nested within any other dataset) for our analysis. This selection resulted in 54 different datasets.

C.3 Statistical analysis In order to test our hypothesis we first identified all datasets for which linear regression analysis revealed statistically significant effect of elevation range on species richness (a total of 43 datasets, SI Appendix Table 7). Each of these datasets was then analyzed in three steps. First, elevation range was regressed against area and the residuals of the model were determined. Second, log species richness was regressed against log area and the predicted values of log species richness were back transformed into units of species richness in order to determine the residuals of species richness over area. In the third step,

the residuals of species richness over area were regressed against the residuals of elevation range over area using both linear and quadratic regression models (note that this analysis is equivalent to the calculation of partial regression coefficients in multiple regression analysis and was necessary to account for logarithmic species-area relationships). A comparison of the results obtained for the quadratic vs. linear models using the sample size-corrected Akaike Information Criterion (AICC) enabled us to determine which of the two models better fits each dataset. The results (SI Appendix Table 7, Fig. 4A) showed that only 6 out of the 43 datasets remained significantly positive after the effect of area was adjusted for. The remaining datasets become significantly unimodal (14 patterns), non-significant (20 patterns), negative (1 pattern) or U-shaped (2 patterns). As a more formal test for the prevalence of unimodal relationships in the data we calculated the average effect size of all area-corrected quadratic models based on their Fvalues. We first transformed the F-statistic into rF, an estimate of effect size, using the formula rF F ( F df ) . The sign of the rF-values was determined by the sign of the quadratic regression coefficients. Values of rF were then transformed to Fisher's z-values for determining the mean and 95% confidence intervals of the effect size (90). The results (mean: -0.1785; confidence intervals: -0.1129 to -0.2441) indicated that the average effect size was significantly smaller than zero, supporting the hypothesis that the general relationship between species richness and elevation range was unimodal. Eighteen of the datasets we found also provided data on habitat diversity. We therefore analyzed these datasets using the same analysis performed for elevation range but with habitat diversity as a measure of environmental heterogeneity. The results were very similar to those obtained for elevation range (Fig. 4B). These overall findings support our hypothesis that previously reported evidence for positive effects of environmental heterogeneity on species richness in insular systems might be biased.

C.4 Criteria used to determine the functional relationship between species richness and the area-corrected elevation range A dataset was classified as unimodal if it satisfied the following requirements: (1) the quadratic model was statistically significant (P < 0.05), (2) the quadratic term was negative, (3) the inflection point was within the range of the data, and (4) the quadratic model showed a smaller value of the sample size-corrected Akaike Information Criterion (AICC) than the corresponding linear model. If requirements 1,3 and 4 were satisfied and the quadratic term was positive, the dataset was classified as U-shaped. Datasets were classified as positive (negative) if: (1) the dataset was not classified as unimodal, (2) the linear model was statistically significant (P < 0.05), (3) the linear term in the linear model was positive (negative) or: (1) the quadratic model was statistically significant (P < 0.05), (2) the quadratic term was negative, (3) the inflection point was to the right (left) of the range of the data or: (1) the quadratic model was statistically significant (P < 0.05), (2) the quadratic term was positive, (3) the inflection point was to the left (right) of the range of the data If the dataset was not classified as either unimodal, U-shaped, positive or negative, it was classified as non-significant.

D. References 1.

Rosenzweig ML (1995) Species diversity in space and time (University of Cambridge press, New York).

2.

Hanski I (1999) Metapopulation ecology (Oxford University Press).

3.

MacArthur RH & Wilson EO (1967) The theory of island biogeography (Princeton University Press, Princeton, NJ).

4.

Bell G (2000) The distribution of abundance in neutral communities. Am. Nat. 155(5):606-617.

5.

S. P. Hubbell, The unified neutral theory of biodiversity and biogeography. Monographs in Population Biology (Princeton University Press, Princeton, NJ, 2001).

6.

Hutchinson GE (1957) Concluding remarks, Cold Spring Harbor Symposium. Quantitative Biology 22:415-427.

7.

Chase JM & Leibold MA (2003) Ecological niches: linking classical and contemporary approaches (University of Chicago Press, Chicago, IL).

8.

Chase JM (2005) Towards a really unified theory for metacommunities. Funct. Ecol. 19(1):182-186.

9.

Haegeman B & Loreau M (2011) A mathematical synthesis of niche and neutral theories in community ecology. J. Theor. Biol. 269(1):150-165.

10.

Diamond JM (1976) Island biogeography and conservation - strategy and limitations. Science 193(4257).

11.

Simberloff D (1988) The contribution of population and community biology to conservation science. Annu. Rev. Ecol. Syst. 19:473-511.

12.

Benton TG, Vickery JA, & Wilson JD (2003) Farmland biodiversity: is habitat heterogeneity the key? Trends Ecol. Evol. 18(4):182-188.

13.

Flather CH, Hayward GD, Beissinger SR, & Stephens PA (2011) Minimum viable populations: is there a 'magic number' for conservation practitioners? Trends Ecol. Evol. 26(6):307-316.

14.

Kadmon R & Allouche O (2007) Integrating the effects of area, isolation, and habitat heterogeneity on species diversity: A unification of island biogeography and niche theory. Am. Nat. 170(3):443-454.

15.

Allouche O & Kadmon R (2009) A general framework for neutral models of community dynamics. Ecol. Lett. 12(12):1287-1297.

16.

McKinney ML (1997) Extinction vulnerability and selectivity: Combining ecological and paleontological views. Annu. Rev. Ecol. Syst. 28:495-516.

17.

Muntaner J, X. Ferrer, and A. Martı´nez-Vilalta (1984) Atles dels ocells nidificants de Catalunya i Andorra [Breeding Birds Atlas of Catalonia and Andorra] (Ketres, Barcelona, Spain).

18.

Estrada J, Pedrocchi, V., Brotons, L. & Herrando, S. (2004) Atles dels ocells nidificants de Catalunya 1999–2002 [Breeding Birds Atlas of Catalonia] (Institut Català d’Ornitologia (ICO)/Lynx Editions, Barcelona, Spain).

19.

Gravel D, Canham CD, Beaudet M, & Messier C (2006) Reconciling niche and neutrality: the continuum hypothesis. Ecol. Lett. 9(4):399-409.

20.

Ricklefs RE & Lovette IJ (1999) The roles of island area per se and habitat diversity in the species-area relationships of four Lesser Antillean faunal groups. J. Anim. Ecol. 68(6):1142-1160.

21.

Moreno Saiz JC & Lobo JM (2008) Iberian-Balearic fern regions and their explanatory variables. Plant Ecol. 198(2):149-167.

22.

Veech JA & Crist TO (2007) Habitat and climate heterogeneity maintain betadiversity of birds among landscapes within ecoregions. Glob. Ecol. Biogeogr. 16(5):650-656.

23.

Tamme R, Hiiesalu I, Laanisto L, Szava-Kovats R, & Partel M (2010) Environmental heterogeneity, species diversity and co-existence at different spatial scales. J. Veg. Sci. 21(4):796-801.

24.

Kerr JT & Packer L (1997) Habitat heterogeneity as a determinant of mammal species richness in high-energy regions. Nature 385(6613):252-254.

25.

Costa GC, Nogueira C, Machado RB, & Colli GR (2007) Squamate richness in the Brazilian Cerrado and its environmental-climatic associations. Divers. Distrib. 13(6):714-724.

26.

Huang Y, et al. (2011) Lizard species richness patterns in China and its environmental associations. Biodivers. Conserv. 20(7):1399-1414.

27.

Rahbek C & Graves GR (2001) Multiscale assessment of patterns of avian species richness. Proc. Natl. Acad. Sci. U. S. A. 98(8):4534-4539.

28.

Jimenez I, Distler T, & Jorgensen PM (2009) Estimated plant richness pattern across northwest South America provides similar support for the species-energy and spatial heterogeneity hypotheses. Ecography 32(3):433-448.

29.

Montigny MK & MacLean DA (2005) Using heterogeneity and representation of ecosite criteria to select forest reserves in an intensively managed industrial forest. Biol. Conserv. 125(2):237-248.

30.

Kati V, et al. (2010) Towards the use of ecological heterogeneity to design reserve networks: a case study from Dadia National Park, Greece. Biodivers. Conserv. 19(6):1585-1597.

31.

Fuhlendorf SD, et al. (2006) Should heterogeneity be the basis for conservation? Grassland bird response to fire and grazing. Ecol. Appl. 16(5):1706-1716.

32.

Tilman D (2004) Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community assembly. Proc. Natl. Acad. Sci. U. S. A. 101(30):10854-10861.

33.

Ruokolainen L, Ranta E, Kaitala V, & Fowler MS (2009) When can we distinguish between neutral and non-neutral processes in community dynamics under ecological drift? Ecol. Lett. 12(9):909-919.

34.

Wasser, D. E. & Sherman, P. W. Avian longevities and their interpretation under evolutionary theories of senescence. Journal of Zoology 280, 103-155 (2010).

35.

Carnicer, J., Brotons, L., Sol, D. & Jordano, P. Community-based processes behind species richness gradients: contrasting abundance-extinction dynamics and sampling effects in areas of low and high productivity. Glob. Ecol. Biogeogr. 16, 709-719 (2007).

36.

Moreno-Rueda, G. & Pizarro, M. Relative influence of habitat heterogeneity, climate, human disturbance, and spatial structure on vertebrate species richness in Spain. Ecol. Res. 24, 335-344 (2009).

37.

Sutherland, G. D., Harestad, A. S., Price, K. & Lertzman, K. P. Scaling of Natal dispersal distance in terrestrial birds and mammals. Conservation Ecology (online) 4, 16 (2000).

38.

Kunin, W. E. Biodiversity at the edge: A test, of the importance of spatial "mass effects" in the Rothamsted Park Grass experiments. Proc. Natl. Acad. Sci. U. S. A. 95, 207-212 (1998).

39.

Ferrer, X., Carrascal, L. M., Gordo, O. & Pino, J. Bias in avian sampling effort due to human preferences: An analysis with Catalonian birds (1900-2002). Ardeola 53, 213-227 (2006).

40.

Legendre, P. & Legendre, L. Numerical Ecology. (Elsevier Science B.V., 1998).

41.

Dormann, C. F. et al. Methods to account for spatial autocorrelation in the analysis of species distributional data: a review. Ecography 30, 609-628 (2007).

42.

Sokal, R. R. Testing statistical significance of geographic-variation patterns. Syst. Zool. 28, 227-232 (1979).

43.

Rangel, T. F., Diniz-Filho, J. A. F. & Bini, L. M. SAM: a comprehensive application for Spatial Analysis in Macroecology. Ecography 33, 46-50 (2010).

44.

Burnham, K. P. & D.R.Anderson. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. (Springer-Verlag, 2002).

45.

McCain, C. M. Global analysis of bird elevational diversity. Glob. Ecol. Biogeogr. 18, 346-360 (2009).

46.

Pulliam, H. R. On the relationship between niche and distribution. Ecol. Lett. 3, 349-361 (2000).

47.

Whittaker, R. H. Vegetation of the Great Smoky Mountains. Ecol. Monogr. 26, 180 (1956).

48.

Allouche, O. & Kadmon, R. A general framework for neutral models of community dynamics. Ecol. Lett. 12, 1287-1297 (2009).

49.

Levins, R. Some demographic and genetic consequences of environmental heterogeneity for biological control. Bull. Entomol. Soc. Am. 15, 227-240 (1969).

50.

Durrett, R. & Levin, S. The importance of being discrete (and spatial). Theor. Popul. Biol. 46, 363-394 (1994).

51.

Berec, L. Techniques of spatially explicit individual-based models: construction, simulation, and mean-field analysis. Ecol. Model. 150, 55-81 (2002).

52.

Abbott, I. Numbers of plant, insect and land bird species on nineteen remote islands in southern-hemisphere. Biol. J. Linn. Soc. 6, 143-152 (1974).

53.

Ackerman, J. D., Trejo-Torres, J. C. & Crespo-Chuy, Y. Orchids of the West Indies: predictability of diversity and endemism. J. Biogeogr. 34, 779-786 (2007).

54.

Cardoso, P., Arnedo, M. A., Triantis, K. A. & Borges, P. A. V. Drivers of diversity in Macaronesian spiders and the role of species extinctions. J. Biogeogr. 37, 10341046 (2010).

55.

Case, T. J. Species numbers, density compensation, and colonizing ability of lizards on islands in gulf of california. Ecology 56, 3-18 (1975).

56.

Case, T. J. & Bolger, D. T. The role of introduced species in shaping the distribution and abundance of island reptiles. Evol. Ecol. 5, 272-290 (1991).

57.

Chen, Y.-H. Combining the Species-Area-Habitat Relationship and Environmental Cluster Analysis to Set Conservation Priorities: a Study in the Zhoushan Archipelago, China. Conserv. Biol. 23, 537-545 (2009).

58.

Duarte, M. C., Rego, F., Romeiras, M. M. & Moreira, I. Plant species richness in the Cape Verde islands - eco-geographical determinants. Biodivers. Conserv. 17, 453-466 (2008).

59.

Dueser, R. D. & Brown, W. C. Ecological correlates of insular rodent diversity. Ecology 61, 50-56 (1980).

60.

Fattorini, S. Biogeography of tenebrionid beetles (Coleoptera: Tenebrionidae) in the circum-Sicilian islands (Italy, Sicily): Multiple biogeographical patterns require multiple explanations. Eur. J. Entomol. 108, 659-672 (2002).

61.

Hamilton, T. H., Rubinoff, I., Barth, R. H., Jr. & Bush, G. L. Species Abundance: Natural Regulation of Insular Variation. Science 142, 1575-1577 (1963).

62.

Hannus, J.-J. & von Numers, M. Vascular plant species richness in relation to habitat diversity and island area in the Finnish Archipelago. J. Biogeogr. 35, 10771086 (2008).

63.

Harris, M. P. Galapagos avifauna. Condor 75, 265-278 (1973).

64.

Johnson, M. P., Simberloff D. S. Environmental Determinants of Island Species Numbers in the British Isles. J. Biogeogr. 1, 149-154 (1974).

65.

Johnson, M. P., Mason, L. G. & Raven, P. H. Ecological parameters and plant species diversity. Am. Nat. 102, 297-& (1968).

66.

Johnson, M. P. & Raven, P. H. Species number and endemism - galapagos archipelago revisited. Science 179, 893-895 (1973).

67.

Kratter, A. W. Montane avian biogeography in southern california and bajacalifornia. J. Biogeogr. 19, 269-283 (1992).

68.

Losos, J. B. Island biogeography of day geckos (phelsuma) in the indian-ocean. Oecologia 68, 338-343 (1986).

69.

McMaster, R. T. Factors influencing vascular plant diversity on 22 islands off the coast of eastern North America. J. Biogeogr. 32, 475-492 (2005).

70.

Millien-Parra, V. & Jaeger, J. J. Island biogeography of the Japanese terrestrial mammal assemblages: an example of a relict fauna. J. Biogeogr. 26, 959-972 (1999).

71.

Nakamura, K., Suwa, R., Denda, T. & Yokota, M. Geohistorical and current environmental influences on floristic differentiation in the Ryukyu Archipelago, Japan. J. Biogeogr. 36, 919-928 (2009).

72.

Newmark, W. D. Species area relationship and its determinants for mammals in western north-american national-parks. Biol. J. Linn. Soc. 28, 83-98 (1986).

73.

Panitsa, M., Trigas, P., Iatrou, G. & Sfenthourakis, S. Factors affecting plant species richness and endemism on land-bridge islands - An example from the East Aegean archipelago. Acta Oecologica-International Journal of Ecology 36, 431437 (2010).

74.

Panitsa, M., Tzanoudakis, D., Triantis, K. A. & Sfenthourakis, S. Patterns of species richness on very small islands: the plants of the Aegean archipelago. J. Biogeogr. 33, 1223-1234 (2006).

75.

Perry, G., Rodda, G. H., Fritts, T. H. & Sharp, T. R. The lizard fauna of Guam's fringing islets: island biogeography, phylogenetic history, and conservation implications. Glob. Ecol. Biogeogr. 7, 353-365 (1998).

76.

Power, D. M. Numbers of bird species on california islands. Evolution 26, 451-& (1972).

77.

Power, D. M. Avifauna richness on california channel islands. Condor 78, 394-398 (1976).

78.

Price, J. P. Floristic biogeography of the Hawaiian Islands: influences of area, environment and paleogeography. J. Biogeogr. 31, 487-500 (2004).

79.

Reed, T. The number of breeding land-bird species on british islands. J. Anim. Ecol. 50, 613-624 (1981).

80.

Rogers, G. & Overton, J. Regional patterns of plant species richness in southern New Zealand. N. Z. J. Bot. 38, 609-627 (2000).

81.

Sfenthourakis, S. The species-area relationship of terrestrial isopods (Isopoda; Oniscidea) from the Aegean archipelago (Greece): A comparative study. Global Ecology and Biogeography Letters 5, 149-157 (1996).

82.

Sfenthourakis, S. & Triantis, K. A. Habitat diversity, ecological requirements of species and the Small Island Effect. Divers. Distrib. 15, 131-140 (2009).

83.

Vanderwerff, H. Species number, area and habitat diversity in the galapagosislands. Vegetatio 54, 167-175 (1983).

84.

Vuilleumier, F. Insular biogeography in continental regions .1. northern andes of south-america. Am. Nat. 104, 373-& (1970).

85.

Weissman, D. B. & Rentz, D. C. Zoogeography of the Grasshoppers and Their Relatives (Orthoptera) on the California Channel Islands. J. Biogeogr. 3, 105-114 (1976).

86.

Welter-Schultes, F. W. & Williams, M. R. History, island area and habitat availability determine land snail species richness of Aegean islands. J. Biogeogr. 26, 239-249 (1999).

87.

White, P. S., Miller, R. I. & Ramseur, G. S. The species-area relationship of the southern Appalachian high peaks - vascular plant richness and rare plantdistributions. Castanea 49, 47-61 (1984).

88.

Wilcox, B. A., Murphy, D. D., Ehrlich, P. R. & Austin, G. T. Insular biogeography of the montane butterfly faunas in the great-basin - comparison with birds and mammals. Oecologia 69, 188-194 (1986).

89.

Yeakley, J. A. & Weishampel, J. F. Multiple source pools and dispersal barriers for Galapagos plant species distribution. Ecology 81, 893-898 (2000).

90.

Borenstein, M., Hedges, L. V., Higgins, J. P. T. & Rothstein, H. R. Introduction to Meta-Analysis. (John Wiley & Sons, 2009).

SI APPENDIX TABLES SI Appendix Table 1. Results of regression models testing the effects of elevation range (H), the squared value of elevation range (H2), area (A), geographical coordinates (X, Y, X2, Y2, XY), mean absolute elevation (E), mean annual temperature (T), mean annual rainfall (R), and human population size (P) on mean species abundance of breeding birds in Catalonia. Analyses were performed using four different models: OLS incorporating the simultaneous effects of the linear term and squared term of elevation range (OLS – basic), OLS incorporating all independent variables with backward variable selection (OLS – extended), Conditional Autoregressive models (CAR) incorporating those variables that were found statistically significant in the OLS models, and corresponding Simultaneous Autoregressive models (SAR). For each model we report the t-values of each variable, the respective significance level (bold values indicate significant effects), the R2 value of the model, the AICC value, and the difference between AICC value of the relevant model and the 'best' model (the one with the lowest value).

H H2 A X Y X2 Y2 XY E R T

OLS (Basic) OLS (Extended) 0.69 -----4.13***

CAR -----

SAR -----

-4.30*** 4.70*** -----4.65*** -7.83*** ----7.66***

-4.57*** 3.53*** -----1.15 3.41*** ----3.27***

-3.91*** 2.07 ----0.29 -1.11 ----0.95

-9.83*** ---------

-10.30*** ---------

-10.89*** ---------

P -4.05*** -4.15*** -4.24*** ----------------------------------------------------------------------------------------------R2 0.32 0.59 0.54 0.63 AICC -1668.6 -1848 -1802.7 -1880.6 ΔAICC 212 32.6 77.9 0 ----------------------------------------------------------------------------------------------* P < 0.05, ** P < 0.01, *** P < 0.001

SI Appendix Table 2. Results of regression models testing the effects of elevation range (H), the squared value of elevation range (H2), area (A), geographical coordinates (X, Y, X2, Y2, XY), mean absolute elevation (E), mean annual temperature (T), mean annual rainfall (R), and human population size (P) on extinction rates of breeding birds in Catalonia. Analyses were performed using four different models: OLS incorporating the simultaneous effects of the linear term and squared term of elevation range (OLS – basic), OLS incorporating all independent variables with backward variable selection (OLS – extended), Conditional Autoregressive models incorporating those variables that were found statistically significant in the OLS models (CAR), and corresponding Simultaneous Autoregressive models (SAR). For each model we report the t-values of each variable, the respective significance level (bold values indicate significant effects), the R2 value of the model, the AICC value, and the difference between AICC values of the relevant model and the 'best' model (the one with the lowest value).

H H2 A

OLS (Basic) OLS (Extended) 3.08** 4.63*** -1.54 -----4.09***

CAR 4.64*** -----2.12*

SAR 3.54*** -----1.36

X

-----

-----

-----

Y X2 Y2 XY E R

7.38*** ---------------------

3.10** ---------------------

3.81*** ---------------------

T 4.56*** 2.47* 2.21* P -------------------------------------------------------------------------------------------------------------------R2 0.32 0.23 0.31 0.40 AICC 2587.4 2528.9 2493.4 2442.0 ΔAICC 145.4 86.9 51.4 0 ----------------------------------------------------------------------------------------------* P < 0.05, ** P < 0.01, *** P < 0.001

SI Appendix Table 3. Results of regression models testing the effects of elevation range (H), the squared value of elevation range (H2), area (A), geographical coordinates (X, Y, X2, Y2, XY), mean absolute elevation (E), mean annual temperature (T), mean annual rainfall (R), and human population size (P) on species richness of breeding birds in Catalonia. Analyses were performed using four different models: OLS incorporating the simultaneous effects of the linear term and squared term of elevation range (OLS – basic), OLS incorporating all independent variables with backward variable selection (OLS – extended), Conditional Autoregressive models incorporating those variables that were found statistically significant in the OLS models (CAR), and corresponding Simultaneous Autoregressive models (SAR). For each model we report the t-values of each variable, the respective significance level (bold values indicate significant effects), the R2 value of the model, the AICC value, and the difference between AICC values of the relevant model and the 'best' model (the one with the lowest value).

H H2 A X

OLS (Basic) OLS (Extended) 6.27*** 5.89*** -5.07*** -3.63*** 7.59*** -----

Y X2 Y2 XY

5.70*** ---------3.18**

CAR 4.86*** -2.63** 6.77*** -----

SAR 5.08*** -3.13*** 6.86*** -----

3.06** ---------2.09**

0.53 ---------1.43

E -4.53*** -3.24*** -3.26*** R ------------T ------------P -------------------------------------------------------------------------------------------------------------------R2

0.13

0.31

0.35

0.36

AICC 2985.2 2904.6 2884.3 2882.6 ΔAICC 102.6 22.0 1.7 0 ----------------------------------------------------------------------------------------------* P < 0.05, ** P < 0.01, *** P < 0.001

SI Appendix Table 4. Results of regression models testing the effects of elevation range (H) and its squared value (H2) on mean species abundance, extinction rates, and species richness, based on 189 grid cells for which sampling effort was similar in the two survey periods (1983 and 2002). Values reported are t-values of each variable and their significance levels (bold values indicate significant effects).

H H2

Mean abundance 0.616 -3.430***

Extinction rates 2.647** -1.708

Species richness 5.452*** -4.735***

----------------------------------------------------------------------------------------------* P < 0.05, ** P < 0.01, *** P < 0.001

SI Appendix Table 5. Results of GLM testing the effects of elevation range (H) and its squared value (H2) on species richness in the first survey (1983) with (left) and without (right) including the effect of sampling effort (SE, five categories based on the number of sampling hours) in the model. Note that sampling effort had a positive effect on species richness but did not affect the qualitative results obtained for the effect of elevation range. See SI Appendix Fig. 2 for a graphical presentation of the relationship between species richness and elevation range.

H H2 SE < 40h 40h < SE < 80h 80h < SE < 120h 120h < SE < 160h SE > 160h

Without correction Coefficient t 43.8 7.21*** -16.8 -5.31***

With correction Coefficient t 29.2 5.66*** -11.66 -4.35*** -36.6 -4.12*** -34.7 -3.95**** -19.7 -2.22* -11.8 -1.32 (Set to zero)

SI Appendix Table 6. Results of GLM testing the effects of elevation range (H) and its squared value (H2) on species richness in the second survey with (left) and without (right) including the effect of sampling effort (SE, five categories based on the number of sampling hours) in the model. Note that sampling effort had a positive effect on species richness but did not affect the qualitative results obtained for the effect of elevation range. See SI Appendix Fig. 2 for a graphical presentation of the relationship between species richness and elevation range.

H H2 SE < 40h 40h < SE < 80h 80h < SE < 120h 120h < SE < 160h SE > 160h

Without correction Coefficient t

With correction Coefficient t

23.8 -11.7

14.2 -6.4 -31.1 -26.3 -16.5 -10.0 (Set to zero)

3.88*** -3.65***

2.70** -2.34* -8.50*** -7.41*** -4.75*** -2.73**

SI Appendix Table 7. Results of linear and quadratic regression analyses for the area-corrected effect of elevation range on species richness in 43 datasets of insular systems. In both models the response variable is the residuals of species richness over area and the independent variable is the residuals of elevation range over area (including its squared term in the quadratic model). The model that better fits the data was determined as the one with the smaller value of the sample size-corrected Akaike Information Criterion (AICC) provided that the pattern was statistically significant (see SI Methods section C.4 for full details).

Linear model

Quadratic model

Area Group Birds

R2

[hectare]

AICC

R2

P-value

AICC

P-value Pattern

Source

9,000-151,000

0.125

-57.3

0.0758

0.323

-60.8

0.0172

Unimodal

Ricklefs & Lovette 1999

9,000-151,000

0.229

-59.7

0.0221

0.322

-60.8

0.0175

Unimodal

Ricklefs & Lovette 1999

Butterflies

9,000-151,000

0.447

-66.4

0.0038

0.428

-64.7

0.0139

Positive

Ricklefs & Lovette 1999

Birds

0.8-10,692

0.0758

-318

0.0105

0.141

-322

0.0018

Unimodal

Reed 1981

Plants

3,500-99,100

0.37

-26.3

0.0367

0.76

-34.1

0.0028

Unimodal

Duarte et al. 2008

Spiders

300-205,800

0.172

-62.5

0.0394

0.184

-61.5

0.069

Positive

Cardoso et al. 2010

Birds

215-46,870

-0.0769

-38.3

0.995

-0.166

-35.6

0.998

NS

Chen 2009

0.05-68.3

0.0147

-362

0.141

0.0044

-360

0.313

NS

Hannus & von Numers 2008

Reptiles and amphibians

Plants, Brunskar archipelago

Plants, Getskar

0.004-14.314

0.0248

-341

0.0897

0.0599

-343

0.0368

Unimodal

Hannus & von Numers 2008

Mammals

1.22-30.51

0.0329

-51.5

0.227

0.0070

-49.6

0.371

NS

Wilcox et al. 1986

Butterflies

1.22-30.51

-0.0351

-50.2

0.524

-0.066

-48.3

0.632

NS

Wilcox et al. 1986

Birds

1.22-30.51

0.525

-41.8

0.0030

0.478

-39

0.0156

Positive

Wilcox et al. 1986

Birds

2,500-348,700

-0.075

-38.3

0.88

-0.119

-36.2

0.779

NS

Vuilleumier 1970

Plants

0.05-50

0.115

-393

0.0008

0.19

-399

5.78E-05

Unimodal

Panitsa et al. 2006

0.0942

-77.5

0.079

0.109

-76.6

0.115

NS

Newmark 1986

archipelago

Mammals

8,3002,073,600

Plants

189-274,063

0.234

-239

9.31E-05

0.239

-239

0.0003

Positive

Rogers and Overton 2000

Terrestrial isopods

0.23-826,120

0.12

-415

0.0005

0.168

-419

0.0001

U-shaped

Sfenthourakis & Triantis 2009

0.225

-132

0.0024

0.204

-130

0.0098

Positive

Veech & Crist 2007

0.556

-218

1.12E-09

0.605

-222

5.11E-10

U-shaped

Ackerman et al. 2007

Birds

Orchids

960,00075,290,000 2,10010,500,700

Lizards

60-143,000

-0.0047

-75

0.355

0.177

-78.6

0.0497

Unimodal

Case 1975

Plants

60-143,000

0.109

-77.9

0.064

0.702

-103

1.17E-06

Unimodal

Case 1975

Rodents

29-2,197

0.337

-22.1

0.0593

0.23

-18.7

0.193

NS

Dueser & Brown 1980

Land birds

259-582,488

0.139

-41.7

0.0944

0.287

-43

0.052

NS

Harris 1973

Plants

259-582,488

-0.0769

-38.3

0.995

0.216

-41.6

0.0918

NS

Harris 1973

Plants

0.5-229,850

-0.0236

-155

0.819

-0.0385

-153

0.788

NS

Johnson & Simberloff 1974

Day geckos

0.4-552,500

0.0781

-189

0.0305

0.192

-194

0.0031

Unimodal

Losos 1986

Plants

52-34,706

0.298

-48.8

0.0169

0.281

-47

0.0462

Negative

Power 1972

Birds

2,600-24,900

0.214

-17.1

0.139

0.209

-14.8

0.24

NS

Power 1976

Plants

2,600-24,900

0.335

-18.5

0.0775

0.212

-14.8

0.237

NS

Power 1976

Plants

52-582,488

-0.0621

-46

0.804

0.299

-51.6

0.0325

Unimodal

Hamilton 1963

Orthoptera

2.6-249

-0.145

-14.1

0.747

-0.368

-10.4

0.943

NS

Weissman & Rentz 1976

Plants

1-466,930

-0.0461

-70

0.863

0.34

-79.3

0.0060

Unimodal

Johnson & Raven 1973

Plants

3-26,668

-0.0498

-65.8

0.953

-0.0893

-63.7

0.871

NS

McMaster 2005

Plants

1-1,636

0.157

-62.2

0.0473

0.144

-60.6

0.103

Positive

Panitsa et al. 2010

Snails

0.58-826,000

-0.0149

-269

0.806

0.0222

-271

0.186

NS

Welter-Schultes & Williams 1999

Plants

1-466,930

-0.0072

-96.4

0.379

0.333

-107

0.0020

Unimodal

Yeakley & Weishampel 2000

Tenebrionid beetles

380-826,000

-0.0312

-109

0.806

-0.0537

-107

0.812

NS

Fattorini 2002

Plants

5-34,706

0.268

-40.1

0.0337

0.2050

-37.4

0.1130

NS

Johnson et al. 1968

Mammals

400-23,051,000

-0.0484

-8.81

0.43

-0.295

-4.26

0.685

NS

Millien-Parra & Jaeger 1999

Plants

2-1,043,300

0.176

-54.4

0.0468

0.272

-55.2

0.0361

Unimodal

Price 2004

Terrestrial Isopods

3-47,620

-0.0239

-160

0.885

-0.0212

-159

0.574

NS

Sfenthourakis 1996

Plants

23-467,000

-0.026

-50.4

0.462

0.345

-57.1

0.0163

Unimodal

Vanderwerff 1983

Plants

3-4,900

0.0251

-22

0.299

-0.0015

-19.8

0.417

NS

White 1984

Extinction rates

Mean abundance

Species richness

SI APPENDIX FIGURES

Environmental heterogeneity SI Appendix Figure 1. Effect of elevation range on mean species abundance, extinction rates, and species richness, in 189 grid cells that had similar sampling efforts in the two survey periods (1983 and 2002). All patterns show the same trends as in the analysis of the full dataset (Fig. 2) and in spite of a 50% reduction in sample size are highly significant (SI Appendix Table 4)

Species richness

• 1983 • 2002

Environmental heterogeneity

SI Appendix Figure 2. Effect of elevation range on species richness in the first (1983) and second (2002) surveys based on the complete dataset (372 grid cells). For both surveys the relationship was significantly unimodal both before and after correction for spatial variation in sampling effort (SI Appendix Tables 5-6).

SI Appendix Figure 3. Effect of environmental heterogeneity on extinction probability of species with different niche widths. Data refer to breeding birds in Catalonia with elevation range as a measure of heterogeneity. Each color represents a different quartile of niche width. Increasing niche width decreases the slope of the regression line. P-values indicate significance levels of the corresponding linear regression models.

Environmental range (regional)

Environmental range (local)

SI Appendix Figure 4. Schematic presentation of species responses to variation in environmental variation in our model. The environment is reduced to a single variable E and each species i is characterized by a Gaussian response with optimum µi and standard deviation ('niche width') σi indicating its establishment probability as a function of E. The range of E available for the local community (ER = Emax-Emin) is embedded within the regional range ERR = ERmax-ERmin.

Species richness Species abundance Extinction rates

R = 30 R = 45 R = 90

Environmental heterogeneity (ER)

SI Appendix Figure 5. Effect of reproduction rates on the relationship between environmental heterogeneity and mean species abundance, extinction rates, and species richness, based on the stochastic model of community dynamics. Each dot represents the mean of 10 independent runs of 100,001 years. Parameters: N = 500, A = 1,000, ERR = 65, σ = 2, M = 1, I = 0.2.

Species richness Species abundance Extinction rates

I = 0.05 I = 0.20 I = 0.50

Environmental heterogeneity (ER)

SI Appendix Figure 6. Effect of immigration rates on the relationship between environmental heterogeneity and mean species abundance, extinction rates, and species richness, based on the stochastic model of community dynamics. Each dot represents the mean of 10 independent runs of 100,001 years. Parameters: N = 500, A = 1,000, ERR = 65, σ = 2, M = 1, R = 45.

Log Species richness

A A = 100 A = 1000 A =50000

Log Species richness

B TE = 1 TE = 5 TE = 20

Environmental heterogeneity (ER) SI Appendix Figure 7. Effect of area (A = number of sites) and temporal extent of data collection (TE = number of years over which the data are recorded) on the relationship between environmental heterogeneity and species richness, based on the stochastic model of community dynamics. Note that extending the spatial and temporal scale of data collection may turn a unimodal pattern into a positive one. Each dot represents the mean of 10 independent runs of 100,001 years. Parameters: N = 500, A = 1,000, R = 45, M = 1, I = 0.2, σ = 2, ERR = 65.

Niche width = 4

I = 0.01 I = 0.05 I = 0.10

Niche width = 12

Niche width = 30

Environmental heterogeneity (ER)

Environmental heterogeneity (ER)

Log Species richness

Log Species richness

Niche width = 1

SI Appendix Figure 8. Effect of niche width (σ) and immigration (I) on the relationship between environmental heterogeneity and species richness based on the stochastic model of community dynamics. Note that increasing niche width may turn a predominantly negative pattern into a positive one with unimodal patterns occurring at intermediate niche widths. Each dot represents the mean of 50 independent runs of 100,001 years. Parameters: N = 1,500, A = 5,000, R = 100, M = 1, ERR = 150.

Species richness Species abundance Extinction rates Environmental heterogeneity (ER)

SI Appendix Figure 9. Simulated responses of ecological communities to environmental heterogeneity after different time intervals. Samples taken after 20,000 years (as in the manuscript, allowing system to reach steady state, blue) and 50 years (orange). Each dot represents the average of 10 independent simulation runs for the 20,000 years period, and 2500 independent runs for the 50 years period. Parameter values are identical to those used in Fig. 3 of the manuscript: A (area) = 1000, R (reproduction) = 45, M (mortality) = 1, I (immigration) = 0.2, and σ (niche width) = 2.