Optimizing speleological monitoring efforts: insights from long-term data for tropical iron caves

Understanding the factors underpinning species abundance patterns in space and time is essential to implement effective cave conservation actions. Yet, the methods employed to monitor cave biodiversity still lack standardization, and no quantitative assessment has yet tried to optimize the amount and type of information required to efficiently identify disturbances in cave ecosystems. Using a comprehensive monitoring dataset for tropical iron caves, comprising abundance measurements for 33 target taxa surveyed across 95 caves along four years, here we provide the first evidence-based recommendations to optimize monitoring programs seeking to follow target species abundance through time. We found that seasonality did not influence the ability to detect temporal abundance trends. However, in most species, abundance estimates assessed during the dry season resulted in a more accurate detection of temporal abundance trends, and at least three surveys were required to identify global temporal abundance trends. Finally, we identified a subset of species that could potentially serve as short-term disturbance indicators. Results suggest that iron cave monitoring programs implemented in our study region could focus sampling efforts in the dry season, where detectability of target species is higher, while assuring data collection for at least three years. More generally, our study reveals the importance of long-term cave monitoring programs for detecting possible disturbances in subterranean ecosystems, and for using the generated information to optimize future monitoring efforts.


INTRODUCTION
Quantifying long-term changes in abundance of cave-dwelling organisms and identifying indicator species, reflecting the health status of subterranean ecosystems, are among the fundamental research goals of modern subterranean conservation biology (Mammola et al., 2020). For instance, the lack of knowledge about the factors underpinning abundance patterns in space and time are among the main impediments to the effective protection of cave fauna (Cardoso et al., 2011). Long-term studies in caves are scarce (Di Russo et al., study are located in two highlands known as Serra Norte and Serra Sul (Fig. 1). These two regions harbor banded ironstone formations known as cangas, unique campo rupestre ecosystems resembling mountain savannas (Zappi et al., 2019), and one of the world's largest deposits of iron ore (Poveromo, 1999).

Database
We used data generated by independent environmental consulting companies, so our study did not involve any field work. Vale S. A., a mining company, began operations in the region more than two decades ago (Souza-Filho et al., 2019), and has conducted numerous caves surveys over the last years as part of a large monitoring program related to environmental licensing processes. All surveys where authorized by Instituto Brasileiro do Meio Ambiente e dos Recursos Naturais Renováveis (IBAMA), under licenses ABIO 455/2014 (Projeto Ferra Carajás S11D n 02001.000711/2009-46) and ABIO 639/2015 (Projeto Ferro Serra Norte-Estudo Global das Ampliações das minas N4 e N5 n 02001. 002197/2002-15). We compiled the data generated in these surveys to collect information from 33 target taxa across 95 caves, surveyed between August 2015 and September 2019. The selection of species included in these monitoring programs was based on the following criteria, as stated in environmental assessment reports (ATIVO AMBIENTAL, 2019; BRANDT, 2019): large body size and easy to identify in the field, abundant and showing a wide distribution range, resolved taxonomic classification (at least to the morpho-species level), and short life cycles allowing the rapid detection of changes in population dynamics (see Table 1 for the full list of target taxa and their ecological classification). All the selected species were actively surveyed during each field trip, so absences represent true absences rather than missing data. In each cave, the absolute abundance of each target taxa was quantified at least once during the rainy and the dry season, and sometimes multiple times in one year. Sampling was performed  (litter, logs, carcasses, guano, etc.). Animals were collected with the aid of tweezers and brushes, and all individuals found in each cave were counted to estimate abundance per species, as performed in other studies (Silva, Martins & Ferreira, 2011;Ferreira et al., 2015;De Bento et al., 2016;Paixão, Ferreira & Paixão, 2017;Ferreira & Pellegrini, 2019;Souza-Silva, Iniesta & Ferreira, 2020).

Environmental conditions and landscape metrics
External and internal environmental conditions were monitored during the entire period across caves. Monitored variables included the deviation in average bimonthly rainfall in relation to the expected from a 20-years series (in mm, retrieved from small weather stations located in nearby mines S11D e N4E), and mean internal temperature ( C) on the date of the surveys (retrieved from portable data loggers placed in the most distant location from cave entrances). We also recorded the Area (meters 2 ) of each studied cave as an additional internal condition widely known to influence biodiversity patterns in these ecosystems (Jaffé et al., 2016(Jaffé et al., , 2018. Using 30 m resolution land-cover maps from 2015 to 2019 , we then quantified a suit of landscape metrics, including the proportional amount of forest, canga and mining land covers surrounding caves, and topographic distance to the nearest mine (see details in Table S1). These were all calculated at two different spatial scales (circular buffers with 500 and 1,000 m radius), using the R packages landscapemetrics (Hesselbarth et al., 2019) and TopoDistance (Wang, 2020). Two of these metrics directly captured possible disturbance of subterranean environments that could account for changes in the abundance of the studied species: mining cover and distance to the nearest mine.

Assessing drivers of community composition across caves
Aiming to quantify how environment, cave and landscape variables influenced overall community composition, we ran a partial redundancy analysis (RDA) controlling for differences between both highlands (Serra Norte and Serra Sul), using the vegan package (Oksanen et al., 2019). The community composition matrix containing relative abundances for each taxa was used as response variable and predictor variables included year, season, microclimate and landscape metrics (Legendre & Legendre, 1998). The highland where caves were located was specified as a conditional variable on the model to control for the effect of cave´s geographical location. Microclimate and landscape variables were standardized, community composition was Hellinger-transformed, and permutation tests were used to assess significance of marginal effects (Legendre & Legendre, 1998).
Assessing the influence of seasonality on the detection of temporal abundance trends One of the main goals of cave monitoring programs was to assess changes in species abundance over time, and thereby identify species with declining or increasing populations in a particular cave. To understand how seasonality influenced the detection of abundance trends over time, we ran linear models containing the total number of observed individuals as the response variable and the interaction between sampling date and season. If seasonality influences temporal abundance trends, we would expect to find significant interaction terms. No significant interactions, on the other hand, would indicate that the trends can be detected regardless of the season when the surveys where performed.
To prevent overfitting, linear models were ran for taxa and caves represented by at least five surveys in each season (final sample size was 16 taxa and 50 caves). Given the large number of models we used the Benjamini & Hochberg approach to adjust p-values, employing the p.adjust function from the stats R package (R Development Core Team, 2020).
Assessing the influence of sampling effort on the detection of temporal abundance trends Given the extensive field exposure of people and elevated costs associated with cave monitoring programs, it is important to quantify how the sampling effort influences the detection of temporal abundance trends. To do so we compared linear model coefficients of models fitted with the full dataset with those of models fitted with reduced datasets. We first split the data by season and ran linear models containing the total number of observed individuals as the response variable and sampling date as predictor. In these full models we included observations for all sampling dates, and excluded taxa and caves represented by less than three surveys per season. We then ran linear models on data subsets containing a reduced number of observations (ranging between two and the maximum number of sampling dates found in each cave and taxa). For each data subset containing a given number of observations (surveys) we performed ten random samplings without replacement, to ensure the sampling of different sampling dates. Finally, we compared coefficients from full models with those of subset models using root mean squared error (rmse), implemented through the rmse function from the Metrics R package (Hamner & Frasco, 2018). Lower values of rmse indicate more similar model coefficients.

Identifying disturbance indicator species
Given the life history variation between species and their different susceptibility to habitat disturbance, it is essential to identify indicator species that show a rapid response to disturbance in order to optimize monitoring programs. By focusing on these indicator species, monitoring programs could survey caves more efficiently, thereby making resources available to study more caves or other aspects of cave biodiversity requiring attention. Here we tried to identify disturbance indicator species by assessing the relationship between disturbance metrics and species abundance patterns. We first modeled patterns of relative abundance across all caves, using the function manyglm from the R package mvabund (Wang et al., 2012). It uses a multivariate generalized linear model (GLM) to make inferences by fitting separate GLMs to a common set of explanatory variables, and testing significance through resampling-based hypothesis testing. We ran negative binomial GLMs containing abundance as the response variable and sampling season nested in year, distance to mine and mining cover as predictor variables.
Significance p-values were calculated using 999 resampling iterations via PIT trap resampling, adjusted for multiple testing using a step-down resampling procedure (Wang et al., 2012). We then used univariate coefficient estimates and significance for individual species, to identify specific responses to disturbance metrics (distance to mine and mining cover). We then assessed the relationship between disturbance metrics and temporal trends in species abundance within each cave. To do so we ran linear models containing the total number of observed individuals as the response variable and sampling date as predictor, excluding taxa and caves represented by surveys spanning less than three years (some caves where surveyed multiple times in a single year but these where only included in this analysis if surveys spanned at least three different years). We used the model coefficients for each species at each cave, representing temporal abundance trends (positive coefficients showing an increase and negative coefficients a decrease in abundance through time), to run a second set of linear models regressing temporal abundance trends on disturbance metrics. These second set of models thus contained as response variable the model coefficients representing temporal abundance trends for each species at each cave, and distance to mine and mining cover (at different spatial and temporal scales) as predictors. To prevent overfitting we excluded species represented by less than ten coefficients (caves), and only constructed models containing a single predictor. We then ran likelihood-ratio tests, where we compared each model with a null model containing no predictors, and selected those predictor variables resulting in a significant increase in the model's log-likelihood. Finally, we retrieved and plotted coefficients and p-values for these best-fitting models. All data and R scripts are available in GitHub (https://github.com/rojaff/cave_monitoring).

RESULTS
Overall community composition was weakly influenced by seasonality, cave size, environmental conditions, and the composition and configuration of landscapes The table shows F-statistics and p-values from permutation tests (adjusted r 2 = 0.13). Significance is highlighted as * (p < 0.05), ** (p < 0.01) and *** (p < 0.001).
surrounding caves, as more than 87% of variance in community composition remained unexplained by these factors (Table 2). Seasonality did not influence the ability to detect species abundance trends over time, since the interaction effect between sampling date and season was not significant in any taxa nor cave (Fig. 2). Increasing the number of samples resulted in more similar model coefficients between full and subset models, and root mean squared errors usually stabilized after three surveys (Fig. 3). However, in most species the dry season datasets allowed a more accurate detection of temporal abundance trends, as revealed by lower root mean squared errors (Fig. 3).
Whereas relative abundance was associated to at least one disturbance metric in 22 species (Fig. 4), temporal trends in abundance were found associated with disturbance metrics in only five species (Fig. 5). Overall, two taxa displayed consistent responses across effects, which makes them potential indicator species for cave monitoring programs: The troglobiont Charinus ferreus, which appeared negatively affected by disturbance, and a species belonging the Theraphosidae family, which seem to be favored by disturbance (Table 3).

DISCUSSION
By analyzing abundance measurements for 33 target taxa surveyed across 95 caves along four years, we found that overall community composition was weakly influenced by seasonality, cave size, environmental conditions, and the composition and configuration of landscapes surrounding caves. Furthermore, our results show that seasonality did not influence the ability to detect abundance trends over time. However, in most species, abundance estimates assessed during the dry season resulted in a more accurate detection of temporal abundance trends, and at least three surveys were required to identify global temporal abundance trends. Finally, we identified a subset of species that could potentially serve as short-term disturbance indicators, some showing consistent responses in different analyses. Subterranean communities have been shown to be affected by seasonality, environmental conditions, cave characteristics, and the structure of surrounding landscapes (Simões, Souza-Silva & Ferreira, 2015;De Bento et al., 2016;Mammola & Isaia, 2018;Salvidio et al., 2019;Pellegrini, Faria & Ferreira, 2020;Rabelo, Souza-Silva & Ferreira, 2020). However, our results reveal that overall community composition was only weakly influenced by these factors, as our model explained merely 13% of total variation in community composition (Table 2). In contrast, previous work have found that cave morphology, microclimate, cave depth, and sampling date explain up to 50% of the variation in community structure in limestone and marble caves (Tobin, Hutchins & Schwartz, 2013;Lunghi, Manenti & Ficetola, 2014). Our results thus suggest that other factors, not considered in our analyses, play an important role structuring subterranean communities of iron caves. Inter-specific interactions, for instance, are known to have a profound influences on community structure (Ferreira & Martins, 1999;. Alternatively, biological samples collected in iron caves may not capture the dynamics of the entire subterranean habitat, comprised by a network of fissures and voids and traditionally referred to as Milieu Souterrain Superficiel (MSS) (Culver & Pipan, 2014;Mammola, 2018). For instance, most of the surveyed caves were larger than 5 × 5 m (Fig. S1), so they did not represent suitable sampling sites for the MSS . Even though seasonality affected overall community composition, it did not influence the ability to detect species abundance trends over time. External climatic conditions are increasingly attenuated at higher cave depths (Tobin, Hutchins & Schwartz, 2013), Rhipidomys sp.
Taxa showing consistent responses (highlighted in bold) are suggested as short-term disturbance indicators. The best sampling season (according to Fig. 3), is indicated for each taxa.
so species occurring in the inner portions of caves appear to have life cycles decoupled from external seasons, whereas species inhabiting the outermost portions of caves seem to be more strongly affected by seasonality (Di Russo et al., 1997;Gunn, Hardwick & Wood, 2000;Bichuette & Trajano, 2003;Ferreira et al., 2015;Lunghi, 2018). Recognizing the impact of seasonality on species detection, the current Brazilian legislation stipulates that cave biodiversity surveys need to comprise at least two sampling events, one during the dry and one during the rainy season (MMA, 2017).
It is worth emphasizing that these sampling requirements targeted a more accurate estimation of species richness, but not the continuous monitoring of focus species in time. Two sampling events are likely insufficient to obtain reliable species richness estimates for highly diverse caves (Auler & Piló, 2015;Wynne et al., 2018), so some authors have argued for the estimation of optimal sample sizes based on species accumulation curves (Trajano & Bichuette, 2010;Trajano, 2013). Our results provide the first evidence-based recommendations to optimize sampling efforts of monitoring programs seeking to assess target species abundance through time. Specifically, our findings suggest that monitoring efforts aiming to detect changes in abundance through time do not need to sample during two different seasons each year (Fig. 2). Sampling efforts of such monitoring programs could thus be optimized by performing more focused surveys and by surveying a larger number of caves during the same period each year. Importantly, restricting sampling to a single season could substantially attenuate the negative impact of cave visitation by researchers on subterranean communities Pellegrini & Lopes Ferreira, 2012;Bernardi, Souza-Silva & Ferreira, 2010).
Although the composition and spatial distribution of subterranean communities can remain constant over periods of several years (Salvidio et al., 2019), our results suggest that sampling during at least three years is necessary to detect temporal changes in abundance patterns in most of our focus species (Fig. 3). We note that our dataset only spans a period of four years (although some caves were sampled multiple times during the same season/year), so it cannot capture longer temporal changes in abundance. We also caution that these results cannot be generalized to all subterranean fauna, as different life histories and generation times will ultimately determine how fast these organisms respond to disturbances (Ferreira, 2005;Culver & Pipan, 2019). Sampling in different seasons did not influence the ability to detect general abundance trends over time, but the dry season datasets allowed a more accurate detection of temporal abundance trends in most species. These results suggest higher detection probabilities in the dry season for the subset of species where RMSE curves show a steeper decrease during the dry season (Fig. 3). Interestingly, this was the case for the troglobitic amblypygid Charinus ferreus, a species that is difficult to detect like other troglobionts (Wynne et al., 2018;Lunghi, 2018). Our results thus suggest that monitoring programs focusing on terrestrial subterranean fauna from our study region could concentrate sampling activity in the dry season, where most species seem to be easier to detect. Likewise, our findings highlight the importance of implementing long-term monitoring efforts spanning at least 3 years.
The concept of indicator species in ecosystem management relies on the idea of identifying taxa responsive to environmental change, that could inform policies, protocols, and best practices (Carignan & Villard, 2002). Such environmental indicators (McGeoch, 1998) seek to provide cost and time effective guidelines to address pressing conservation issues, such as those faced by large-scale mining projects (Sonter, Ali & Watson, 2018). Assessing the response of subterranean fauna to anthropogenic disturbance nevertheless requires access to long-term cave monitoring datasets, which are remarkably rare for tropical caves (McGeoch, 1998;Carignan & Villard, 2002;Mammola et al., 2020). Here we identified 20 taxa where overall abundance responded to cave disturbance, and five where temporal abundance trends where associated with disturbance. Only two taxa displayed consistent responses across effects, which makes them candidate indicator species for cave monitoring programs: Charinus ferreus and a species belonging the Theraphosidae family (Table 2). Both are arachnids, a group that was recently identified as biodiversity indicator for iron caves (Trevelin et al., 2019). Being a top predator restricted to cave ecosystems, the first species is a well-known troglobitic Amblypygi (De Lao Giupponi & De Miranda, 2016). Its strong and consistent response to disturbances (Figs. 4 and 5) suggest the species is associated with pristine and undisturbed ecosystems, which makes it an ideal disturbance indicator. Theraphosidae spiders, on the other hand, are sedentary sit-and-wait predators from the epigea, rarely occupying subterranean environments for reproduction or shelter (Fonseca- Ferreira, De Zampaulo & Guadanucci, 2017). Our results suggest that they apparently benefit from disturbance to opportunistically colonize caves, or alternatively, that disturbances in the surrounding external habitats are forcing them to look for shelter inside the caves. The species nevertheless awaits formal taxonomic description, which currently limits its usefulness as an indicator species.
Effect sizes of disturbance on overall abundance and temporal abundance trends where generally small, suggesting that some effects could have remained undetected because they would require sampling over longer time periods (Di Stefano, 2001;Legg & Nagy, 2006). For instance, the ability to detect trends in tropical bat population abundance was shown to be dependent on the duration of the monitoring efforts, and only long programs (>20 years) showed sufficient statistical power to reliably detect abundance trends (Meyer et al., 2010). This could explain why some of our focus species did not exhibit coherent responses across analyses, like the troglobionts Pyrgodesmidae sp. and Escadabiidae sp., or opportunistic colonizers like the anuran Leptodactylus pentadactylus or Pristimantis fenestratus. Although empirical evidence from long-term cave monitoring efforts focusing on invertebrates is scarce (Faille, Bourdeau & Deharveng, 2015;Cajaiba, Cabral & Santos, 2016;Owen et al., 2016), our results suggest that longer monitoring efforts are needed to detect disturbance responses in most cave-dwelling species.

CONCLUSIONS
Our study reveals the importance of long-term cave monitoring programs for detecting possible disturbances in subterranean ecosystems, and for using the generated information to optimize future monitoring efforts. Results show that iron cave monitoring programs implemented in our study region could focus sampling efforts in the dry season, where detectability of target species is higher, while assuring data collection for at least three years. Charinus ferreus was identified as the most promising short-term disturbance indicator species.

Data Availability
The following information was supplied regarding data availability: All data and R scripts are available in GitHub: https://github.com/rojaff/cave_monitoring

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.11271#supplemental-information.