Earthquake impacts on microcrustacean communities inhabiting groundwater-fed springs alter species-abundance distribution patterns

Earthquakes are important natural events, yet their impacts on animal communities are poorly known. Understanding earthquake impacts on groundwater communities is essential to assess their resilience and hence to perform conservation actions. We investigated how a 6.3 Mw earthquake that occurred in 2009 altered the community structure (diversity, evenness, dominance, species abundance distributions and beta-diversity) of microcrustaceans (Crustacea Copepoda) inhabiting springs fed by the Gran Sasso Aquifer (Central Italy). Sampling was done in low-discharge (1997), high-discharge (2005), and post-seismic (2012) hydrological years. Stygobites (obligate groundwater species) and non-stygobites (non-obligate groundwater species) showed different patterns. A high-water discharge in 2005 altered abundance patterns of non-stygobites. The earthquake re-established former abundance patterns. Stygobites were less affected by high-water discharge in 2005, and showed strong increases in diversity and evenness after the earthquake. This effect was due to the fact that the earthquake induced a strong population decline of previously dominant stygobites (especially of Nitocrella pescei) in the aquifer, and subsequently at the main spring outlets, thus allowing a more equitable species-abundance distribution. These results highlight the importance of considering species ecology to understand the effects of a significant earthquake event on animal communities.

Comparison of pre-and post-seismic communities for stygobites. For the stygobites, the two pre-seismic years differed only for the Berger-Parker dominance (higher in 2005, P = 0.019), whereas the post-seismic community showed higher diversity (0.001 < P < 0.004) and evenness (0.001 < P < 0.023), and lower dominance (P < 0.001), in comparison with the two pre-seismic communities (Table 1).
Species-abundance analysis. The copepod community before the earthquake was best fitted by a lognormal model in 1997 and by a geometric series in 2005; in the post-seismic year (2012) the geometric series and the lognormal series fitted the copepod species abundance distribution equally well (Table 2). Thus, the copepod community shifted from the 1997 rather equitable species abundance distribution (lognormal) to a distribution characterised by lower evenness (geometric series) in 2005, and finally returned to a more balanced distribution of species abundances in 2012, after the earthquake. Because the 2005 and the 2012 communities are adequately fitted by the geometric series, rank-abundance distributions were modelled using a regression approach (Fig. 1a). The slope of the 2012 line was significantly lower than the slope of the 2005 line (equality of slopes: F = 9.881, P = 0.003), which indicates that the post-seismic community was less influenced by the most dominant species.
When stygobites and non-stygobites are analysed separately, the geometric series gave the best fit in all cases ( Table 2). When rank-abundance distributions were modelled using a regression approach, no significant differences between slopes were found (1997 vs. 2005: F = 9.081, P = 0.779; 1997 vs. 2012: F = 3.146, P = 0.091; 2005 vs. 2012: F = 1.736, P = 0.203), which indicates that non-stygobiotic species showed no changes in the dominance pattern (Fig. 1b). For the stygobiotic species (Fig. 1c), slopes were significantly different between the two pre-seismic communities (F = 10.010, P = 0.007) and marginally different between the 2005 and the post-seismic community (F = 5.149, P = 0.044), but not between the 1997 and the post-seismic community (F = 0.030, P = 0.865), which indicates that the dominance of the stygobiotic species after the earthquake returned to be similar to the 1997 situation after the increase in 2005.
Beta-diversity analysis. Average beta-diversity values were similar among years for the whole community (arithmetic mean ± SE: 0.419 ± 0.059 for 1997; 0.467 ± 0.070 for 2005; 0.444 ± 0.043 for 2012) and for the stygobiotic species (0.526 ± 0.065 for 1997; 0.632 ± 0.072 for 2005; 0.633 ± 0.049 for 2012; difference between 1997 and 2012 is in fact marginally significant), whereas increased significantly after the earthquake for the non-stygobites (0.340 ± 0.048 for 1997; 0.323 ± 0.044 for 2005; 0.535 ± 0.043 for 2012) ( and for the non-stygobiotic species, but not for the stygobites (Table 3). Beta-diversity values observed in 1997 were significantly correlated with those recorded in 2012 only for the stygobites and no correlation was detected between 2005 and 2012 beta-diversity patterns ( Table 3).

Discussion
Data about groundwater communities of unconsolidated (porous) or fractured aquifers are usually gathered by sampling drilling wells because they are directly connected with the aquifer [41][42][43][44] . However, groundwater-fed springs may represent a more interesting, albeit more complex, source of information in the case of karstic aquifers, because the structure of their animal communities is influenced by the full variety of habitats that occur across the entire aquifer and transport history 45 . Groundwater-fed springs host species belonging to different ecological categories, including crenic species (which dwell exclusively in the spring [46][47][48] but that are rare among copepods), stygobiotic species (which colonize the spring from the aquifer 49 ) and non-stygobiotic species (species that live in spring habitats being cold stenothermic or generalists, species coming from downstream, and species flushed out from surface waters of the recharge area of the aquifer 30,50,51 ). The same disturbance event may affect non-stygobiotic and stygobiotic species in different ways. The high discharge that occurred in 2005 increased the dominance (and hence decreased diversity and evenness) of non-stygobiotic species (but not of the stygobiotic ones). The earthquake that occurred in 2009 hit this already perturbed community, by inducing a significant decrease in dominance, and a significant increase in diversity and evenness of non-stygobites, which led this group of species to return to the 1997 values. Therefore, the effect of the earthquake on non-stygobiotic species was to re-equilibrate a situation altered by the high discharge in 2005. This pattern can be explained by the colonisation dynamics of these species within the TS system. Non-stygobiotic species occurring in benthic habitats of the TS system (which is an aquitard, with the carbonate bedrock with several fractures and strong upwelling in the sediment matrix overlying the bedrock 30 ) may reach the spring both via surface-water dispersal from downstream and via the aquifer when drifting from the surface recharge area 52,53 . Therefore spring colonisation by non-stygobites is not strictly dependent on groundwater dynamics 52 but is primarily regulated by the sediment texture of the spring-bed characteristics 49 . Thus, the earthquake might have induced variations in diversity, evenness and dominance of non-stygobites by re-shaping the sediment texture of TS through the effect of an increased discharge, more than by a direct effect of the aquifer dynamics. However, re-arrangements of sediment texture could have been equally produced by the high discharge that occurred in 2005. Thus, the non-stygobiotic species were affected by both the 2005 anomalous discharge and the mainshock-induced high discharge. By contrast, stygobiotic species, being only affected by changes in the karst groundwater flow that is focused to spring outlets, were much more sensitive to the effects of the earthquake than to the anomalous 2005 discharge. In fact, the 2005 anomalous discharge did not change diversity, dominance and evenness of stygobites, probably because it was within the range of the hydrological changes to which these groundwater-dweller species are used to. Post-seismic values in diversity and evenness of stygobites were higher not only in comparison with the values recorded in 2005 (and characterised by a high dominance effect), but also in comparison with 1997 values. Differently from non-stygobiotic species, stygobites reflected more directly the ecological processes that occurred into the aquifer feeding the spring system 49,52,53 . The 2009 mainshock markedly changed the Gran Sasso groundwater flow 34,54,55 as well as water isotope and chemical composition 37 , and these changes were mirrored by a variation in the stygobiotic assemblage composition 30,32 that was observed also in this study. In fact, the increase in stygobiotic diversity recorded in 2012 occurred as a result of a reduction in the abundance of most species 30 and was associated with an increase in niche overlap due to the redistribution of animals caused by the earthquake-triggered discharge 32 . This scenario is paralleled by variations in species abundance distribution patterns. For the whole community, the lognormal series was identified as the best fit model in the 1997 community, whereas the 2005 community followed the geometric series, and the 2012 community was equally best fitted by both models. The geometric model predicts very uneven abundances, broken stick predicts very even abundances, while lognormal is intermediate  and assumes a small number of very rare species 56 . Low productivity systems were claimed to have uneven species abundance distributions and be well fitted by a geometric series, while high productivity systems are well fitted by lognormal curves and exhibit the highest evenness 57 . Species abundance distributions was observed to change along a successional gradient in deciduous forest plots in Illinois, USA, with more lognormal, more even communities occurring late in succession just as for productivity, whereas early stages of succession (e.g. following the felling of trees for timber) followed a geometric series 58 . Since these pioneer observations, lognormal shapes have mainly been associated with "fully censused" communities 59 regulated by a large number of biotic and abiotic factors that together produce a lognormal abundance distribution according to the central limit theorem of statistics 40,[60][61][62] .
In the copepods of TS, the 2005 above-average discharge conditions induced a shift in the species abundance distribution from relatively even abundances (expressed by the lognormal model) to uneven abundances (expressed by the geometric series). The effect of the earthquake was to level the strongest differences in species abundances and hence to re-establish a species abundance distribution similar to that recorded in 1997, especially to the detriment of stygobiotic species. In other words, the earthquake decimated the whole community, but the impact was particularly severe on the species-abundance distribution of stygobites, hence allowing an increase in the diversity and evenness of the non-stygobiotic species, thus re-establishing a more balanced species abundance distribution in the whole community. Non-stygobiotic species showed similar patterns of abundance distribution in all three years, whereas stygobiotic species showed a significant increase in the slope of the geometric series between 1997 and 2005 and then a reduction after the earthquake. Thus, changes in the species abundance distribution were mostly driven by the stygobiotic component.
Stygobites and non-stygobites showed important differences also for changes in beta-diversity patterns between years. For the whole community and for the non-stygobiotic species, beta-diversity patterns remained similar between the two pre-seismic years, but changed between pre-seismic and post-seismic years. Thus, non-stygobites were not affected by the higher discharge in 2005, but the earthquake disrupted the previous patterns by increasing beta-diversity. In contrast, for the stygobiotic species, beta-diversity pattern of 2005 was different from those of 1997 and 2012, whereas 2012 beta-diversity pattern was similar to that of 1997. Thus, the effect of the earthquake for the stygobites was that of reconstructing a beta-diversity pattern similar to that observed in 1997. It is counter-intuitive that the earthquake determined an increase in diversity and evenness. However, this unexpected "positive" effect can be explained in consideration of the type of changes in the community structure determined by the seismic event. The earthquake changed species composition (three species disappeared) 30 and lowered species abundances that were increased by the 2005 high discharge; this induced a strong population decline of the species that dominated copepod communities, thus allowing a more equitable species-abundance distribution, and hence higher diversity and evenness.
Evidence that natural disasters may severely threaten biodiversity typically refer to population decline due to destruction of resources 63,64 whereas effects on communities may be much more complex and, in certain ecosystems, periodic disaster events may be necessary for maintaining or introducing variability in community structure [65][66][67] .
For example, communities of small mammals inhabiting areas hit by earthquakes, fires, clearcutting, and floods have low relative abundances and high species diversity [68][69][70][71] . Also, it has been observed that a reduction in food availability immediately after the disturbance might be a reason for low relative abundances 69,72 . In the case of spring copepods, it is difficult to invoke a reduction in food availability as a main reason for population decline, because particulate organic matter (POM), which constitutes a consistent food supply for copepods along with bacteria 73 , increased with the earthquake 30 . Rather, decline in species abundances after the earthquake was mainly due to a strong increase in hydraulic conductivity and to the consequent aquifer dewatering, which massively flushed out individuals of fracture-dwelling species 30,32 . Moreover, one of the disruptive consequences of the earthquake was that of redistributing the pre-seismic stygobiotic species and causing new co-occurrence patterns and interspecific interactions 32 . Although the post-seismic community showed structural parameters similar to those of the 1997, species abundances at the level of individual spring of the TS system and their distribution across springs were altered by the earthquake more than by the increased discharge in 2005.
The community dynamics in disturbed areas is highly dependent on the re-colonisation processes from source populations 25,74,75 . However, copepod re-colonisation did not occur at the TS system for the duration of our study. Changes in the aquifer structure due to the earthquake led to a change in the pre-seismic species patterns. These patterns could not be re-established by colonisation from source pools for two reasons. The first reason is that most stygobites present in the storage subsystems of the aquifer were flushed out during the mainshock; this led to a strong reduction of the populations living in the primary habitat, thus preventing subsequent recolonisation of the springs 32 . The second reason is related to the fact that discharge remained above average throughout 2012. During a high discharge event, a karst aquifer feeding a spring works like a hydraulic system under pressure 76 . The pressure in the main drain pushes the groundwater into the less permeable parts of the aquifer, thus producing what is called a "piston effect", which consists in the propagation of the hydraulic pressure at large distances 76 . In this aquifer type, which includes large conductive systems with high flow rate and current velocity, most copepods P < 0.0001. Panel c: comparison for stygobites. Regression statistics for the 1997 pre-seismic community: log(abundance) = (−0.299 ± 0.023) × rank + (2.654 ± 0.127), R² = 0.961, live in the small fractures of annex capacitive subsystems, and, even if good swimmers, they remain somewhat confined to this habitat 30,32 .
Because of this piston effect, stygobiotic copepods that survived fracture cleaning during the dewatering phase were transported to, and remained trapped in, the less permeable part of the saturated zone 30,77,78 . Thus, the 2012 post-seismic community was only marginally influenced (if any) by recolonisation processes, and changes in its structure can be substantially attributed only to the effects of the earthquake.

Conclusions
Contrary to expectations, the mainshock of L' Aquila earthquake on 6 April 2009 did not impact negatively on structural parameters of the copepod community, but re-established a more balanced species abundance distribution after the changes induced by the anomalous discharge occurred in 2005. This apparently paradoxical situation is a consequence of the different processes that characterised the stygobiotic and non-stygobiotic species, and highlights the importance of considering species ecology to understand the effects of a catastrophic event, especially when it hits a community comprised of species that differ markedly in their response to environmental changes. However, sorting groundwater taxa into ecological categories is not an easy task and involves a remarkable taxonomic effort, which complicates the study of groundwater ecosystems, when compared with surface-water ecosystems. For example, although identification to species may be important to study the response of Ephemeroptera-Plecoptera-Tricoptera (EPT) to floods 79 , it seems that, in general, identification to the genus level may be sufficient for defining ecological categories in most surface-water invertebrates 80 . By contrast, in the groundwater fauna, the same genus may include both stygobiotic and non-stygobiotic species, which have completely different adaptations, trophic roles and colonisation dynamics 81 . We are aware that ecological studies requiring taxonomic identifications to the species level are onerous, time consuming and can be performed only by trained people; yet our study demonstrates that it is necessary not to jump to misleading conclusions.
Groundwater communities are well known for their low resilience 82 and disturbance events that negatively affect their populations may easily lead to local extinction 30,82 . Most groundwater species are phylogenetic and/ or distributional relicts, thus they are species of high conservation concern 83 . Groundwater habitats are generally considered stable environments, but our study demonstrates that, in fact, they can be suddenly modified by natural changes and that copepod communities can be subject to profound alterations due to occasional, but strong disturbance events represented by earthquakes. Groundwater environments are under a variety of severe anthropogenic pressures, such as pollution and water extraction 43 . Thus, anthropogenic disturbances occurring in a community already stressed by an earthquake might have extremely negative consequences. Of course, we cannot avoid earthquakes, but we should address any effort to avoid, or at least to reduce, the impact of anthropogenic stressors.

Study area and sampling procedures. The TS area is a spring complex at the boundary of the Gran Sasso
Aquifer (GSA) located in the Gran Sasso Massif in central Italy (Apennines mountain range), featuring the highest peak south of the Alps (Corno Grande, 2922 m a.s.l.) and characterised by a high-to moderate-altitude montane landscape with low human impact. The GSA is a karstic aquifer with fast-flowing sections (karstic conduits) and interconnected low-flowing water small chambers 30 . The TS is the largest GSA-fed spring system, receiving ~65% of the GSA discharge 35 . The short (~15 km) Tirino River originating from TS joins the Aterno-Pescara River before eventually emptying into the Adriatic Sea.
Mean annual discharge at TS was relatively low in the first sampling year (1997: mean ± SD: 5.68 ± 0.21 m 3 s −1 ); it was above-average in the second sampling year (2005: 6.02 ± 0.26 m 3 s −1 ) and was well above average in the third sampling year (2012: 7.14 ± 0.26 m 3 s −1 ) due to a 3-yr rising in discharge caused by the 6.3-M w 2009 earthquake, before slowly returning to pre-seismic discharge values in summer 2013 30 .
Copepods were collected at eight sampling sites at the TS system adopting a random sampling method with four temporal replicates of three samples of 20 L in each of the eight sampled sites, for a total of 96 samples (and hence 1920 L) in each year. Subsurface samples of water were collected from the springbed (sediment patches and karstic fractures) with a hand-made Bou-Rouch pump 32 and mobile pipes hammered at each sampling site. For each replicate, a standardised sample size of 20 L was withdrawn, a volume of water/sediments that is sufficient to obtain reliable estimates of abundance of rare species 84 . The meiofauna was extracted by filtering the 20-L samples through a hand net (mesh size = 60 µm). Samples were preserved in 80% ethyl alcohol. Individuals were later counted and identified to species level. Species were assigned to two ecological categories: stygobites and non-stygobites. Stygobites (obligate groundwater species) are strictly dependent on groundwater to complete their life cycles, and are drifted or washed periodically to the aquifer outlets following the groundwater flow.  Non-stygobites found in the springs are freshwater species that live on the springbed surface, or in sediment interstices (e.g., to avoid predation) or are habitat generalists. Some of them are drifted from surface waters of the recharge area of the aquifer, others can be defined as crenobionts, i.e. species that complete their life cycle in the stable and relatively cold thermal regime of surface spring waters; other are generalist species that can colonise the spring sediments from downstream via the surface hydrological continuum. Further details about the study site, discharge patterns and sampling procedures are given elsewhere 30,36,37 . Primary data used in this study have been published in previous papers 30,32 . Statistical analysis. Because no single diversity index encompasses all the characteristics of an ideal index 85 , a combination of them that reflects richness, dominance, evenness, and relative abundance was used. Thus, the following community parameters were calculated to compare pre-and post-seismic communities: where n i is the abundance of species i and n is the overall abundance (total number of individuals); H′ ranges from 0 (one species dominates the community completely) to high values for communities with many species, each with few individuals.

H S /
where H′ is the Shannon index and S is the total number of species. This index varies from 0 (highest dominance by a single species) to 1 (all species have the same abundance). Ninety-five percent confidence intervals for all these indices were computed with a bootstrap procedure with 9999 randomizations. To compare diversity indices of pre-and post-seismic communities, we generated 9999 random matrices with two columns (samples), each with the same row and column totals as in the original data matrix. The probability of obtaining the observed difference by random sampling from a unique parental population was calculated as the number of times that the absolute difference of the indices of a replicate pair exceeded or equalled that of the original samples. Calculations were done using PAST v. 3.0 89 .
We also investigated if the earthquake modified the species abundance distributions (SADs) because the study of SADs allows inferences about patterns in the commonness and rarity of species in a community beyond those that flow from many simple diversity indices and can therefore provide insights into the effects of disturbance on ecological communities 90 . We modelled SADs using rank-abundance curves 40,85 . In the abundance-rank representations, all the species in a community are ranked from the most to the least abundant. Each species has a rank, which is plotted on the horizontal axis, while its abundance is plotted on the vertical axis: the abundance for the most abundant species is plotted first, then the next most common and so on until the array is completed by the rarest species.
Several a priori established distributions can be used to model empirical rank-abundance curves 59 . We compared pre-and post-seismic SADs using a selection of widely applied SAD models that are appropriate for discrete distributions 56,91 : the geometric series, the broken stick model, the lognormal series, and the Zipf model.
In the geometric series, also known as the niche preemption model, each species takes a constant fraction (α) of the remaining resources and the expected abundance of a species at rank r is: The only estimated parameter is the preemption coefficient α, which gives the decay rate of abundance per rank, whereas J is the total number of individuals. The geometric series is the mathematical model used to express the niche preemption hypothesis, in which the sizes of the niche hypervolumes (measured by species relative abundances) are sequentially preempted by the most abundant to the least abundant species. If in the rank-abundance plot a log scale is used for abundance, the species exactly fall along a straight line, according to the equation: where a is the species abundance, r is the respective rank, and b 0 (the intercept) and b 1 (the slope) are optimised fitting parameters 39 . With this approach, it is possible to use the regression slope to compare different species assemblages that follow the same rank-abundance distribution 39 . Among all proposed SAD models, the geometric series represents the least equitable distribution and it is known to provide a good fit to simple communities characterised by the high dominance of a few species 40,85,92 . On the opposite, most equitable empirical distributions should be modelled by the broken stick (BS) model 93 . The BS model is theoretically questionable and communities rarely are correctly characterised by such model 88,94 . Yet, the BS model is useful in comparative analyses because it represents a simple benchmark in opposition to the geometric series.
In the broken stick model, the expected abundance of species at rank r is: where J is the total number of individuals and S is the total number of species in the community 95 . In the BS there are no fitted parameters. Note that another species abundance distribution model widely used in community ecology for communities dominated by few species is the log-series [96][97][98] . However, the geometric series and the log-series abundance distributions are interrelated and are two representations of, essentially, the same underlying abundance distribution 97,99 .
The lognormal is one of the most commonly used models for describing SADs 100 . It has been derived as a null form of the distribution resulting from the central limit theorem 97 , and it is classified among the purely statistical models 56 , but can be the limit of population dynamics 101 , or niche partitioning 102,103 .
The lognormal model assumes that the logarithmic abundances are distributed normally: r where µ and σ are, respectively, the mean and standard deviation of the variable's natural logarithm, and Φ is a normal deviate. The Zipf distribution (which is a type of power law probability distribution based on branching processes 56 ) is: where p 1 is the fitted proportion of the most abundant species, and γ is a decay coefficient. Following current best practices in the study of species abundance distributions 98,104 , we used maximum likelihood estimation to fit models [105][106][107] and likelihood-based model selection to compare models 108 . The lognormal and Zipf models were fitted using generalized linear models with logarithmic link function. The preemption model was fitted as a non-linear (quasi Newton) algorithm. Since species abundances were expressed as count data, we used the Poisson error. We used the Akaike Information Criterion (AIC) to compare the fits of the different models 98,104,108 . All models were fitted and compared using the R package 'vegan' version 2.4-3 109 .
For communities that followed the geometric series, we also used the regression approach described above and tested the equality of slope with an ANCOVA approach using R 110 . We conducted all the analyses for all species and for stygobites and non-stygobites separately for the whole TS system.
Finally, we investigated how beta-diversity varied between years. To express between-spring beta-diversity we used the Morisita index, which is suggested as the most appropriate for quantitative data 111 . We tested for differences in average beta-diversity values between years using paired t-tests. Then, we correlated matrices of between-spring beta-diversity values using Mantel tests (Pearson correlation coefficient, 10000 permutations) to assess if patterns were similar between years. Calculations were done with the R package 'vegan' version 2.4-3 109 .