It’s All about the Zone: Spider Assemblages in Different Ecological Zones of Levantine Caves

: Caves possess a continuum of ecological zones that differ in their microhabitat conditions, resulting in a gradient of nutrients, climate, and illumination. These conditions engender relatively rapid speciation and diverse assemblages of highly specialised spider fauna. It is unclear, however, how zonation of these caves affects spider assemblage composition and structure. Surveys of 35 Levantine caves were conducted to compare the assemblages of spiders between their different ecological zones. The diverse spider assemblages of these caves differed between the entrance, twilight, and dark zones, with troglophiles and accidental species occupying the cave entrance, endemic troglobites occupying the dark zones, and hybrid assemblages existing in the twilight zones. The progression of assemblage composition and divergence throughout cave zones is suggestive of processes of ecological specialisation, speciation, and adaptation of cave-endemic troglobites in the deepest zones of caves, while cave entrance assemblages are composed of relatively common species that can also be found in epigean habitats. Moreover, the cave entrance zone assemblages in our study were similar in the different caves, while the cave dark zone assemblages were relatively distinct between caves. Cave entrance assemblages are a subset of the regional species pool ﬁltered by the cave conditions, while dark zone assemblages are likely a result of adaptations leading to local speciation events.


Introduction
Subterranean habitats such as caves are home to species with adaptations and preadaptations to darkness and nutrient limitation. Given their unique environments and reduced connectivity, caves resemble islands for obligate troglobite (obligated to life in caves) inhabitants, with little gene flow occurring between populations [1]. These conditions lead to a formation of unique assemblages of highly specialised invertebrates. The cave structure-particularly the characteristics of cave openings-has a direct effect on the abiotic conditions within the cave, such as light intensity, climatic and air conditions, and energy and nutrient influx. These abiotic conditions establish several ecological zones-entrance, twilight, and dark [2][3][4][5]-that are primarily characterised by decreasing illumination with progression into the cave.

Cave Surveys
Arachnids were collected from 35 caves located in Israel and Palestine (West Bank) ( Figure 1) from three cave ecological zones: the cave entrance (inside the cave near its entrance, with high incidence of light); the twilight zone (in the intermediate part of the cave, when present, with low incidence of light); and the inner dark zone (beyond the twilight zone, when present, with no light), as well as outside each cave (but near to its entrance; Figure 2). The caves are distributed along the climatic gradient from the mesic Mediterranean climate in the north and centre of Israel and Palestine (12 caves in each region), to the arid and hyper-arid climates in the south of Israel (11 caves). A blue overlay is given to represent the mean minimum temperature in January (the coldest month of the year) over the 20 years preceding 2021. Dark blue and light blue denote low and high minimum temperatures, respectively. Caves with both twilight and dark zones are denoted by black dots, caves with twilight zones but not dark zones are denoted by grey dots, and caves lacking both twilight and dark zones are denoted by yellow dots. Numbers denote specific cave identities related to numbers assigned by Gavish-Regev et al. [19]. In each cave ecological zone, as well as outside the cave, temperature and illumination were measured. The temperature was measured using a PicoLite 16-K and a single-trip USB temperature logger (FOURTEC, Rosh Ha'Ayin, Israel), with measurements taken once every hour for 74-77 days. The illumination was recorded at the time of each survey using an Extech 401025 Lux Light Meter (Extech, Nashua, NH, USA). The light meter was positioned on the ground until the reading stabilised for a minimum of 1 min. Illumination measurements inside caves ranged between 0 and 420 lux, while measurements outside caves ranged between 60 and 70,000 lux. Temperature and illuminance were also measured outside the caves. Cave length was estimated from cave maps when available (via the Israel Cave Research Centre), or in the field by measuring the distance from the cave opening to the darkest region of the cave accessible during the surveys. Elevation and geological data were provided by the GIS (geographic information system) Centre at the Hebrew University of Jerusalem. Each of the 35 caves was sampled twice according to a specific protocol during 2014 (6 March-6 April 2014 and 22 May-22 June 2014). Due to differences in cave morphology (including microhabitat, fractal shape of the substrates, size, and volume), we standardised our sampling effort by time. Our protocol included a 20 min thorough visual search by one of three experienced arachnologists in 3-10 m long zones using headlamps and UV lights in each of the ecological zones of each cave. In the first visit to each cave, most of the arachnids observed were collected by hand for further identification in the lab. In further visits to each cave, species considered sensitive or common were identified to the species level in situ, without the need to collect specimens (see Gavish-Regev et al. (2021) for additional information and species lists [19]). Voucher specimens were deposited in the National Arachnid Collection in the National Natural History Collections of The Hebrew University of Jerusalem (NNHC, HUJI).

Statistical Analysis
All analyses were conducted in R v4.0.3 [21]. To ascertain the diversity of the ecological zones of the caves surveyed, and the completeness of those surveys, coverage-based rarefaction and extrapolation were carried out, and Hill diversity was calculated as a robust estimate of species diversity [22,23]. This was performed using the "iNEXT" package, with species represented by their frequency of occurrence across samples [22,24]. These Hill numbers include the three most widely-used diversity measures-species richness, Shannon diversity, and Simpson diversity-each varying in their diversity order, q (q = 0, 1, and 2, respectively), which indicates the sensitivity of the measure to relative abundance [24]. Species richness counts species equally, irrespective of their relative abundance; Shannon diversity counts individuals equally, thus representing species proportional to their abundance; and Simpson diversity gives greater weight to the dominant species in the assemblage. These were visually represented by plotting the cumulative diversity (according to all three indices) against the number of detections (occurrences), the cumulative sample coverage (i.e., completeness of sampling) against the number of detections, and the cumulative diversity against the sample coverage. Visual inspection of these figures effectively facilitates the drawing of conclusions regarding the suitability of sampling efforts, and how diversity is structured in cave zones.
Spider assemblages were compared between cave zones using multivariate generalised linear models (MGLMs) via "manyglm" in the "mvabund" package [20], with a Poisson error family and Monte Carlo resampling. To account for the presence of multiple ecological zones in the same cave, any entrance zones in caves containing twilight zones, and any twilight zones in caves containing dark zones, were removed for this analysis; this left 10 entrance zones, 13 twilight zones, and 11 dark zones, all from different caves. Model independent variables included cave zone (entrance, twilight, or dark), zone minimum temperature, cave elevation, geographical region, and all two-way interactions between cave zone and the remaining variables. Geographical region was included, since it was found to be an important factor in spider assemblage composition in a previous study, and may disentangle any regional effects of zone assemblage differences [19]. An alternative analysis, presented in Appendix A, was carried out without excluding zones from the same cave. This enables the reader to inspect results from the complete dataset. This analysis is less conservative and should be viewed with caution; since we were unable to include cave identity as a random effect in these models, we could not discount autocorrelation in the data. The analysis, however, is supportive of the findings of the more conservative analysis presented in the main text.
Coarse differences between assemblages were visualised by non-metric multidimensional scaling (NMDS) via metaMDS in the "vegan" package [25], with Bray-Curtis distance in 2 dimensions and 999 tries. For visualisation of the effect of categorical variables against the NMDS, spider plots were created using "ordispider" with "ggplot2" [26]. For visualisation of the effect of continuous variables against the NMDS, surf plots were created with scaled coloured contours using "ordisurf" with "ggplot2".
Pairwise co-occurrence analysis was carried out to identify spider species that occurred together more, or less, than expected by chance. The "cooccur" function in the "cooccur" package [27] was used to calculate the expected frequencies of co-occurrence between each pair of taxa via null models, which were then compared against the observed patterns of co-occurrence to identify deviations from random. Nestedness-the ordered loss of species-was assessed across cave zone spider assemblages using the binary-matrix nestedness temperature calculator (BINMATNEST; [28,29]. A binary presence-absence matrix was used to calculate "temperature" (0-100 • C; the deviation from perfect nestedness, represented by 100 • C). The nestedness of the binary matrix was compared against 100 matrices generated randomly using null models [30].

Spider Assemblage Diversity in Cave Zones in the Southern Levant
Across the surveys, a total of 1054 spiders were collected or identified in the field, belonging to 62 species and morphospecies across 38 genera and 22 families. Many species were found across all three cave zones, but some in only one or two ( Figure 3). The cave assemblages were highly diverse (Hill richness = 63.00 ± 7.79; Hill-Shannon diversity = 32.21 ± 4.92; Hill-Simpson diversity = 19.53 ± 4.38; Figure 4) and largely composed of rare species, with many species found only in a single cave. The surveys identified an estimated 89.57% (± 2.8%) of the total spider diversity in caves ( Figure 4). Limited nestedness was observed across cave zones (3.98 • C, vs. 16.79 ± 2.79 • C in the null models; p < 0.05).

Zone Assemblage Comparison
Specific spider assemblages were significantly related to cave zones (MGLM: Dev = 426.0, d.f. = 31, p = 0.001; Figure 5), but also to the interaction between zone and minimum temperature (MGLM: Dev = 124. 8 Fourteen species were significantly associated with cave zones or interactions between cave zone and other variables; of these, six spider species were significantly associated with interactions between cave zone and temperature, cave length, or geographical region (Table 1, Figure 7). A further five spider species were significantly associated with minimum temperature, elevation, or geographical region (Table S1, in supplementary material). Specifically, Artema nephilit (Aharon, Huber, and Gavish-Regev, 2017), Filistata insidiatrix (Forsskål, 1775) (interacting with minimum temperature), and Tegenaria pagana (C.L. Koch, 1840) (interacting with geographical region, minimum temperature, and elevation) were almost equally common in both twilight and entrance zones; Filistata sp., Holocnemus pluchei (Scopoli, 1763), Loxosceles rufescens (Dufour, 1820) (interacting with geographical region and minimum temperature), Micronetinae sp. (interacting with minimum temperature), Oecobius sp., Steatoda triangulosa (Walckenaer, 1802) (interacting with geographical region and minimum temperature), and Tegenaria angustipalpis (Levy, 1996) in entrance zones; Pholcus sp. in twilight zones; Tegenaria sp. from the Galilee (interacting with minimum temperature) almost equally in both twilight and dark zones; and Hoplopholcus cecconii (Kulczynski, 1908) and Tegenaria sp. from Zavoa cave in dark zones ( Figure 7). Table 1. Significant univariate MGLM results for the 14 species with significant associations with cave zones, or with interactions between cave zones and other variables; deviance and probability are given for each. Other significant variables are listed, as well as interactions between zone and those variables. For the other associations, the nature of the association is given as + or − for positive or negative associations, respectively, and "N", "C", and "S" are given for prevalence in north, central, or southern geographical regions, respectively. The "category" column denotes whether that species is troglobite, troglophile, or accidental (see Gavish-Regev et al., 2021 [19] for the assignment of species to categories). The "dominant zone" column denotes the ecological zones(s) in which the species was mostly commonly found. The abundance of these spiders across different cave zones is visualised in Figure 6. For spiders with significant associations with interactions between cave zone and other variables, but not with the cave zones themselves (n = 3), the deviance and probability are given in italics.  Figure 5. Spider plots of assemblages by cave zone derived from non-metric multidimensional scaling. Colours denote cave zones (black, grey, and yellow denoting dark, twilight, and entrance zones, respectively). Axes represent a two-dimensional variation in spider assemblages. Each smaller point in both images represents a single assemblage present in a specific cave and a specific zone, with the distance between them indicating their dissimilarity (i.e., proximate points are similar, distant points are dissimilar). Average species coordinates are overlaid in Figure S1. Stress = 0.1002019. (A) Points are joined by the centroids of assemblages in each group (larger points = mean coordinates in that group); a high degree of overlap is observable, particularly between twilight and entrance assemblages, but dark zone assemblages are more distinct and highly variable. (B) Joined points belong to the same cave; these are typically more similar (i.e., points are closer) than those of other caves, with the exception of many of the dark zone assemblages (generally around the outside of the central cluster of points), which show substantial variation. Figure 6. Surf plots of assemblages by (A) minimum cave temperature and (B) cave elevation, both derived from non-metric multidimensional scaling. Each point represents a single assemblage present in a specific zone of a specific cave, with the distance between them indicating their dissimilarity (i.e., proximate points are similar, distant points are dissimilar). Point colours denote cave zones (black, grey, and yellow denoting dark, twilight, and entrance zones, respectively). Contour colours denote (A) temperature (purple-yellow denoting low-high temperatures) and (B) elevation (dark-light blue denoting low-high elevation). Axes represent a two-dimensional variation in spider assemblages. Most of the assemblages that existed in cooler caves appear to be more similar than those in warmer caves. A complex nonlinear trend is observable between assemblage composition and elevation, with assemblages in intermediate-elevation caves appearing more similar. Average species coordinates are overlaid in the spider plot constructed with the same sample coordinates in Figure S1. Stress = 0.1002019. from Zavoa cave (entrance)), but due to subsampling of the data and the use of only one zone per cave for the MGLM, these data are absent from this graph (see Figure 3 and Appendix A for the complete dataset, and Materials and Methods for more details regarding the MGLM). Dark, twilight, and entrance zone abundances are denoted by black, grey, and yellow, respectively. Occurrence is calculated as the percentage of surveys of a zone type in which at least one individual of each species was found.

Co-Occurrence of Spiders in Cave Zones
Of the 1953 possible pairwise species co-occurrence combinations, 1872 (95.85%) were not carried forward for co-occurrence analysis, since they were expected to co-occur less than once. Of the 81 pairs analysed, 68 co-occurred randomly, and 13 (16%) non-randomly ( Figure 7, Figure S2, Table 2). Of these non-randomly co-occurring pairs, 11 co-occurred more than expected by chance-Chaetopelma sp. with Loxosceles rufescens, F. insidiatrix with H. pluchei, F. insidiatrix with L. rufescens, F. insidiatrix with Tegenaria pagana, L. rufescens with H. pluchei, H. pluchei with Micronetinae sp. 1, H. pluchei with Steatoda triangulosa, H. pluchei with U. plumipes, L. rufescens with Micronetinae sp. 1, L. rufescens with T. pagana, and Pholcus sp. with T. pagana-and 2 less than expected by chance: A. nephilit with Hoplopholcus cecconii, and Artema nephilit with S. triangulosa (Table 2, Figure 8). Table 2. Species co-occurrence across cave zones. The probability corresponds to the probability that the respective species co-occur more, or less, than expected (listed as "relationship", with "+" and "−" denoting positive and negative co-occurrences, respectively). Negatively co-occurring species are in bold.  8. Species co-occurrence matrix for those species identified from cave spider assemblages from distinct cave zones. Yellow, grey, and blue points denote significantly negative, random, and significantly positive co-occurrences, respectively.

Discussion
We have shown that the different ecological zones of caves host unique assemblages of spiders, including many rare species. The diverse cave spider assemblages within these environments are structured by a high degree of endemicity, with many species occurring only in a single cave or cave zone. The fact that this study identified an estimated 89.57% of the total assemblages in these caves, based on rarefied sample coverage, was considered promising and unsurprising given the endemicity of many of the extant fauna; each additional cave sampled is likely to host a unique assemblage, possibly including its own endemic species-especially where those caves contain multiple distinct ecological zones. Indeed, such caves have recently yielded taxa new to Israel and new to science [31][32][33][34].
Cave assemblages show a graded pattern of dissimilarity, with dark and entrance assemblages appearing the most dissimilar, and twilight assemblages existing as an intermediate. The relative similarity of entrance assemblages between different caves, compared with the dissimilarity of different dark assemblages (as exemplified in Figure 5), suggests that dark zone assemblage composition is not directly derived from the species pool at the cave entrance, but may have diverged from ancestral entrance assemblages over time. This is supported by a low nestedness of the assemblages within the caves, contrary to our original hypothesis. These patterns of dissimilarity and low nestedness may be due to the existence of relict species in some of the caves sampled, or to speciation within caves as a consequence of low gene flow [1]. Weak gene flow could result in a transition from entrance assemblage species to distinct dark zone troglobite species through adaptive radiation [1,35]. Following the adaptive shift hypothesis (ASH; [36]), troglobitic dark zone spider assemblages could have developed from entrance and twilight assemblages, but differentiated from them over time through these speciation events. These assemblages may have included species with pre-adaptations to caves, such as those adapted to shallow depressions or hollows [37], thus facilitating further adaptation to dark zones. Following the climatic relict hypothesis (CRH), dark zone assemblages may include ancestral or relict species that did not possess pre-adaptations, but survived in caves after changes in climatic conditions rendered the epigean habitat unsuitable for them [38]. The two hypotheses are not mutually exclusive, however, and our study-an ecological snapshot in time-was not designed to test these hypotheses. Further study is required in order to fully elucidate these evolutionary processes.
Our models explaining the assemblage structures in different caves revealed interactions between ecological zone and other variables such as minimum temperature, elevation, and region, which may reflect the particular conditions found in caves hosting multiple zones. The differences in assemblage structure between zones were mostly affected by those taxa individually associated with different zones (highlighted in Table 1). These species belong to just a few families (Agelenidae, Filistatidae, Linyphiidae, Pholcidae, Sicariidae, and Theridiidae). The two most common families among these-Pholcidae and Agelenidae-had consistent and generally opposite cave zone occupancies. Pholcidae (troglophile species) were most often prevalent in the entrance, with a consistent presence in the twilight zone. One species (Hoplopholcus cecconii), however, occurred often in the dark zone, as seen in the data used in the MGLM (Table 1, but for complete data, see Figure 3 and Appendix A). When considering these additional observations, H. cecconii was most abundantly found in the twilight and entrance zones (Appendix A), indicating that this species was in fact found in all three zones, but of the same caves, and was exclusively found in caves large enough to contain all three ecological zones.
Agelenidae appeared across all three zones, with four species that differed in abundance between zones (T. angustipalpis, T. pagana, and two undescribed species found in the Galilee caves and in the Zavoa cave) and occurred in different permutations of cave zones, indicating a high degree of habitat generalism. Most of the species that differed significantly in their zone occupancy were troglophiles. In the MGLM, only three of these troglophiles were not found in the cave entrance: Pholcus sp., Hoplopholcus cecconii, and Tegenaria sp. Zavoa. All of these species, however, were found in entrance zones, which were excluded from the MGLM because they were from the same cave as a zone already represented in the model (as explained in the Materials and Methods section and Appendix A). The only troglobite species associated with a specific cave zone was Tegenaria sp. Galilee, found exclusively in the dark and twilight zones. Congeneric species such as Tegenaria pagana occupied both the cave entrance and twilight zones, indicating that zone occupancy is not taxon-specific, with species within a genus differing in their adaptation to cave systems. The only accidental species with significant association with a specific cave zone (Oecobius sp.) also occurred almost exclusively in the entrance, suggesting that this accidental species may lack the pre-adaptations necessary to thrive in twilight and dark zones.
The 13 co-occurrence relationships identified as deviations from random highlight some consistencies in the caves' regional spider assemblage structuring. The positive co-occurrence relationships prominently include the most common species-namely, Holocnemus pluchei (5/13 significant co-occurrences, found in 28% of zone surveys) and Filistata insidiatrix (3/13 significant co-occurrences, found in 31% of zone surveys). These cooccurrences may have been more frequent than expected simply because these species are widespread and common representatives of the regional cave assemblages in all but the dark zones.
Only two co-occurrences were identified as being significantly less common than expected: Artema nephilit with Hoplopholcus cecconii, and A. nephilit with Steatoda triangulosa. Artema nephilit, which was encountered in 30% of zone surveys, was not present in any positive co-occurrence relationships, despite being in both negative co-occurrence relationships. Artema nephilit was found mainly in warmer caves (Appendix A), and was often among a very small number of species found in hot and dry caves (Gavish-Regev, personal observations). The specific negative co-occurrence of A. nephilit and Hoplopholcus cecconii is probably a result of different habitat requirements or, as highlighted by the MGLM results (Table 1), differences in elevation. Additionally, A. nephilit was found mainly in the Rift Valley to the east, whereas H. cecconii occurred mainly toward the Mediterranean region in the west; this would not be identified by the MGLMs, despite their inclusion of region, since the analysis was focused on a north/central/south division. Where the distributions of A. nephilit and H. cecconii do overlap, in the Karmel, A. nephilit was found exclusively in the warmer, south-facing cave, while H. cecconii was found only in the cooler, north-facing cave, just a few meters away.
The distribution of spiders in caves can be affected by niche dynamics and competition for resources, with dominant species forcing species with pre-adaptations away from entrance zones and into the relatively prey-poor twilight and dark zones [7]. This "step back effect" may be a significant driver of diversity in the innermost regions of caves [4]. Such effects could theoretically be responsible for the negative co-occurrence relationships between species such as A. nephilit and Steatoda triangulosa; although in this case these two species were not present together in the same caves, they occupy a similar spatialand likely trophic-niche in the caves, and their negative association could result from competitive exclusion of one another. A similar process may have occurred in Pholcidae, since most caves contain only a single species of Pholcidae [39]. Broader analyses of cave spider prey communities, interaction networks, and trophic ecology may further disentangle additional spatial patterns in, and drivers of, spider cave occupancy.

Conclusions
Cave ecological zones host divergent yet diverse assemblages of spiders. This diversity largely reflects an adaptive process of speciation, from relatively unspecialised and widely distributed entrance zone assemblages, to highly specialised dark zone assemblages comprising many endemic species. Our results indicate that differences in assemblages between cave zones are further affected by temperature, elevation, and geographical region, but the patterns of co-occurrence in these assemblages also suggest that competition is a driver of cave spider spatial dynamics. The high degree of endemicity in these systems lends credence to their importance as hotspots of biodiversity, but their relative climatic stability does not render them impervious to the effects of climate change and other drivers of species loss [40]. A greater understanding of the drivers of diversity in cave spider populations is necessary in order to ensure sufficient knowledge to protect these important and unique contributors to the global arachnofauna. zones from the same cave, and is presented herein, but should be viewed with caution and scepticism. Since MGLMs do not allow the inclusion of random factors-in this case, cave identity-this analysis will be affected by the autocorrelation of the data collected from the same caves. We have decided to present it nevertheless, in order to highlight the overall consistency of the analyses, but also the effects of some of the additional observations (e.g., the large number of observations of Pholcus sp. in entrance zones).

Appendix A1. Materials and Methods
Spider assemblages were compared between cave zones using multivariate generalized linear models (MGLMs) via "manyglm" in the "mvabund" package [20], with a Poisson error family and Monte Carlo resampling. Model independent variables included cave zone (entrance, twilight, or dark), zone minimum temperature, cave elevation, geographical region, and all two-way interactions between cave zone and the other variables. Ten species were found to have significant associations with cave zones; three of these species were also significantly associated with temperature and geographical region (Table A1, Figure A1). Specifically, Artema nephilit (Aharon, Huber, and Gavish-Regev, 2017) and Pholcus sp. were most common almost equally in both twilight and entrance zones; Filistata insidiatrix (Forsskål, 1775) (depending on an interaction between zone and geographical region), Filistata sp., Holocnemus pluchei (Scopoli, 1763) (depending on an interaction between zone and geographical region), and Uloborus plumipes (Lucas, 1846) in entrance zones; Tegenaria sp. from cave in Judea in twilight zones; Tegenaria sp. from caves in the Galilee (depending on an interaction between zone and temperature) almost equally in both twilight and dark zones; and Steatoda sp. from Zavoa cave and Tegenaria sp. from Zavoa cave in dark zones ( Figure A1). Figure A1. Bar plot of the percentage of surveys in which the 10 species with significant associations to cave zones were identified across cave zones. Dark, twilight, and entrance abundances are denoted by black, grey, and yellow, respectively. Table A1. Significant univariate MGLM results for the 10 species with significant associations to cave zones; deviance and probability are given for each. Other significant variables are listed, as well as interactions between zone and those variables. For the other associations, the nature of that association is given as "+" or "−" for positive or negative associations, respectively, and "N", "C", and "S" are given for prevalence in north, central, or southern geographical regions, respectively. The "category" column denotes whether that species is troglobite, troglophile, or accidental (see Gavish-Regev et al., 2021 [19] for the assignment of species to categories). The "dominant zone" column denotes which ecological zones(s) a species was mostly commonly found in. The abundance of these spiders across different cave zones is visualised in Figure A1.