A sampling optimization analysis of soil‐bugs diversity (Crustacea, Isopoda, Oniscidea)

Abstract Biological diversity analysis is among the most informative approaches to describe communities and regional species compositions. Soil ecosystems include large numbers of invertebrates, among which soil bugs (Crustacea, Isopoda, Oniscidea) play significant ecological roles. The aim of this study was to provide advices to optimize the sampling effort, to efficiently monitor the diversity of this taxon, to analyze its seasonal patterns of species composition, and ultimately to understand better the coexistence of so many species over a relatively small area. Terrestrial isopods were collected at the Natural Reserve “Saline di Trapani e Paceco” (Italy), using pitfall traps monthly monitored over 2 years. We analyzed parameters of α‐ and β‐diversity and calculated a number of indexes and measures to disentangle diversity patterns. We also used various approaches to analyze changes in biodiversity over time, such as distributions of species abundances and accumulation and rarefaction curves. As concerns species richness and total abundance of individuals, spring resulted the best season to monitor Isopoda, to reduce sampling efforts, and to save resources without losing information, while in both years abundances were maximum between summer and autumn. This suggests that evaluations of β‐diversity are maximized if samples are first collected during the spring and then between summer and autumn. Sampling during these coupled seasons allows to collect a number of species close to the γ‐diversity (24 species) of the area. Finally, our results show that seasonal shifts in community composition (i.e., dynamic fluctuations in species abundances during the four seasons) may minimize competitive interactions, contribute to stabilize total abundances, and allow the coexistence of phylogenetically close species within the ecosystem.


Introduction
Biological diversity analysis is frequently used to describe both communities and regional species compositions (Magurran 1988(Magurran , 2013 and is at the basis of many ecological models (MacArthur and Wilson 1967;Connell 1978;Stevens 1989). Evaluating biodiversity is not only important for drawing comparisons among sites (Cornell 1999), but also to estimate population sizes and to understand community compositions (Sanders and Entling 2010). Moreover, knowing which species occur in a given region is sometimes fundamental for conservation purposes and requires an accurate understanding of their distributions (Cardoso 2009). Maximizing species richness within a reserve is often a goal of conservation efforts (May 1988). Arthropods are a megadiverse group, and their abundance and diversity make it almost impossible to fully assess their richness, their functions in the ecosystem, and their geographical patterns (Ramos et al. 2001). The Oniscidea, or terrestrial isopods, are abundant and widespread components of the soil's fauna and play significant roles in soil ecology. They contribute in the regulation of organic matter and nutrients (Hassall and Sutton 1978;Sutton 1980;Zimmer et al. 2003) and are important elements of soil food webs, being themselves food sources for other arthropods (Vetter and Isbister 2006) as well as for vertebrates (Ben Hassine and Nouira 2009;Covaciu-Marcov et al. 2012). Moreover, due to their biological and ecological characteristics, terrestrial isopods are used as biological indicators of heavy metal pollution (Paoletti and Hassall 1999) as well as of grassland habitats quality (Souty-Grosset et al. 2005). Methods of data collection, including the sampling period, are highly important for studies focusing on terrestrial isopods diversity. Usually, collection of samples takes place over long periods of time, in order to assure an appropriate sample size for further analyses. As a consequence, the sampling procedure is generally time consuming and expensive in terms of materials and personnel (Sutherland 2006), while collecting very large numbers of specimens may reveal detrimental for conservation purposes. Even if a wide range of techniques is available for sampling isopods, these techniques need to be standardized to varying degrees to enable catches to be comparable over time or between areas (Sutherland 2006). This study tries to assess the optimum sampling effort and the seasonal timing necessary to maximize detected species diversity during data collection, while investing in fewer resources in terms of time, money, personnel, etc. Moreover, our research analyzes the spatiotemporal shifts in community composition in order to understand mechanisms allowing the coexistence of such a high number of Isopoda species in a relatively small area.

Sampling area
This study was carried out in the Natural Reserve "Saline di Trapani e Paceco" (Fig. S1), which includes a Special Protection Area ("SPA", for the purposes of the European "Birds Directive" 79/409/CEE), an "Important Bird Area", a RAMSAR-Convention protected wetland, as well as a EU Site of Community Importance ("SCI", ITA 01007). The reserve is located in the south of the Province of Trapani, on the west coast of Sicily (Italy); it covers a surface of 960 ha and consists of a flatland with sandy coasts and large wetlands.

Sampling method
Data were obtained following a standardized sampling method. Pitfall traps were randomly distributed all over the area and consisted of 29 sampling units, grouped in five spatially disjoined replicates (Koivula et al. 1999;Lassau et al. 2005). Traps were monitored monthly all through a 24-month period, from February 2008 to February 2010. The traps (diameter 10 cm, height 14 cm) were sunk into the ground, with their rims level with the soil surface, and were half-filled with a saturated water/sodium chloride solution, in order to avoid the attractive effects of formalin and vinegar. Pitfall trapping is a widely used method (Pek ar 2002), as it is time efficient, easy to use, and inexpensive, while it produces large species-rich samples suitable for statistical analyses (Spence and Niemel€ a 1994). It is broadly recognized as a valid sampling technique for the soil invertebrate fauna in general (New 1999;Brandmayr et al. 2005) as well as for isopods in particular (Becker 1975;Fleugge and Levens 1977;Al-Dabbagh and Block 1981;Caruso and Zetto Brandmayr 1983).

Statistics
Data of species richness and of individual species abundances obtained by each sampling unit were pooled separately for the 2 years and the four seasons of each year (see Table 1). We analyzed data of aand bdiversity by the software EstimateS (Colwell 2013). More in detail, for each species, we recorded the number of sampled individuals and evaluated the species richness of the community by Chao 1 index. For each assemblage, we calculated the relative dominance and McNaughton's dominance (which accounts for three most dominant species), as well as Shannon's and Simpson's inverse diversity and evenness. Classical Jaccard and Sorensen similarity measures were computed to better understand the complementarity of each coupled seasonal assemblage independently of each species' abundance (i.e., only on incidence data) and to disentangle the effects of dominance on the communities. We also calculated Whittaker's bW measure, which is widely recognized as the less affected by error among indexes of a-diversity (Cazzolla Gatti 2014). We computed abundance-based b-diversity measures (such as Bray-Curtis, Chao-Jaccard, and Chao-Sorensen), to understand the influence on similarity of the distribution of individuals among species. Finally, we calculated the Morisita-Horn index, which is known as the best performing abundance-based index of complementarity (Magurran 2013). For distributions of species abundances and accumulation and rarefaction curves, we followed the approach proposed by Magurran (2004) or by Cazzolla Gatti (2014) to analyze changes in biodiversity over time.
We checked the statistical significance of differences in the distribution of abundances between the rank-abundance plots of the 2 years by Kolmogorov-Smirnov twosample test.

Accumulation and rarefaction indexes
SACs (Species-area curves) are widely used in biodiversity research and can provide useful information to optimize sampling efforts and reduce resource wasting in future research. Although species accumulation curves, such as the abovementioned SACs, can be used to draw inferences about the diversity of a more fully censused assemblage (Gotelli and Colwell 2001), rarefaction curves allow to estimate richness at the abundance level of the smallest sample. For these reasons, we used rarefaction curves to better compare seasonal assemblages and speculate about their "real" richness, irrespectively of the sampled area. We used Coleman's type rarefaction curve, which estimates the number of species in samples, on the assumption that all individuals in all samples were randomly mixed (Chazdon et al. 1998). We used SACs also to compare the four seasons of each year.

ECDF (Empirical Cumulative Distribution Function)
As, when comparing samples with high differences in richness, most plot types tend to overemphasize differences in richness, and curves of the richest sites become stretched in the low right corner of the graph (i.e., plotting data on k-dominance graph), which makes them apparently more even, we represented abundance data with ECDF plots, which can better discriminate between assemblages by rescaling their ranks according to richness (Magurran 2013;Cazzolla Gatti 2014).

Species viability (population dynamics)
Species viability (i.e., population dynamics according to abundances) of the eight seasons of both years was used to describe patterns of species compositions during time.

SAD (Species abundances distributions)
Species abundances distributions are generally used to shed light on processes that determine the biological diversity of a species assemblage. Frequency distribution plots, where the number of species is displayed in relation to the number of individuals per species, permit to understand and compare dominance/evenness patterns. In these graphs, the mode usually falls on the lowest abundance class.

Rank-abundance plot
Rank-abundance plots (also known as Whittaker's plots) are among of the most informative methods to evidence contrasting patterns of species richness and to highlight differences in evenness among assemblages.

Results
A total of 25,690 individuals belonging to 24 species were collected. Table 1 shows a list of all the species sampled in each season.
The aand b-diversity indexes are, respectively, summarized in Tables 2 and 3, where the basic indicators (such as the number of species and the number of individuals) are shown together with dominance (relative dominance and McNaughton dominance), diversity and evenness indexes (Shannon's and Simpson's inverse diversity and evenness), and Chao 1. The SACs for each season of the 2 years are shown in Figure 1 and were calculated by plotting sampling units' increments against the number of species sampled, with 100 randomizations of sample units to obtain a smoothed curve. Accumulation and rarefaction curves are shown in Figure 2. The ECDF independence to richness is shown in Figure 3 for both years. For a clearer and more detailed analysis of population dynamics according to abundances, we separately plotted all the species (Fig. 4A), all but the most abundant one (i.e., Armadillidium granulatum, Fig. 4B), only the rare species (those with an annual total species abundance of n ≤ 100, Fig. 4C), and only the dominant species (having annual total abundance of n > 100) but excluding A. granulatum (Fig. 4D).
Frequency distribution plots for all, rare (n ≤ 200 individuals/year collected) and dominant (n > 200) species are presented in Figure 5A-C, respectively. These plots show the highest frequency in the rare species class, as well as differences between the 2 years.
As in the case of SADs, we plotted rank-abundance data for all, rare, and dominant species to better understand differences in species distributions occurring in each annual assemblage (Fig. 6A-C, respectively). On the base of these distributions, we also calculated Kolmogorov-Smirnov statistical tests.

Discussion and Conclusion
Our research focuses on the study of the Oniscidea, a taxon already widely studied in the Mediterranean region, most particularly in Sicily (Caruso 1968(Caruso , 1973a(Caruso ,b, 1974(Caruso , 1976Caruso and Lombardo 1976;Caruso et al. 1987;Messina et al. 2011), north Africa (Hamaїed-Melki et al. 2010Khemaissia et al. 2012a,b) and Greece (Alexiou and Sfentourakis 2013). The present work analyzes both the isopod diversity of the study area and the optimal approach capable to maximize our understanding of the contribution of species diversity while allowing to spend fewer resources for data sampling.

a-Diversity
As is shown in Table 2, in both sampling years, the highest species richness (S) occurs in spring (17 and 21 species, respectively). The highest total abundance (as value of N) is also shown to occur in the spring of both years (10,157 and 3406 specimens, respectively). In general, spring is the best growth season for most of the 24 species observed. According to the Chao 1 estimator (Table 2), in both years, the number of species expected    Table 2). Cases of population explosions (Warburg 1993) of Armadillidium vulgare (Latreille, 1804) in North America (Hatch 1947), A. granulatum on Panarea (Caruso 1968), and Armadillidium decorum at Collesano (Italy) are well documented in the literature. In 2009, on the contrary, no particularly high dominance of one or a few species is observed, and the maximum diversity level is reached between summer (Simpson inv.) and autumn (Shannon). This latter result is due to the well-known effects of highest abundances on Simpson index (which, in summer 2009, is influenced by a triple total number of individuals compared to that of the autumn 2009). Summarizing, in terms of richness, spring is the best season for the analyzed species, while, considering evenness values, diversity was maximum between summer and autumn (when dominance of the most abundant species decreases).

b-Diversity
Analyzing changes in b-diversity values along time, by pairwise summing in turns two of the four seasons of each year (Table 3), in parallel with incidence indexes (presence/absence data only and ignoring abundances), for 2008 the best combination of the coupled seasons is spring-autumn (where both Jaccard and Sorensen indexes are lowest, which means less similarity and more diversity, Table 3). In 2009, in contrast, the best combination is slightly anticipated by one season (spring-summer), although b-diversity remains high also in spring-autumn. Anyway, all these measures are strongly influenced by species richness and by sample sizes. To increase the reliability of our results, we also calculated Whittaker's (Bw) index (a b-diversity incidence measure, Table 3), which provides values that are less influenced by errors and has fewer restrictions. This index confirms the best combinations of seasons shown by Jaccard and Sorensen indexes. Because our sampling cannot be considered totally complete (SACs were close to saturation but, apart from the case of summer, did not reach plateau Fig. 1) and differences in the number of individuals of the dominant species play a major role in abundance distributions ( Fig. 5A-C), the abundance-based measures of b-diversity (Chao-Sorensen, Chao-Jaccard, and Bray-Curtis) gave contradictory results, so that they cannot be considered reliable indexes of diversity, in this case. The same considerations apply to the Morisita-Horn index, which is strongly influenced by the abundance of the dominant species in the sample. In an attempt to provide the best possible representation of c-diversity, and eventually optimize time and resources, if we compare the highest number of species (S) collected in the coupled seasons with results shown by the indexes (Table 3), it appears that only the three incidence measures really reflect the highest values of S in each couple (20 in spring-autumn 2008 and 22 in spring-summer 2009).In summary, b-diversity is maximized if samples are collected first during the spring and then between summer and autumn. The

Accumulation and rarefaction
The analysis of species areas curves (Fig. 1) shows that our sampling (five replicates of 29 sampling units) is very close to saturation in every season, and reaches complete plateau in the summer of both years. Also in both years, spring is the season that shows the highest values of species richness (as shown also by aand b-diversity analyses) irrespectively of the number of collected samples, whereas the winter of 2009/2010 shows the lowest increments. Indeed, three of the four species not collected during the latter season (Tylos ponticus, Halophiloscia hirsuta, and Stenophiloscia glarearum) live near the shoreline and their occurrence depends on the weather and the conditions of the sea (Vandel 1962), which were particularly harsh in the winter 2009/2010. The analysis of Coleman's rarefaction curves (Fig. 2) shows that if, to optimize the sampling effort, we consider the minimum number of individuals as our reference value, we get controversial results. This may be due to the large fluctuations of abun- dances during the year and to the frequency of oscillations in different years (see Figs. 5,6). For instance, the low number of individuals recorded in the autumn 2008 was due to a strong decrease in the populations of the dominant species, which was not observed in 2009. This is the reason why, if we rarefy the number of individuals to a minimum, autumn shows the highest richness in 2008 but not in 2009 (when dominant species populations are more evensee populations dynamics).In summary, while accumulation curves (SAC) show that the sampling dimension utilized is adequate for a reliable representation of c-diversity, rarefaction in this case is inadequate to suggest an optimized sampling effort, being strongly influenced by abundances and, in our case, by those of only one or a few dominant species. Anyway, rarefaction offers valuable information on the minimum overall number of individuals needed to sample the highest richness in each season. What emerges, in fact, is that while less than 2000 individuals are enough to represent the richness of the community during summer, autumn, and winter, this number has to be increased in the case of spring.

ECDF
The ECDF affords a better understanding of the distribution of species abundances, it is less influenced by richness, and it is also more reliable than Whittaker plots. ECDF shows that in 2008 (Fig. 3A), spring accounts for the highest number of rare species (bottom left part of the curve), but at the same time also for the lowest (central part of the curve with a gradual slope, with less verticality) and more evenness, as shown by a-diversity analysis. In 2009 (Fig. 4B), instead, winter evenness is the lowest because of the absence of some species, while total abundance remains high, as also shown by aand b-diversity analyses, but the number of rare species is lower than in spring and summer.

Species viability (population dynamics)
Analyzing species population dynamics (Fig. 4), a general trend of population decrease emerges in the autumn of both years and for all species, although, apart from Porcellio laevis (in the dominant group), the rare species seems to be less prone to a reduction, if compared with dominant species. These fluctuations can explain the observed patterns of aand b-diversity: The dominant species generally decline in autumn, whereas rare species are stable in this season, and even if their populations fluctuate, often disharmoniously so, during the rest of the year. Species that dominate the growing season of the first year seem to reduce their dominance in the next year, and vice versa. At the same time, the abundance of the dominant species in the first year (A. granulatum) is higher than that in the second year, whereas, apart from A. granulatum and as also shown by aand b-diversity, the abundances of both rare and dominant species are higher in the second year, but with a lower number of species. In agreement with Shimadzu et al. (2013), our results show that spatiotemporal shifts in community composition can minimize competitive interactions, increase some biodiversity values, and help stabilize total abundances.

SAD
To gain better understanding of the patterns shown by aand b-diversity and by population dynamics, we analyzed SADs also as species abundance histograms and Whittaker plots, by splitting the data into three categories: total abundances, dominant abundances, and rare abundances. The simple plotting of species and abundances (Fig. 5) shows a low number of dominant species (8) compared to the high number of rare species (16) in both years (Fig. 5A). In 2008, the three most abundant species show even higher dominance than in 2009 (Fig. 5C). The other five dominant species may suffer the dominance of the three most abundant ones, as they show reduced populations in the first year, but not in the second, when the three most dominant species accounted for a lower number of individuals (Fig. 5C). Figure 5B shows that, in contrast, the abundances of rare species (n ≤ 200) are higher when populations of the three most dominant species are lower (second year). Moreover, the number of rare species increases when dominance is reduced (second year, Fig. 5B). Dividing species and their relative abundances into three groups (total, dominant, and rare), Whittaker plots show that while both communities maintain the same general distribution in both years (Fig. 6A), the dominant species (n > 200) plot (Fig. 6B) confirms a less even distribution and an higher dominance (steeper slope) in the first year (2008) as compared to the second one (2009). The rare species (n ≤ 200) plot (Fig. 6C), instead, shows no visual difference in distribution between years. Anyway, the cumulative two-year distributions show no significant statistical difference (D21,23 = 0.201 P > 0.01) when checked with Kolmogorov-Smirnov test and this confirms that the distribution patterns of the two communities in 2008 and 2009 were similar and that our results can be considered valid for both years. This closer look at the distributions of dominant and rare species allowed us to get a detailed picture of the assemblages of the analyzed communities. It also confirms observation derived from aand b-diversity measures, which we made with the aim to provide advice on how to reduce the sampling effort, as well as to save time and economic resources, without losing information when analyzing soil-bugs biological diversity.