Evaluating the correlation between area, environmental heterogeneity, and species richness using terrestrial isopods (Oniscidea) from the Pontine Islands (West Mediterranean)

Area and environmental heterogeneity influence species richness in islands. Whether area or environmental heterogeneity is more relevant in determining species richness is a central issue in island biogeography. Several models have been proposed, addressing the issue, and they can be reconducted to three main hypotheses developed to explain the species-area relationship: (1) the area-per se hypothesis (known also as the extinction-colonisation equilibrium), (2) the random placement (passive sampling), and the (3) environmental heterogeneity (habitat diversity). In this paper, considering also the possible influence of geographic distance on island species richness, we explore the correlation between area, environmental heterogeneity, and species richness by using faunistic data of Oniscidea inhabiting the Pontine Islands, a group of five small volcanic islands and several islets in the Tyrrhenian Sea, located about 60 km from the Italian mainland. We found that the colonisation of large Pontine Islands may occur via processes independent of geographic distance which could instead be an important factor at a much smaller scale. Such processes may be driven by a combination of anthropogenic influences and natural events. Even in very small-size island systems, environmental heterogeneity mostly contributes to species richness. Environmental heterogeneity could influence the taxocenosis structure and, ultimately, the number of species of Oniscidea via direct and indirect effects, these last mediated by area which may or may not have a direct effect on species richness.


Introduction
The positive relationship between species diversity and area is certainly one of the most solid paradigms in island biogeography. Over the years, biogeographers have been addressing such an issue in detail, investigating what ecological factors might ultimately contribute to the species-area relationship. However, environmental heterogeneity also influences species richness in islands (Hortal et al., 2009) and whether area or environmental heterogeneity is more relevant in determining species richness is also a central issue in island biogeography. Despite that the expression "habitat diversity" has been and is still being widely used to refer to "environmental heterogeneity," the use of the first expression received substantial criticism, well summarized and further addressed in Fattorini et al. (2015). We agree with such criticism, and in particular, we argue that whereas "habitat" is properly used when referring to the resources needed by a single species, or set of species requiring similar resources, it is improperly used when referring to general environmental features used as proxies of ecological requirements of species.
Several models have been proposed, addressing the issue, and they can be reconducted to three main hypotheses developed to explain the species-area relationship: (1) the areaper se hypothesis (known also as the extinction-colonisation equilibrium), (2) the random placement (passive sampling), and the (3) environmental heterogeneity (habitat diversity). The first predicts a higher probability for colonists to colonise larger islands than small ones. The possibility of setting larger populations would limit extinction. The second hypothesis is based on the assumption that larger islands host a larger number of individuals, which would increase the probability of observing more species. The third hypothesis is based on the observation that larger islands are generally more structured, with the occurrence of different types of habitat that would ultimately determine higher diversity. In this regard, it has been theorised that area and environmental heterogeneity may multiplicatively interact to predict species richness better than by area alone (Triantis et al., 2003). Similarly, studies based on path analysis, a method that implies causal relationships, have shown that the influence of area on species richness has a direct component and an indirect component through the effect of area on environmental heterogeneity (Kohn & Walsh, 1994;Power, 1972;Triantis et al., 2006).
At small spatial scales, environmental heterogeneity and area can become uncorrelated, and area alone is no longer an efficient proxy of the combined effects of environmental heterogeneity and island size on species richness (Triantis et al., 2003(Triantis et al., , 2005. Indeed, below a certain threshold value, the number of species can vary independently of area generating the "Small Island Effect" (SIE). Whether and to which extent environmental heterogeneity can still predict species richness when area alone cannot is ground for further debate and investigation.
Oniscidea are particularly suitable for biogeographical and ecological studies because of their low dispersal ability (Vandel, 1960) and high stenoecy (Argano & Manicastri, 1996), which may favour strong genetic differentiation (Gentile & Sbordoni, 1998;Gentile et al., 2010). Past studies of Oniscidea from Mediterranean islands have indeed indicated that the number of species is directly proportional to environmental heterogeneity, which may also influence community structure (Cazzolla Gatti et al., 2018;Sfenthourakis, 1996a). Additionally, in a refined analysis that accounted for the fact that species exploit different ranges of habitats (Hart & Horwitz, 1991;Ricklefs & Lovett, 1999), Sfenthourakis and Triantis (2009) investigated the relation area/environmental heterogeneity/species richness in the frame of SIE. They provided evidence of a strong effect of "habitat availability" on the SIE, indicating that generalists dominate islands within the SIE range while specialists contribute more in larger islands. Interestingly, in their investigation, they also used path analysis to address SIE following an approach conceptually different from others. Such an approach was criticized by Dengler (2010; but see Triantis & Sfenthourakis, 2012 for a response) mainly because (i) it did not include islands with no species when investigating SIE, (ii) it relied on an arbitrary definition of "habitat," not transferable among island datasets, and (iii) it would not account for the proportion of "habitats" on the islands.
Whereas more recent papers would suggest that inclusion of islands with no species may have little effect on the evidence for a SIE, mostly evident at a local scale (Morrison, 2014;Schrader, 2019), little can be done to improve the transferability of "habitats" among island datasets, especially considering, as said earlier, that "habitat" is a taxon-dependent concept, as it refers to the resources needed by a single species, or set of species requiring similar resources. However, keeping in mind the limitations of the use of biotope types as proxies of species' requirements, it could be possible to incorporate information on the abundance of island biotope types into the model proposed by Triantis et al. (2006).
Recently, it has also been suggested that the occurrence of a SIE should be investigated for robustness using a wide range of piece-wise models (Gao et al., 2019). Whereas such an approach is sound, aiming at introducing further rigorousness in the detection of SIE, it still does not allow the evaluation of SIE in the light of environmental heterogeneity. For this reason, it seems that path analysis, aware of possible limitations, may deserve further test, especially at a very small geographic scale.
In this paper, we use path analysis as well as other statistical approaches, to analyse faunistic data in the light of their possible dependence of area, geographic distance, and environmental heterogeneity. We used data of Oniscidea inhabiting the Pontine Islands, a group of five small volcanic islands and several islets in the Tyrrhenian Sea, located about 60 km from the Italian mainland ( Fig. 1). Additionally, for comparison purposes, we also used path analysis to address the relationship among species richness, area, and environmental heterogeneity as well as the existence of SIE using data of terrestrial isopods from the Sicilian islands, recently published in Cazzolla Gatti et al. (2018), who did not apply such an approach. We used faunal and environmental data, these last considered both in terms of presence/ absence of biotope types and their abundance to incorporate information of abundance of island biotope types into the path analysis model. The limited number of islands composing the Pontine archipelago and the variance of species richness in the islands imposed some limits to our analytical approach. For these reasons, we did not conduct a refined generalist/ specialist analysis as in Sfenthourakis and Triantis (2009). Nevertheless, this work may still provide further insights on the role of geographic distance, area, and environmental heterogeneity as predictors of the species richness and taxocenosis structure of organisms with limited dispersal ability, inhabiting a continental archipelago composed only by islands of very small size.

Methods
Oniscidea data used for statistical analyses were from Gentile and Argano (2005), as revised in Gentile et al. (2019). Further data are shown in the Supplementary Material (Table S1).
Aware that the categorization of environmental heterogeneity is problematic and not univocal, we estimated H in two ways: as the number of island biotope types (B and LogB) and accounting for their abundance in terms of integer dam 2 . For this last purpose, two different statistics have been estimated: the Shannon index (Shannon & Weaver, 1949) and 1-D, where D is dominance as defined by Simpson (1949). Necessary ecological data were acquired from maps published by Veri et al. (1979). These maps were used because they are proximate to the time during which the samples were collected. Thirteen different types of biotope were considered, reported in the Supplementary Material (Table S2). Their abundance in the different islands and islets is reported in Table S3. The image analysis and coverage quantification have been performed by using IMAGEJ (Schneider et al., 2012). This software is freely available at the website https:// imagej. nih. gov/ ij/. We investigated S/A/H relationships considering these variables in their logarithmic transformation, except when H was estimated by the Shannon and 1-D indexes.
We used standard multiple regression to focus on how much area and environmental heterogeneity uniquely contributed to a single linear equation that predicts species richness. In this analysis, we verified that the assumption of normality of the residuals was met. Partial correlation coefficients (pc) were used to measure the degree of association between two random variables, controlling for the third. Tolerance values were estimated as a measure of collinearity between variables. The threshold of tolerance beyond which the resulting estimates (regression coefficients) were considered increasingly unreliable was set at 0.1. Statistical significance of t statistic was estimated at a one-tailed test because only a positive relationship between S and A/H has biological meaning. In fact, a negative relationship is comprised in the null hypothesis of no effect of A/H on S.
Additionally, following an approach similar to Triantis et al. (2006) we performed a path analysis to seek support to the a priori model according to which island area directly affects environmental heterogeneity and both area and environmental heterogeneity directly affect species number per island. Path analysis was run multiple times using each time a different estimator of H according to a model shown in Fig. 2. Collinearity between A, H, and S was measured by the variance inflation factor (VIF), with VIF < 10 considered as acceptable. In separate path analyses, we also applied the rationale as per Triantis et al. (2006) to estimate the upper limit of the SIE and calculate the net effects of area. This was done by successive elimination of data, removing islands from larger to smaller, until the direct contribution of area to species richness (partial standardised regression coefficient b AS ) was zero. In agreement with Triantis et al. (2006), we think that drastically testing for a departure from zero may be not optimal and could lead to inappropriate conclusions in tests that lose statistical power when successive elimination of data is performed. For this reason, we considered the SIE limit to be reached when the estimated partial standardised regression coefficient was zero, regardless of statistical significance.
The Mantel (1967) test for matrix comparison was used to study the relationship between faunal similarity (as a proxy of the taxonomic structure) and environmental heterogeneity. In this case, we followed two approaches. In both approaches, the Jaccard (1908) index was applied to the faunal data. However, in the first approach, the Jaccard index was also used to estimate the similarity between islands in terms of presence/absence of biotope types on the islands. In the second approach, biotope-type-abundance data were used to calculate the Bray-Curtis (Bray & Curtis, 1957) which takes into account the quantitative aspect of data. Despite that the Mantel test is appropriate when testing hypotheses based on distances, it has received criticism mainly because it shows lower power than other analyses as linear correlation, regression, and canonical analysis (Legendre & Fortin, 2010). This implies that the Mantel test has a lower ability to detect a relationship when one is present in the data. We used it here for comparative purposes, as the test has been previously used in ecological studies (Malhotra & Thorpe, 1997;Sfenthourakis, 1996b). Additionally, following a different approach, we used the canonical correspondence analysis (CCA, ter Braak, 1986) to provide insights into the role played by environmental heterogeneity in shaping the ecological structure of the terrestrial isopod taxocenosis of the Pontine archipelago. Such a method may prove useful to study the role of ecological and geographical factors in explaining arthropod species composition on islands (Fattorini, 2011). Ecological categories of species per island as in Gentile et al. (2019) were used as response variables (see caption of Fig. 3 for an explanation). Biotope-type abundances were the environmental descriptors used as explanatory variables (Table S3). Both response and explanatory variables were used as percentages. Given the opportunistic sampling design underlying the present study, we did not use CCA in a rigorous hypothesisdriven context that implies rejecting the null hypothesis that ecological categories of species are unrelated to environmental factors (biotope types). We rather aimed to explain the variation and explore patterns to highlight possible relationships between ecological categories of species and environmental factors. Even if CCA is used for exploratory purposes, it requires some attention to avoid the introduction of poorly informative or redundant variables. For this reason, in this approach, halophilous species were excluded as their records could be incomplete because access to rocky shores and cliffs may be difficult or impossible in many cases. Consequently, associated environments (biotope 7, beaches with halopsammophilous vegetation; biotope 13: unvegetated areas of cliffs) were also removed. Myrmecophilous species were also excluded due to likely incomplete sampling. Before performing CCA, we removed biotopes present only in two islands or less and occurring at a percentage lower than 1% (biotopes 3, 4, 11, and 12). We also investigated the correlation structure of explanatory variables, to identify and remove superfluous variables. Biotope 10 (Built-up areas and human activities) was removed as it was strongly correlated with biotopes 8 (r = 0.773; p < < 0.01) and 9 (r = 0.857; p < < 0.001). We chose to use Type 2 scaling to graphically report CCA results. Type 2 scaling emphasises the relationships among response variables (ecological categories of species) and between response and explanatory variables (these last biotope-type abundances). Type 1 scaling would emphasise the relationships among objects (islands), which is not the purpose of this work.
We investigated the correlation between geographic distance and species richness in the two different sets of islands: large (5) and small (7). We used only distances measured as the minimum distance from the mainland (in case of large islands) or from the closest larger island (in the case of islets). We then followed a correlative and regression approach, using log-transformed data.
Statistical analyses were carried out using PAST ver 3.23 (Hammer et al., 2001), freely available at the website

Results
Whereas we found no correlation between geographic distance and species richness in the first set (large islands), we found a strong and significant negative correlation (r = − 0.826; p < 0.01) in the second set (islets). In this case, data fitted well a linear regression model (LogS = − 3.78 LogDistance + 0.428; R 2 = 0.683; p = 0.01). The faunal similarity was not statistically significantly correlated to environmental heterogeneity when Jaccard index was applied, whereas it was when biotope abundance was used (r Bray-Curtis = 0.431; p < < 0.01).
Partial correlation coefficients and standardised regression coefficients of the standard multiple regression are summarised in Table 1. In the regression, area did not contribute indicates the maximum value (abundance) of the response variable along that explanatory variable. The vector length is multiplied by a factor of 3, for better reading. Endogeous: EDG; euriecious: EUR; humicolous: HUM; littoral: LIT; synanthropic: SYN; xerophilic: XER Table 1 Standard multiple regression. The multiple correlation coefficient R is in the first column. Partial correlation coefficients for A and H are in the second and third columns, respectively. Standardised regression coefficients for A and H are in the fourth and fifth columns, respectively. Intercept is in the sixth column. Tolerance resulting from including both A and H in the regression is in the last column. As the regression includes only two variables, tolerance is identical for both of them. *** indicates p < < 0.001; ** indicates p < 0.01; * indicates p < 0.05. Coefficients with marginal statistical probability values (0.07 > p > 0.05; one-tailed test) are reported in italic Results of path analyses run on the complete island datasets of the Pontine and Sicilian archipelagos are shown in Table 2 (A and B, respectively). For the Pontine Islands, in no case, the partial standardised regression coefficient b AS (not statistically significant) was greater than b AH or b HS (Table 3). VIFs indicated not negligible collinearity in some cases, yet always under the established threshold. For the Pontine Islands, the upper limit of the SIE (the T value beyond which b AS become negative) was found at 1.3 and 1 km 2 , depending on the indicator used for environmental heterogeneity. Interestingly, no positive effect of A on H or of H on S was detected when the indicator of environmental heterogeneity was 1-D. We could not find a SIE for the Sicilian archipelago. In this case, b AS resulted negative when all islands were included (see Table 2, B). Figure 3 shows the pattern suggested by CCA. Most of the variance was explained by the first two axes (63.65% by axis 1 and 19.33% by axis 2). Complete results of CCA are reported in Table S4.

Discussion
In island biogeography, the occurrence of a positive relationship between species richness (S) and area (A) is a paradigm. Similarly, the occurrence of an inverse relationship between among-island geographic distance and species richness could be reasonably expected. The influence of geographic distance is in general difficult to disentangle. For terrestrial isopods, geographic distance seems to play little or no role in species richness of Oniscidean faunas of circum-Sardinian (Balletto, 1996), Tuscanian (Gentile et al., unpublished data), Sicilian archipelagos (Cazzolla Gatti et al., 2018), and Aegean islands (Sfenthourakis, 1996b). No evidence of a correlation between β-diversity and geographic distance was found for the Sicilian archipelagos (Cazzolla Gatti et al., 2018), whereas an inverse relationship between geographic distance and faunal similarity has been observed (Sfenthourakis, 1996b). In our case, no correlation between distance and species richness was found for the larger islands. Besides the obvious consideration that the limited number of large islands would not provide enough power to detect minor correlations, a role could be played by several factors such as the geographic proximity to the mainland, the short distances between islands, and the possible passive transport mediated by ferries connecting islands and the mainland. This last could contribute to an accelerated dispersal for synanthropic species. Indeed, a combination of anthropogenic influences (ballast water; tourism; soil, fertilisers; and ornamental plant transportation) and natural events (rafting on wood or floating seeds, island hopping) has been invoked by Cazzolla Gatti et al. (2018) as factors homogenizing species diversity. Indeed, for very small islets, which host a very low number of species and are all very close to a larger island with higher species richness, we found a strong negative correlation well modelled by linear regression of log-transformed data. Perhaps, under these circumstances, the effect of geographic distance may become evident, without the bias introduced by other confounding factors.
Apparently, these results could be consistent with the random placement model by which stochastic processes of colonisation would lead to the absence of a recognisable relationship between geographic distance and species number. Furthermore, these processes would also cause the species composition of a certain island to be independent of the geographic distance from other lands from which new colonisation can start. This model implies that the rate of species turnover would be erratic and unpredictable, and no equilibrium between extinction-colonisation (as emphasized by MacArthur & Wilson, 1967) is necessarily reached. The applicability of this model to Oniscidea of the archipelagos here mentioned could be a reflection of the rather low rate and passive nature of dispersal characterising many of the species belonging to this group. Species number is predicted to be much more affected by geographic distance in species with high and active dispersal (Taylor, 1987). In agreement with this view, species richness was not found related to geographic distance in tenebrionid beetles (Fattorini, 2002), Amphibians, and Reptiles (Parlanti et al., 1988) from several Mediterranean archipelagos. However, the random placement model, which can generate a positive correlation between species number and area, also denies the importance of environmental heterogeneity as a determinant of species richness (Connor & McCoy, 1979). This is in contrast with our results that showed a strong correlation between environmental heterogeneity and both area and species richness.
Indeed, in a model that included both environmental heterogeneity and area linear combination to better predict species richness, area did not contribute significantly to the model, although probability values were marginal. Anyway, partial correlations and standardised regression coefficients would still indicate that environmental heterogeneity would mostly contribute to species richness, except when the descriptor of environmental heterogeneity is 1-D that is more sensitive to the evenness of the island, whereas the Shannon index is also sensitive to the number of biotope types. Thus, not surprisingly, the contribution of H to S is predominant when H was described by B, LogB, and Shannon index.
Interestingly, for both Pontine and Sicilian archipelagos, a direct effect of A on S (b AS > 0) was not statistically supported by the path analysis. In fact, at the geographic scale considered in this study, a direct effect of area on species richness, if it exists, would be much less important than both the direct effect of area on environmental heterogeneity and the indirect effect on species richness. In general, our outcomes are in agreement with Cazzolla Gatti et al. (2018) who found the same results for the Sicilian archipelagos.
When the path analysis was performed by successive elimination of data, the SIE was detected for the Oniscidean fauna of the Pontine archipelago and the threshold (T) was found at A ~ 1.3-1 km 2 , depending on the indicator used to estimate environmental heterogeneity. Such value is of the same order of magnitude and only slightly lower than values observed by Triantis et al. (2006) and Sfenthourakis and Triantis (2009), between 3 and 5 km 2 , for terrestrial isopods from islands of the Aegean Sea. It is possible that the difference in T values between the Pontine and Aegean archipelagos, if any, may reflect different topography, palaeogeographic histories and proximity from the mainland of the two archipelagos. For the Pontine Islands, T, as estimated by path analysis, is also slightly different from the estimate provided by Gentile and Argano (2005). However, it has to be reminded that path analysis allowed us to calculate and extract the net effect of area. Such analysis indicated that species richness of islands beyond the threshold is sensitive to environmental heterogeneity, as mediated also through indirect area effects. This is conceptually different from the SIE threshold as found by Gentile and Argano (2005), based only on the species/area relation. A different indication is obtained for the Sicilian islands for which a SIE was not found, contrary to outcomes of the study by Gentile and Argano (2005), who found T ranging between − 0.3 and 0.5, which would correspond to 0.5 and 3.2 km 2 , respectively. However, their breakpoints were found in an area range where there were no islands. Cazzolla Gatti et al. (2018) questioned the existence of SIE for the Sicilian islands, and our results would support their conclusion. Perhaps a more refined definition of environmental heterogeneity which accounts also for the percentage of abundance of different biotopes could shed more light in this regard.
Mantel test and CCA provide an insight into the role of environmental heterogeneity in shaping the taxocenosis and ecological structure of Oniscidea from the Pontine Islands. The similarity matrices based on the faunal and environmental data of the islands were significantly and positively correlated but only when biotope-type abundance was taken into account. This is interesting because Cazzolla Gatti et al. (2018), although they did not use Mantel statistics, did find a positive correlation between faunal and environmental similarities, both being estimated by the Jaccard index, for Oniscidea from Sicilian islands. Indeed, we confirmed such a relationship after running the Mantel test on their data (r Jaccard = 0.388; p < < 0.01). Thus, although it has to be used with caution, the Mantel test proved conservative in this case. A clear and expected pattern emerged from the CCA, illustrating that littoral and euriecious species are mostly related to arid and poorly vegetated environments, which are most abundant in small islets. On the contrary, endogeous and humicolous species are less abundant in such environments and are more linked to environments with richer soil horizons, most abundant in larger islands. Xerophilic species are related to arid or semi-arid herbaceous environments and abandoned cultivations and, to some extent, to Mediterranean maquis that, in the Pontine archipelago and in particular in Ponza, was subject to human impact represented by cutting and fire. Perhaps this last consideration may explain the relationship of synanthropic species with this biotope type.
We are aware of the limitation of the use of biotope type data as a proxy of environmental heterogeneity relevant to the species investigated (see Dengler, 2010 andSfenthourakis, 2012). Actually, there is evidence indicating that terrestrial isopods are not evenly distributed in biotope types similar to those here considered. Isopods may vary in the number of species as well as of specimens (Achouri et al., 2008). When trying to describe environmental heterogeneity, it should be considered that, in human-inhabited islands, biotope types may be reduced and disappear more rapidly than the invertebrates species do, causing a potential interpretation bias. Such a bias would be reduced when just the number of biotope types is used. However, we can also speculate that abundance data, in the definition of the contribution of different types of biotopes to environment heterogeneity might provide additional information otherwise lost in case only presence/ absence environmental data are considered and the size of islands is really small. In fact, a biotope type in small islands, yet present, could be not sufficiently large to support one or more species. Considering biotope abundance in terms of percentage could be a reasonable compromise, at least for organisms with limited dispersal ability. Scale effects and scale dependency of the species/area relationship have been recently described, and "habitat quality" in islands was also found important at local scales (Schrader et al., 2019). The large difference in size between the Sicilian and Pontine archipelagos (Sicilian mean rank = 12.3 km 2 ; Pontine mean rank = 3.2 km 2 ; U Mann-Whitney = 5; z = 4.394; p < < < 0.01) would be in agreement with such an interpretation. Such prediction remains to be investigated.
In conclusion, it is conceivable that the colonisation of large Pontine Islands may occur via processes independent of geographic distance which could instead be an important factor at a much smaller scale. Such processes may be driven by a combination of anthropogenic influences and natural events. Even in very small-size island systems, environmental heterogeneity mostly contributes to species richness. Environmental heterogeneity could influence the taxocenosis structure and, ultimately, the number of species of Oniscidea via direct and indirect effects, these last mediated by area which may or may not have a direct effect on species richness. Our data are in agreement with other studies regarding other Mediterranean islands and, on the whole, they might be indicative of a general pattern, although the same approach remains to be applied to other Mediterranean archipelagos.