Comparison of terrestrial isopod (Crustacea: Oniscidea) assemblages from two preserved areas (Bouhedma and Chambi) in arid regions

Abstract A total of 1290 specimens belonging to 11 species and four families were captured in two national parks, Bouhedma and Chambi. Among these species, seven belonged to Porcellionidae, and the remaining species represented Agnaridae, Armadillidae and Armadillidiidae families. Five species were common and recorded in both parks. However, Armadillidium tunisiense, Hemilepistus reaumurii and Porcellio djahizi were recorded only in Chambi while Agabiformius lentus, Armadillo officinalis and Porcellio albinus were collected only in Bouhedma. The distribution structure of the collected species was analyzed according to altitude and plant assemblages. Seasonal sampling showed that the highest abundance and species richness were recorded in spring. In both parks, the species richness decreased as the altitude increased. Arid regions sheltered specific species such as H. reaumurii and P. albinus, which were often the dominant component of the arthropod macrodecomposer guild in some habitats. The similarity analysis showed a quantitative and qualitative difference between the two parks. The two parks Bouhedma and Chambi shared five species (Leptotrichus panzerii, Porcellio laevis, P. variabilis, Porcellionides pruinosus, Armadillidium sulcatum) with areas studied in the north of Tunisia, Kroumirie, supralittoral zones and around the wetlands.


Introduction
The diversity of terrestrial isopod species is directly connected with the number of microhabitats, which allow the coexistence of several species (Lopes et al. 2005). Terrestrial isopods, as members of the soil fauna, are expected to respond to environmental changes given their limited dispersal abilities. They are probably very sensitive to changes in soil chemical and physical properties (Paoletti & Hassall 1999), thus responding rapidly to environmental changes (Nakamura et al. 2003). Many studies have shown the influence of the general characteristics of the area on the terrestrial isopods Mureşan et al. 2003;Tomescu et al. 2008Tomescu et al. , 2011. Hence, knowledge of the isopod faunal richness of a region is one of the factors responsible for success in monitoring and conserving biodiversity (Lopes et al. 2005).
In Tunisia, some studies were performed on reproductive behavior, systematics and geographical distribution of species of the genera Hemilepistus Porcellio, Porcellionides and Armadillidium (Rezig & Nasri 1992;Medini & Charfi-Cheikhrouha 1998;Medini-Bouaziz et al. 2006. In the last few decades, other investigations have been performed as part of a broader study aiming to assess the biodiversity of terrestrial isopods in humid zones in Tunisia (Fraj et al. 2010;Hamaïed-Melki et al. 2011;Hamaied et al. 2013).
Protected areas have restrictions on human activities aiming at preserving biotic and abiotic components of the landscape. The present study contributes to our knowledge of the diversity of terrestrial isopod communities in two national parks, Bouhedma and Chambi, in arid regions in Tunisia. Both Kremen et al. (1993) and Massa and Ingegnoli (1999) reported that "the diversity and abundance of terrestrial arthropods can provide a rich base of information to aid efforts in the conservation of biodiversity and the planning and management of nature reserves". In addition, the isopods in the remaining forests and steppes facing desertification, if well preserved, may play an important role in the local food webs.
Thus, we intend to present data on the community of terrestrial isopods in arid regions, which are currently insufficiently studied, and threatened by drought and desertification due to global climate change. Specific objectives were (1) to make a list of terrestrial isopod species recorded in both Bouhedma and Chambi parks; (2) to analyze their distribution according to season, altitudinal gradient and plant associations; and (3) to determine which species are present in these areas, and also to describe any clear eco-morphological characteristics (Schmalfuss 1984) associated with different habitat types.

Sampling areas
The first site, the Bouhedma National Park (Figure 1), is located in southern Tunisia (34°39°N, 9°48°E). It stretches 85 km east of Gafsa, 100 km south of Sidi Bouzid, 105 km west of Sfax, 100 km northwest of Gabes. It covers an area of 16,488 ha of which 8804 ha are fully protected. It is considered among the last savanna gum in North Africa covered by Acacia tortilis subsp. raddiana with Retama raetam, Cenchrus ciliaris and Digitaria commutate (Zaafouri et al. 1994). In 2010 it became a part of the "World Network of Biosphere Reserves" (WNBR) after having already been a part of "Man and the Biosphere Programme" (MAB) since 1971. Bouhedma Park is located in an arid bioclimatic zone (Chaieb & Boukhris 1998) where the mean annual rainfall varies between 100 and 200 mm (Abdallah et al. 2008). In this area, three habitats were selected according to an elevation gradient ranging from 120 to 210 m and also according to vegetal associations ( Table I).
The second study area, the National Park of Chambi (Figure 1), which has been a part of the WNBR since 1977 and a part of the MAB since 1971, belongs to the governorate of Kasserine. Its total area is 6723 ha. This park has the distinction of being representative of all the conditions of the Tunisian Atlas. It belongs to the area of the central highlands and has the highest peak in Tunisia, namely "Jebel Chambi", with an altitude of 1544 m.
The vegetation of Chambi consists of a remarkable assemblage of plant groups ranging from the green oak (Quercus ilex) in the summit rocks to the steppe alfa (Stipa) in the piedmont. The undergrowth of the pine forests of Aleppo consists of Juniperus phoenicea, J. oxycedrus, Cistus, Globularia and Rosmarinus. The National Park of Chambi is divided into seven zones according to altitudinal gradient and plant associations (Table I). In this park, seven habitats were investigated, one in each zone.  (Table II). In each habitat, specimens were collected only once per season (Table II), using 30 (0.5 × 0.5 m) quadrates randomly distributed. The rocks, soil and dead leaves were lifted using a small pallet, while  (*) January = winter; April = spring; October = autumn; July = summer.
the tree barks and palm trees were removed manually. In order to compare the results of this study to those obtained previously (Achouri et al. 2008a;Hamaïed-Melki et al. 2011;Hamaied et al. 2013), and to get the same probability of capture, individuals were collected by three authors during a period of 1 hour. This time-effort procedure offered semi-quantitative information (Souty-Grosset et al. 2008). To complete the list of terrestrial isopods in the two parks, specimens were also collected by hand, but these were not included in the quantitative data.

Laboratory analysis
Collected specimens, preserved in 70% ethanol, were identified and sexed in the laboratory according to Vandel (1962) using a Leica M5 stereomicroscope (Leica Microsystems, Bensheim, Germany).

Data analysis
Different ecological indices estimating the diversity in the studied habitats were used, such as the Berger-Parker dominance index (May 1975), with d = Nmax⁄N, (Nmax, number of individuals of the most abundant species; and N, the total number of individuals). Total abundance of terrestrial isopods (N), number of species (S), Shannon-Wiener diversity index (H′) (Shannon & Weaver 1963), and Pielou's evenness index (J′) (Pielou 1966) were calculated for each sample. The difference between diversity indices (H′ and J′), across season and different altitude, was tested using one-way analysis of variance (ANOVA). Non-metric multidimensional scaling (MDS; Kruskal & Wish 1978) was performed based on mean abundance of species at each sample. Thus, MDS could be considered an alternative to factor analysis. In general, the goal of the analysis was to detect meaningful underlying dimensions that allowed us to explain observed similarities or dissimilarities between the investigated variables. In factor analysis, the similarities between variables were expressed in the correlation matrix. With MDS, we could analyze any kind of similarity or dissimilarity matrix in addition to correlation matrices. Furthermore, Similarity Percentages (SIMPER) broke down the contribution of each species (or other variable) to the observed similarity (or dissimilarity) between samples. It allowed us to identify the species that were most important in creating the observed pattern of similarity. The PRIMER 5.0 (Plymouth Routines in Multivariate Ecological Research) software package (Clarke & Warwick 1994) was used to construct a complete linkage cluster analysis. Canonical correspondence analysis (CCA) was developed to allow ecologists to correlate the abundance of species to environmental variables (Ter Braak Cajo 1986). In our study, CCA was applied to correlate the species abundance to the altitudinal gradient and vegetal associations. This analysis was performed using the package Excel XLSTAT v. 7.5.2. Table III. Spatial distribution of isopods in the two parks (Bouhedma and Chambi), relative abundance, and values of Shannon-Wiener (H′) and evenness indexes (J′).

Species composition
In the studied areas, 1290 individuals were collected: 347 from Bouhedma Park and 943 from Chambi Park. They belong to 11 species in the families Agnaridae, Porcellionidae, Armadillidiidae and Armadillidae (Table II) (Table III).
Species richness, abundance and diversity Spatial variation. The total numbers of each species recorded are summarized in Table III. Species richness ranged between 1 and 7. The habitats Bouh 1 and Cham 2 exhibited the highest species richness (7; Figure 2), while the habitat Cham 7 exhibited the lowest (1). Analyzing the spatial variation data by nested ANOVA, we found that there were no statistically significant differences between habitats in overall species richness of isopods recorded in the two parks (F = 1.627, df = 16.28, p = 0.22). In fact, the species richness was negatively correlated with altitude (r 2 = 0.6452), and no relationship was found between the Shannon-Wiener index (H′) and altitude (r 2 = 0.1543). However, significant differences were found between species abundance related to altitude gradient (one-way ANOVA, F = 7.348, df = 9.008, p = 0.02395). Temporal variation. In terms of abundance, the ANOVA revealed significant differences between the parks throughout the years (F = 24.57, df = 2, p = 0.03) and between seasons (F = 12.31, df = 3, p = 0.03). In both studied parks, H. reaumurii was the most abundant species (37.59%), followed by Porcellio albinus (16.43%). Armadillidium tunisiense and P. djahizi were less abundant than the two previous species (12.36% and 11.53%, respectively). Despite its high abundance, H. reaumurii was only found in the first habitat in Chambi Park which was characterized by a sandy loam soil with Atriplex and other halophytes. Armadillidium tunisiense and P. djahizi occupied almost the same habitats, Cham 2 to Cham 7, which were characterized by high elevation (900 to 1500 m). However, P. variabilis and L. panzerii were the most frequent, and they were collected in five and six habitats, respectively, in the two parks (Table III). The Shannon-Wiener (H') and evenness (J') indexes showed seasonal fluctuations between habitats, with significant differences (F = 17.78, df = 10.65, p = 0.00; and F = 24.55, df = 9.23, p = 0.00, respectively).

Assemblages and species affinities
The MDS ordination analysis from quantitative data (Figure 3a) showed the existence of three groups (I,  II  Quantitatively, a significant difference was found between seasons in the global result (global R = 0.76, p = 0.01). The pairwise test results showed no significant differences between autumn-spring (R = 0.083, p = 0.4) and autumn-winter (R = 0.18, p = 0.1), but significant differences were found between winter-summer (R = 0.29, p = 0.013), winter-spring (R = 0.37, p = 0.011) and autumn-summer (R = 0.35, p = 0.045). The SIMPER analysis indicated that in autumn and summer, the species contributing most to the difference between seasons were P. albinus, P. djahizi and Ar. tunisiense, while in winter and spring, the most contributing species was H. reaumurii. The highest dissimilarities of SIMPER analysis were obtained between winter-summer (85.48%), winter-spring (57%) and autumn-summer (56.05%).
The MDS ordination analysis in Figure 3b, based on the abundance of collected species, suggested that the isopod assemblages could be grouped in two categories. The first, Group A, contains the seven habitats of Chambi Park, which revealed discrimination among three subgroups. The first subgroup (Cham 1) was characterized by the lowest altitude with halophyte vegetation and sheltered only H. reaumurii. The second was made of Cham 2, 3, 4, 5 and 6 where the diversity was inversely proportional to the altitude. The third one (Cham 7), with the highest altitude, showed a plant association of Quercus ilex and Aleppo pine, where only A. tunisiense was collected. The second, Group B, is composed of the three habitats of Bouhedma Park. Bouh 1, a reforestation area of Eucalyptus and Acacia, exhibited the highest species richness (7), compared to Bouh 2 and 3, dominated by Gramineae, C. ciliaris and D. commutate. Nonetheless, the difference in abundance among these habitats was not significant.

Distribution of isopods according to altitude, plant association and soil type
The CCA was performed to correlate the distribution of isopod species to vegetation cover and altitude. Our analysis showed that 71.49% of the information (Figure 4) was obtained by adding the first two axes together, of which 53.38% was displayed by the first axis alone. Cham 1 was positively correlated to the first axis (F1) and habitats Cham 3 to Cham 7, characterized by high altitude and vegetal associations composed of Aleppo pine forest, Juniperus phoenicea, J. oxycedrus and Quercus ilex, were also positively correlated with F1. However, Bouh 1 was negatively correlated with F1. Bouh 2, a degraded meadow formed by Acacia tortilis, and Bouh 3, an open area dominated by two Gramineae, Cenchrus ciliaris and Digitaria commutata, were positively correlated with F2, whereas Cham 2, characterized by a steppe of Stipa, Globularia and Rosmarinus on limestone soil, was negatively correlated with the axis F2. Hence, H. reaumurii was clearly associated with Cham 1 characterized by a vegetal association of Atriplex and other halophytes on a sandy loam soil. Armadillidium tunisiense and P. djahizi were clearly related with Cham 3-7. Armadillo offcinalis and P. albinus were related to Bouh 2 and Bouh 3, while A. lentus, A. sulcatum, P. laevis, P. variabilis and P. pruinosus were related to Bouh 1 and Cham 2.

Discussion
The composition of isopod assemblages differed in the two parks. Despite their relatively equal species richness (8), Bouhedma and Chambi showed a quantitative difference in collected specimens -347 and 943, respectively. Within the studied parks Bouhedma and Chambi, the five common species were represented by a different number of specimens depending on habitat and sampling period. Leptotrichus panzerii and P. variabilis were the most common species, collected in five and six habitats from the two parks, respectively. The presence of P. variabilis in both parks suggests that this species is able to use a variety of habitats, as previously reported (Medini & Charfi-Cheikhrouha 1998;Medini-Bouaziz et al. 2015). Armadillidium sulcatum, collected in all habitats in Bouhedma and in the second habitat of Chambi, characterized by morphological and ecological variability in Tunisia (Hamaied-Melki 2008), was capable of adapting to habitats with harsh conditions. As they were usually collected in places linked to human activities, the presence of P. pruinous and P. laevis in Bouh 1 and Cham 2, more exposed to anthropogenic impact, can be explained by the adjacent motels and other buildings which are favorable habitats for synanthropic species.
Besides the species common to both parks, there were characteristic species of each one. In Chambi Park, the exclusive existence and high abundance of H. reaumurii in the first habitat (Cham 1) suggests that this species is linked to soil type (sandy loam soil) and halophytes (Linsenmair 1984(Linsenmair , 2007Nasri & Morgan 2005). Armadillidium tunisiense and P. djahizi were the most common isopod species in this park; their presence in large number of specimens, from Cham 2 to Cham 7, suggests that these species are able to use a variety of habitats, as previously reported (Medini & Charfi-Cheikhrouha 2001;Hamaïed & Charfi-Cheikhrouha 2007). In Bouhedma Park, the small species A. lentus was collected mainly in the first habitat (Bouh 1) and sporadically in the two other habitats (Bouh 2 and 3). In fact, for smaller species, the preferences of microhabitats are satisfied by the existence of small places of shelter, regardless of the general characteristics of the habitat (Vilisics et al. 2007(Vilisics et al. , 2011. Thus, the small number of captured individuals indicates the scarcity of favorable microhabitats in the Bouhedma Park. The presence of A. officinalis, reported as xerophilic (Vandel 1962;Messina et al. 2011), was expectable. In North Africa this species was very common in Benghazi, Lybia (Aljatlaoui & Nair 1994), in Egypt (Kheirallah 1980), and also in Morocco (Taiti & Figure 4. Canonical correspondence analysis (CCA) of isopod species, altitude and habitat type. Abbreviations: Bouh 1, 2, 3 and Cham 1, 2, 3, 4, 5, 6, 7habitats in Bouhedma Park and Chambi Park, respectively. Rossano 2015). Porcellio albinus's high abundance in Bouhedma Park suggests its preference for arid and presaharian areas. Some authors have attributed that to its preference for desert regions (Medini-Bouaziz et al. 2016). It has been reported that P. albinus (Linsenmair 1984;Fraj et al. 2008) and H. reaumurii (Shachak & Yair 1984;Baker 2004Baker , 2005 dig burrows. Moreover, for survival and reproduction, H. reaumurii in Chambi Park and P. albinus in Bouhedma Park shelter from the extreme heat and dryness in burrows. Porcellio albinus digs in loose soil and sandy dunes (Medini-Bouaziz et al. 2016, 2017 while H. reaumurii digs in only solid soil (Nasri et al. 1996;Linsenmair 2007).
Regarding the diversity of the isopod communities, relative abundance and species richness values decreased with altitude in both parks. Additionally, the highest diversity values were recorded in open habitats (Bouh 1, 2 and Cham 2). These are likely to be more favorable for isopods than forests with dense vegetation (Cham 3,4,5,6) or dense maquis (Cham 7). However, it was expected that protected areas with dense vegetation should provide more favorable conditions for isopods as they offer abundant food resources. Sfenthourakis et al. (2005) and Hamaïed-Melki et al. (2011) reported that open habitats with a mixture of shrubby vegetation are more favorable for isopods than mono-dominant forests, or very dense vegetation. Moreover, H. reaumurii, present only in Cham 1, exhibited the highest relative abundance. This result supported the habitat specialist hypothesis (Magura et al. 2004). In other cases, species richness was negatively correlated with altitude but no relationship was found between altitude and the Shannon-Wiener index H' (R 2 = 0.1543). These results are, in part, in agreement with previous works on altitudinal gradients and terrestrial isopod diversity which reported a tendency of species richness to decrease with altitude in the mountains of Greece (Sfenthourakis 1992;Lymberakis et al. 2003;Sfenthourakis et al. 2008). Furthermore, the sparsity of species in Cham 3, 4, 5 and 6 (mixed forests) and in Bouh 3 is the result of many factors, namely altitude, dryness, mixed forests with conifers and degree of degradation. The explanation for the low diversity observed in mixed forests of Aleppo and Juniperus is that isopods are generally in need of Ca 2+ to build up their exoskeleton. This suggestion is supported by the high diversity in the open steppe of Stipa and Rosmarinus in Chambi Park and the reforestation of Eucalyptus and Acacia in Bouhedma Park.
The temporal distribution of isopod assemblages studied by analyses of similarity between samples and based on the abundance of species exhibited three major groups whose distribution showed a clear seasonal pattern. Similar distribution patterns for isopods have been reported in other studies (i.e. Hamaïed-Melki et al. 2011;Khemaissia et al. 2012Khemaissia et al. , 2016. The high relative abundance of representative species recorded in spring may be a consequence of the increased surface activity of the adults which are looking for a sexual partner or the gravid females seeking shelter for their descendants in this period (Dangerfield & Hassall 1994). Moreover, the comparison of species composition among habitats based on altitude and plants assemblages allowed for the identification of three habitat groups: the first (Cham 1) dominated by halophilous vegetation, the second composed of the remaining habitats in Chambi Park and the third including the habitats in Bouhedma Park. CCA summarized the overall situation concerning the abundance of specimens for each species and habitat and confirmed the previous two results. Indeed, according to CCA, the altitude, plant assemblage and soil type are the most frequently reported factors in determining the distribution and composition of terrestrial isopods. For instance, several studies concluded that isopods were highly influenced by variations in habitat structure (Davis 1984). Additionally, the community composition of isopods was linked to the degree of openness of the habitat (David et al. 1999).
Compared to Achouri et al. (2008a,b), Hamaïed-Melki et al. (2011), Khemaissia et al. (2012Khemaissia et al. ( , 2016 and Taiti and Rossano (2015), the number of isopod species recorded here, eight species in each park, is lower than that found in humid regions in North Africa. Thus, in the Kroumirie region (northwestern Tunisia) Hamaïed-Melki et al. (2011) found 11 species of terrestrial isopods in the Wadi Moula-Bouterfess area, and in the Berkoukech catchment area, Achouri et al. (2008a) collected 12 species belonging to five families. However, in the supralittoral zones of wetlands, 13 species were recorded (Khemaissia et al. 2012). More recently, studies conducted on the diversity of terrestrial isopods around the wetlands revealed a higher species richness (20 species) (Khemaissia et al. 2016). Further away, in the Oued Laou basin, in the Rif area of northeastern Morocco Achouri et al. (2008b) and Taiti and Rossano (2015) collected 20 and 34 species, respectively. In the Wadi Tahaddart catchment area of northwestern Morocco, Hamaied et al. (2013) identified 18 species belonging to six families. Even farther away, on the north shore of the Mediterranean, in five mountains in Greece, 13 species were collected across 1000 m, 14 species were found at elevations above 1600 m in Tymphi (northwestern Greece), 16 species were found within 1700 m and 11 species were found across 2300 m (Sfenthourakis 1992).
Finally, our results show that Bouhedma and Chambi, located in arid regions, host characteristic species such as H. reaumurii and P. albinus that are the dominant component of the arthropod macrodecomposer guild in certain habitats. These species have developed specific adaptations to the arid environment that enhance their survival and reproduction, such as finding shelter from the extreme heat and dryness in burrows. Usually collected in sites linked to human activities, the occurrence of P. pruinosus and P. laevis is affected by the presence of adjacent motels and other buildings which are favorable habitats for such synanthropic species. Finally, small-sized species such as A. lentus and L. panzerii can find shelter in small crevices and similar places, almost regardless of the general characteristics of the habitat.